nx_step_nurbs.nx source

1// nx_step_nurbs.nx -- RATIONAL B-SPLINE SURFACE evaluation from STEP (cadtwin P3a: the ROUND parts, exact). 2// OpenCascade exports cylinders/holes as complex-instance rational B-spline surfaces, degree (1,3): linear along 3// the axis x rational-CUBIC around -- weights (1,1/3,1/3,1) = an EXACT circle arc (NURBS circle). This organ 4// parses the complex instance (control net 2x4 + weights) and evaluates the rational Bezier in fixed point: 5// Bernstein Q12, weights Q20 (fx256 would corrupt 1/3 -> Q20 keeps radius error ~ppm), coords fx256 (1 unit = 6// 3.90625 um). ENGINEERING-TWIN doctrine: dimensional truth first -- the gate measures roundness/radius/positions 7// against the file's own declared values, in integer-exact teeth where the math is exact (t=0, 1/2, 1) and 8// measured-um bounds where it rounds. Units: as1 declares SI_UNIT(.MILLI.,.METRE.) -> mm confirmed at runtime, 9// never assumed. license_tier: ORIGINAL 10import "nx_step_tess.nx" 11const SN_MAGIC_1000000000000: i64 = 1000000000000 12 13const SN_Q12: i64 = 4096 14const SN_QW: i64 = 1048576 15 16// naive substring find of nul-terminated pat in buf[off..off+len); returns ABSOLUTE index or -1 17func sgn_find(buf: *u8, off: i64, len: i64, pat: *u8) -> i64 { 18 var pl: i64 = 0 19 while pat[pl] != (0 as u8) { pl = pl + 1 } 20 if pl == 0 { return 0 - 1 } 21 var i: i64 = off 22 let e: i64 = off + len - pl 23 while i <= e { 24 var j: i64 = 0 25 var ok: i64 = 1 26 while j < pl { 27 if buf[i + j] != pat[j] { ok = 0; j = pl } else { j = j + 1 } 28 } 29 if ok == 1 { return i } 30 i = i + 1 31 } 32 return 0 - 1 33} 34// parse a non-negative int at p; out[0]=val; returns next pos 35func sgn_int(buf: *u8, p: i64, e: i64, out: *i64) -> i64 { 36 var v: i64 = p 37 var val: i64 = 0 38 var going: i64 = 1 39 while going == 1 { 40 if v < e { let d: i64 = buf[v] as i64; if sg_isdigit(d) == 1 { val = val * 10 + (d - 48); v = v + 1 } else { going = 0 } } else { going = 0 } 41 } 42 out[0] = val 43 return v 44} 45// parse a STEP real at p scaled by q (like sg_parse_real but arbitrary Q); out[0]=val; returns next pos. 46// fraction digits stop accumulating once the denominator passes 1e12 (guard vs overflow; precision far beyond Q20). 47func sgn_realq(buf: *u8, p: i64, e: i64, q: i64, out: *i64) -> i64 { 48 var v: i64 = p 49 var sign: i64 = 1 50 if v < e { if buf[v] == (45 as u8) { sign = 0 - 1; v = v + 1 } else { if buf[v] == (43 as u8) { v = v + 1 } } } 51 var intval: i64 = 0 52 var g1: i64 = 1 53 while g1 == 1 { 54 if v < e { let d: i64 = buf[v] as i64; if sg_isdigit(d) == 1 { intval = intval * 10 + (d - 48); v = v + 1 } else { g1 = 0 } } else { g1 = 0 } 55 } 56 var val: i64 = intval * q 57 if v < e { if buf[v] == (46 as u8) { 58 v = v + 1 59 var fnum: i64 = 0 60 var fden: i64 = 1 61 var g2: i64 = 1 62 while g2 == 1 { 63 if v < e { let d: i64 = buf[v] as i64; if sg_isdigit(d) == 1 { 64 if fden <= SN_MAGIC_1000000000000 { fnum = fnum * 10 + (d - 48); fden = fden * 10 } 65 v = v + 1 66 } else { g2 = 0 } } else { g2 = 0 } 67 } 68 if fden > 1 { val = val + (fnum * q) / fden } 69 } } 70 var hase: i64 = 0 71 if v < e { if buf[v] == (69 as u8) { hase = 1; v = v + 1 } } 72 if hase == 0 { if v < e { if buf[v] == (101 as u8) { hase = 1; v = v + 1 } } } 73 if hase == 1 { 74 var esign: i64 = 1 75 if v < e { if buf[v] == (45 as u8) { esign = 0 - 1; v = v + 1 } else { if buf[v] == (43 as u8) { v = v + 1 } } } 76 var exp: i64 = 0 77 var g3: i64 = 1 78 while g3 == 1 { 79 if v < e { let d: i64 = buf[v] as i64; if sg_isdigit(d) == 1 { exp = exp * 10 + (d - 48); v = v + 1 } else { g3 = 0 } } else { g3 = 0 } 80 } 81 var k: i64 = 0 82 while k < exp { if esign > 0 { val = val * 10 } else { val = val / 10 } k = k + 1 } 83 } 84 out[0] = val * sign 85 return v 86} 87// rounded signed division 88func sgn_rdiv(num: i64, den: i64) -> i64 { 89 if num >= 0 { return (num + den / 2) / den } 90 return 0 - ((0 - num + den / 2) / den) 91} 92// cubic Bernstein basis at t (Q12) -> b[4] (Q12; exact at t=0, Q12/2, Q12) 93func sgn_bern(t: i64, b: *i64) -> i64 { 94 let u: i64 = SN_Q12 - t 95 b[0] = ((u * u) / SN_Q12) * u / SN_Q12 96 b[1] = 3 * (((u * u) / SN_Q12) * t / SN_Q12) 97 b[2] = 3 * (((t * t) / SN_Q12) * u / SN_Q12) 98 b[3] = ((t * t) / SN_Q12) * t / SN_Q12 99 return 0 100} 101 102// ---- nt ctx: [0]=st [1]=mesh [2]=counters(*i64 c2[16]) [3]=ctrl(*i64 24: 2x4 pts xyz fx256) [4]=wts(*i64 8, Q20) 103// [5]=out2 [6]=bern(*i64 4) [7]=refs(*i64 16) 104// c2: 0=spline-faces 1=parsed 2=tessellated 3=skip-degree 4=skip-refs 5=skip-weights 6=tris 7=BOUNDARY-EXACT-FAILS 105 106// parse the complex-instance rational surface (entity index sidx) into nt ctrl/wts. 1 ok. 107func sgn_parse_surface(nt: *i64, sidx: i64) -> i64 { 108 let st: *i64 = nt[0] as *i64 109 let c2: *i64 = nt[2] as *i64 110 let buf: *u8 = st[0] as *u8 111 let aao: *i64 = st[5] as *i64 112 let aal: *i64 = st[6] as *i64 113 let off: i64 = aao[sidx] 114 let len: i64 = aal[sidx] 115 let out2: *i64 = nt[5] as *i64 116 var p: i64 = sgn_find(buf, off, len, "B_SPLINE_SURFACE(" as *u8) 117 if p < 0 { c2[3] = c2[3] + 1; return 0 } 118 p = p + 17 119 let e: i64 = off + len 120 p = sgn_int(buf, p, e, out2) 121 let udeg: i64 = out2[0] 122 p = p + 1 // ',' 123 p = sgn_int(buf, p, e, out2) 124 let vdeg: i64 = out2[0] 125 if udeg != 1 { c2[3] = c2[3] + 1; return 0 } 126 if vdeg != 3 { c2[3] = c2[3] + 1; return 0 } 127 // control net refs: exactly 8 for a single Bezier span (2 rows x 4) 128 let refs: *i64 = nt[7] as *i64 129 let wk: i64 = sgn_find(buf, p, e - p, "B_SPLINE_SURFACE_WITH_KNOTS" as *u8) 130 var refend: i64 = e 131 if wk >= 0 { refend = wk } 132 let nr: i64 = sgt_refs(buf, p, refend - p, refs, 9) 133 if nr != 8 { c2[4] = c2[4] + 1; return 0 } 134 let ctrl: *i64 = nt[3] as *i64 135 var i: i64 = 0 136 while i < 8 { 137 if sgt_tuple_of(st, refs[i], ((ctrl as i64) + i * 24) as *i64) == 0 { c2[4] = c2[4] + 1; return 0 } 138 i = i + 1 139 } 140 // weights: 8 reals at Q20 after RATIONAL_B_SPLINE_SURFACE( 141 var q: i64 = sgn_find(buf, off, len, "RATIONAL_B_SPLINE_SURFACE(" as *u8) 142 if q < 0 { c2[5] = c2[5] + 1; return 0 } 143 q = q + 26 144 let wts: *i64 = nt[4] as *i64 145 i = 0 146 while i < 8 { 147 q = sg_skip_to_num(buf, q, e) 148 if q >= e { c2[5] = c2[5] + 1; return 0 } 149 q = sgn_realq(buf, q, e, SN_QW, out2) 150 if out2[0] <= 0 { c2[5] = c2[5] + 1; return 0 } 151 wts[i] = out2[0] 152 i = i + 1 153 } 154 return 1 155} 156 157// evaluate rational cubic row (0/1) at t (Q12) -> out3 fx256 158func sgn_eval(nt: *i64, row: i64, t: i64, out3: *i64) -> i64 { 159 let ctrl: *i64 = nt[3] as *i64 160 let wts: *i64 = nt[4] as *i64 161 let b: *i64 = nt[6] as *i64 162 sgn_bern(t, b) 163 let base: i64 = row * 4 164 var den: i64 = 0 165 var nx: i64 = 0 166 var ny: i64 = 0 167 var nz: i64 = 0 168 var j: i64 = 0 169 while j < 4 { 170 let bw: i64 = b[j] * wts[base + j] 171 den = den + bw 172 nx = nx + bw * ctrl[(base + j) * 3] 173 ny = ny + bw * ctrl[(base + j) * 3 + 1] 174 nz = nz + bw * ctrl[(base + j) * 3 + 2] 175 j = j + 1 176 } 177 if den <= 0 { out3[0] = 0; out3[1] = 0; out3[2] = 0; return 0 } 178 out3[0] = sgn_rdiv(nx, den) 179 out3[1] = sgn_rdiv(ny, den) 180 out3[2] = sgn_rdiv(nz, den) 181 return 1 182} 183 184// tessellate one spline face (surface entity sidx) into the mesh: N=16 segments x 2 rims. BOUNDARY-EXACTNESS 185// invariant: eval(row,0)==ctrl[row][0] and eval(row,1)==ctrl[row][3] integer-exact, else counted + face skipped. 186func sgn_tess_face(nt: *i64, sidx: i64) -> i64 { 187 let c2: *i64 = nt[2] as *i64 188 if sgn_parse_surface(nt, sidx) == 0 { return 0 } 189 c2[1] = c2[1] + 1 190 let ctrl: *i64 = nt[3] as *i64 191 let mesh: i64 = nt[1] 192 let ev: *i64 = ((nt[6] as i64) + 64) as *i64 // eval scratch after bern block 193 // boundary invariant both rows 194 var row: i64 = 0 195 while row < 2 { 196 sgn_eval(nt, row, 0, ev) 197 let c0: i64 = row * 4 * 3 198 if ev[0] != ctrl[c0] { c2[7] = c2[7] + 1; return 0 } 199 if ev[1] != ctrl[c0 + 1] { c2[7] = c2[7] + 1; return 0 } 200 if ev[2] != ctrl[c0 + 2] { c2[7] = c2[7] + 1; return 0 } 201 sgn_eval(nt, row, SN_Q12, ev) 202 let c3: i64 = (row * 4 + 3) * 3 203 if ev[0] != ctrl[c3] { c2[7] = c2[7] + 1; return 0 } 204 if ev[1] != ctrl[c3 + 1] { c2[7] = c2[7] + 1; return 0 } 205 if ev[2] != ctrl[c3 + 2] { c2[7] = c2[7] + 1; return 0 } 206 row = row + 1 207 } 208 // emit 2 x 17 verts + 32 tris 209 let h: *i64 = m3_hdr(mesh) 210 let vbase: i64 = h[0] 211 row = 0 212 while row < 2 { 213 var i: i64 = 0 214 while i <= 16 { 215 sgn_eval(nt, row, (i * SN_Q12) / 16, ev) 216 m3_add_vert(mesh, ev[0], ev[1], ev[2]) 217 i = i + 1 218 } 219 row = row + 1 220 } 221 var i2: i64 = 0 222 while i2 < 16 { 223 let a: i64 = vbase + i2 224 let bb: i64 = vbase + i2 + 1 225 let cc: i64 = vbase + 17 + i2 226 let dd: i64 = vbase + 17 + i2 + 1 227 m3_add_tri(mesh, a, bb, cc) 228 m3_add_tri(mesh, bb, dd, cc) 229 c2[6] = c2[6] + 2 230 i2 = i2 + 1 231 } 232 c2[2] = c2[2] + 1 233 return 1 234} 235 236// sweep every ADVANCED_FACE whose surface is a complex-instance spline; tessellate into nt mesh 237func sgn_tessellate_all(nt: *i64) -> i64 { 238 let st: *i64 = nt[0] as *i64 239 let c2: *i64 = nt[2] as *i64 240 let anl: *i64 = st[4] as *i64 241 let out2: *i64 = nt[5] as *i64 242 let cnt: i64 = st[7] 243 var idx: i64 = 0 244 while idx < cnt { 245 if sp_name_is(st, idx, "ADVANCED_FACE" as *u8) == 1 { 246 let sid: i64 = sgt_arg_ref(st, idx, 2, out2) 247 if sid >= 0 { 248 let sidx: i64 = sp_find(st, sid) 249 if sidx >= 0 { if anl[sidx] == 0 { 250 c2[0] = c2[0] + 1 251 sgn_tess_face(nt, sidx) 252 } } 253 } 254 } 255 idx = idx + 1 256 } 257 return c2[2] 258} 259 260// runtime unit confirmation: some complex instance declares LENGTH_UNIT + SI_UNIT(.MILLI.,.METRE.) -> mm 261func sgn_unit_mm(st: *i64) -> i64 { 262 let buf: *u8 = st[0] as *u8 263 let anl: *i64 = st[4] as *i64 264 let aao: *i64 = st[5] as *i64 265 let aal: *i64 = st[6] as *i64 266 let cnt: i64 = st[7] 267 var idx: i64 = 0 268 while idx < cnt { 269 if anl[idx] == 0 { 270 let f1: i64 = sgn_find(buf, aao[idx], aal[idx], "LENGTH_UNIT()" as *u8) 271 if f1 >= 0 { 272 let f2: i64 = sgn_find(buf, aao[idx], aal[idx], "SI_UNIT(.MILLI.,.METRE.)" as *u8) 273 if f2 >= 0 { return 1 } 274 } 275 } 276 idx = idx + 1 277 } 278 return 0 279}