1/**
2 * shape -> SurfaceModel.
3 *
4 * Two things happen per B-rep face:
5 * - **Triangulation** (always): OpenCASCADE's BRepMesh produces a triangle
6 * approximation that respects the face's trimming. This is what three.js
7 * draws, and it never depends on NURBS extraction succeeding.
8 * - **NURBS extraction** (best-effort): after converting the shape's faces to
9 * B-spline form, we read each face's poles / weights / knots / degree — the
10 * actual polynomial that defines the surface — plus a few sampled isocurves.
11 *
12 * The mesh density is controlled by a dimensionless `quality` in [0, 1] fed to
13 * BRepMesh in *relative* mode, so it is independent of the model's scale.
14 */
15import type { NurbsSurface, Patch, SurfaceModel, TriMesh } from '../model/types'
16import type { OpenCascade, Shape } from './types'
18export interface BuiltModel {
19 model: SurfaceModel
20 /** The (NURBS-converted, if possible) shape kept for re-tessellation. */
21 meshShape: Shape
22}
24const EMPTY_TRI: TriMesh = {
25 positions: new Float32Array(0),
26 normals: new Float32Array(0),
27 indices: new Uint32Array(0),
28}
30/** Map quality in [0,1] to BRepMesh relative-linear and angular deflections. */
31function deflections(quality: number): { lin: number; ang: number } {
32 const q = Math.min(1, Math.max(0, quality))
33 return {
34 lin: 0.02 - (0.02 - 0.0007) * q, // relative to edge size
35 ang: 0.8 - (0.8 - 0.15) * q, // radians
36 }
37}
39function runMesh(oc: OpenCascade, shape: Shape, quality: number): void {
40 const { lin, ang } = deflections(quality)
41 try {
42 // second arg (forceFaceDeflection) is required in this OCCT build; Clean
43 // wipes the old triangulation so a coarser deflection actually coarsens.
44 oc.BRepTools.Clean(shape, false)
45 } catch {
46 /* Clean unavailable — IncrementalMesh still recomputes when deflection differs */
47 }
48 new oc.BRepMesh_IncrementalMesh_2(shape, lin, true, ang, false)
49}
51/** Per-vertex normals from the triangle connectivity (fallback if OCCT's fail). */
52function computeNormals(positions: Float32Array, indices: Uint32Array): Float32Array {
53 const normals = new Float32Array(positions.length)
54 for (let t = 0; t + 2 < indices.length; t += 3) {
55 const a = indices[t] * 3
56 const b = indices[t + 1] * 3
57 const c = indices[t + 2] * 3
58 const ux = positions[b] - positions[a]
59 const uy = positions[b + 1] - positions[a + 1]
60 const uz = positions[b + 2] - positions[a + 2]
61 const vx = positions[c] - positions[a]
62 const vy = positions[c + 1] - positions[a + 1]
63 const vz = positions[c + 2] - positions[a + 2]
64 const nx = uy * vz - uz * vy
65 const ny = uz * vx - ux * vz
66 const nz = ux * vy - uy * vx
67 for (const i of [a, b, c]) {
68 normals[i] += nx
69 normals[i + 1] += ny
70 normals[i + 2] += nz
71 }
72 }
73 for (let i = 0; i < normals.length; i += 3) {
74 const len = Math.hypot(normals[i], normals[i + 1], normals[i + 2]) || 1
75 normals[i] /= len
76 normals[i + 1] /= len
77 normals[i + 2] /= len
78 }
79 return normals
80}
82/** Read the BRepMesh triangulation stored on a face (world coordinates). */
83function readFaceTriangulation(oc: OpenCascade, face: Shape): TriMesh {
84 const loc = new oc.TopLoc_Location_1()
85 const handle = oc.BRep_Tool.Triangulation(face, loc, 0)
86 if (handle.IsNull()) {
87 loc.delete()
88 return EMPTY_TRI
89 }
90 const tri = handle.get()
91 const trsf = loc.Transformation()
92 const nNodes = tri.NbNodes()
94 const positions = new Float32Array(nNodes * 3)
95 for (let i = 1; i <= nNodes; i++) {
96 const node = tri.Node(i)
97 const p = node.Transformed(trsf)
98 positions[(i - 1) * 3] = p.X()
99 positions[(i - 1) * 3 + 1] = p.Y()
100 positions[(i - 1) * 3 + 2] = p.Z()
101 node.delete()
102 p.delete()
103 }
105 const forward = face.Orientation_1() === oc.TopAbs_Orientation.TopAbs_FORWARD
106 const triangles = tri.Triangles()
107 const nTri = tri.NbTriangles()
108 const indices = new Uint32Array(nTri * 3)
109 for (let nt = 1; nt <= nTri; nt++) {
110 const t = triangles.Value(nt)
111 let n1 = t.Value(1)
112 let n2 = t.Value(2)
113 const n3 = t.Value(3)
114 if (!forward) {
115 const tmp = n1
116 n1 = n2
117 n2 = tmp
118 }
119 indices[(nt - 1) * 3] = n1 - 1
120 indices[(nt - 1) * 3 + 1] = n2 - 1
121 indices[(nt - 1) * 3 + 2] = n3 - 1
122 t.delete()
123 }
125 let normals: Float32Array = new Float32Array(nNodes * 3)
126 let haveNormals = false
127 try {
128 const pc = new oc.Poly_Connect_2(handle)
129 const nrm = new oc.TColgp_Array1OfDir_2(1, nNodes)
130 oc.StdPrs_ToolTriangulatedShape.Normal(face, pc, nrm)
131 for (let i = nrm.Lower(); i <= nrm.Upper(); i++) {
132 const d0 = nrm.Value(i)
133 const d = d0.Transformed(trsf)
134 const s = forward ? 1 : -1
135 normals[(i - 1) * 3] = s * d.X()
136 normals[(i - 1) * 3 + 1] = s * d.Y()
137 normals[(i - 1) * 3 + 2] = s * d.Z()
138 d0.delete()
139 d.delete()
140 }
141 nrm.delete()
142 pc.delete()
143 haveNormals = true
144 } catch {
145 /* fall back below */
146 }
147 if (!haveNormals) normals = computeNormals(positions, indices)
149 triangles.delete()
150 trsf.delete()
151 handle.delete()
152 loc.delete()
153 return { positions, normals, indices }
154}
156/** Try to convert every face of the shape to B-spline form. */
157function nurbsConvert(oc: OpenCascade, shape: Shape): { shape: Shape; ok: boolean } {
158 try {
159 const conv = new oc.BRepBuilderAPI_NurbsConvert_2(shape, true)
160 return { shape: conv.Shape(), ok: true }
161 } catch {
162 return { shape, ok: false }
163 }
164}
166/** Extract a face's B-spline surface data plus sampled isocurves, or null. */
167function extractNurbs(oc: OpenCascade, face: Shape): { nurbs: NurbsSurface; isoLines: Float32Array[] } | null {
168 let handle: Shape | null = null
169 let bsHandle: Shape | null = null
170 try {
171 handle = oc.BRep_Tool.Surface_2(face)
172 const raw = handle.get()
173 if (raw.$$?.ptrType?.name !== 'Geom_BSplineSurface*') return null
175 bsHandle = new oc.Handle_Geom_BSplineSurface_2(raw)
176 const bs = bsHandle.get()
178 const uDegree: number = bs.UDegree()
179 const vDegree: number = bs.VDegree()
180 const nu: number = bs.NbUPoles()
181 const nv: number = bs.NbVPoles()
182 const rational = bs.IsURational() || bs.IsVRational()
184 const poles = new Float32Array(nu * nv * 3)
185 const weights = rational ? new Float32Array(nu * nv) : null
186 for (let i = 1; i <= nu; i++) {
187 for (let j = 1; j <= nv; j++) {
188 const p = bs.Pole(i, j)
189 const idx = ((i - 1) * nv + (j - 1)) * 3
190 poles[idx] = p.X()
191 poles[idx + 1] = p.Y()
192 poles[idx + 2] = p.Z()
193 p.delete()
194 if (weights) weights[(i - 1) * nv + (j - 1)] = bs.Weight(i, j)
195 }
196 }
198 const uKnots: number[] = []
199 const uMults: number[] = []
200 for (let i = 1; i <= bs.NbUKnots(); i++) {
201 uKnots.push(bs.UKnot(i))
202 uMults.push(bs.UMultiplicity(i))
203 }
204 const vKnots: number[] = []
205 const vMults: number[] = []
206 for (let i = 1; i <= bs.NbVKnots(); i++) {
207 vKnots.push(bs.VKnot(i))
208 vMults.push(bs.VMultiplicity(i))
209 }
211 const isoLines = sampleIsoLines(bs, uKnots, vKnots)
212 const nurbs: NurbsSurface = { uDegree, vDegree, nu, nv, poles, weights, uKnots, uMults, vKnots, vMults }
213 return { nurbs, isoLines }
214 } catch {
215 return null
216 } finally {
217 bsHandle?.delete()
218 handle?.delete()
219 }
220}
222/** Sample constant-u and constant-v curves on the surface for the isocurve view. */
223function sampleIsoLines(bs: Shape, uKnots: number[], vKnots: number[]): Float32Array[] {
224 const NLINES = 5
225 const NS = 40
226 const uMin = uKnots[0]
227 const uMax = uKnots[uKnots.length - 1]
228 const vMin = vKnots[0]
229 const vMax = vKnots[vKnots.length - 1]
230 const lines: Float32Array[] = []
231 const evalTo = (line: Float32Array, s: number, u: number, v: number) => {
232 const p = bs.Value(u, v)
233 line[s * 3] = p.X()
234 line[s * 3 + 1] = p.Y()
235 line[s * 3 + 2] = p.Z()
236 p.delete()
237 }
238 for (let a = 0; a < NLINES; a++) {
239 const u = uMin + ((uMax - uMin) * a) / (NLINES - 1)
240 const line = new Float32Array(NS * 3)
241 for (let s = 0; s < NS; s++) evalTo(line, s, u, vMin + ((vMax - vMin) * s) / (NS - 1))
242 lines.push(line)
243 }
244 for (let a = 0; a < NLINES; a++) {
245 const v = vMin + ((vMax - vMin) * a) / (NLINES - 1)
246 const line = new Float32Array(NS * 3)
247 for (let s = 0; s < NS; s++) evalTo(line, s, uMin + ((uMax - uMin) * s) / (NS - 1), v)
248 lines.push(line)
249 }
250 return lines
251}
253/** Iterate the faces of a shape, calling `fn` with each `TopoDS_Face`. */
254function forEachFace(oc: OpenCascade, shape: Shape, fn: (face: Shape, index: number) => void): void {
255 const exp = new oc.TopExp_Explorer_1()
256 let index = 0
257 for (
258 exp.Init(shape, oc.TopAbs_ShapeEnum.TopAbs_FACE, oc.TopAbs_ShapeEnum.TopAbs_SHAPE);
259 exp.More();
260 exp.Next()
261 ) {
262 const face = oc.TopoDS.Face_1(exp.Current())
263 fn(face, index++)
264 face.delete()
265 }
266 exp.delete()
267}
269/** Build a full SurfaceModel from a freshly created / imported shape. */
270export function buildModel(
271 oc: OpenCascade,
272 rawShape: Shape,
273 quality: number,
274 source: SurfaceModel['source'],
275 raw?: SurfaceModel['raw'],
276): BuiltModel {
277 const { shape: meshShape, ok } = nurbsConvert(oc, rawShape)
278 runMesh(oc, meshShape, quality)
280 const patches: Patch[] = []
281 forEachFace(oc, meshShape, (face, index) => {
282 const tri = readFaceTriangulation(oc, face)
283 const extracted = ok ? extractNurbs(oc, face) : null
284 patches.push({
285 kind: 'nurbs',
286 id: index,
287 tri,
288 nurbs: extracted?.nurbs ?? null,
289 isoLines: extracted?.isoLines,
290 })
291 })
293 return { model: { patches, source, raw }, meshShape }
294}
296/** Re-tessellate an existing model at a new quality, reusing its NURBS data. */
297export function retessellate(
298 oc: OpenCascade,
299 meshShape: Shape,
300 quality: number,
301 patches: Patch[],
302): Patch[] {
303 runMesh(oc, meshShape, quality)
304 const tris: TriMesh[] = []
305 forEachFace(oc, meshShape, (face) => {
306 tris.push(readFaceTriangulation(oc, face))
307 })
308 return patches.map((p, i) => ({ ...p, tri: tris[i] ?? p.tri }))
309}