`nx_csg.nx` source

↩ module page · 245 lines · 10833 B
1// nx_csg.nx -- CSG solid -> watertight mesh via marching tetrahedra over an SDF.
2//
3// Samples a two-primitive boolean field (nx_sdf) on a grid, splits every cube into
4// 6 tetrahedra (Kuhn subdivision along the 0-7 diagonal -> shared faces split
5// identically, no cracks), and emits the zero-isosurface.  Edge crossings are
6// computed in a CANONICAL endpoint order so a grid edge shared by two tets yields
7// the IDENTICAL point (watertight, the same fix the slicer needed).  Triangle
8// winding is forced outward via an inside-reference test.  Output is a pipeline-
9// native Q14-mm NxMesh -> flows to nx_stl_write_mesh -> slicer, like every other
10// solid.  Robust by construction: SDF booleans never fail the way B-rep does.
11// license_tier: ORIGINAL
12
13import "nx_syscalls.nx"
14import "nx_mesh.nx"
15import "nx_sdf.nx"
16
17// canonical edge crossing: order endpoints by (z,y,x) so both tets sharing this
18// grid edge interpolate the SAME point.  out[0..2] = crossing (Q14).
19func mt_cross(out: *i64, P: *i64, i: i64, j: i64, D: *i64) -> i64 {
  var xi: i64 = P[i*3+0]; var yi: i64 = P[i*3+1]; var zi: i64 = P[i*3+2]; var di: i64 = D[i]
  var xj: i64 = P[j*3+0]; var yj: i64 = P[j*3+1]; var zj: i64 = P[j*3+2]; var dj: i64 = D[j]
  var swap: i64 = 0
  if zj < zi { swap = 1 }
  if zj == zi { if yj < yi { swap = 1 } }
  if zj == zi { if yj == yi { if xj < xi { swap = 1 } } }
  if swap == 1 {
      let tx: i64 = xi; let ty: i64 = yi; let tz: i64 = zi; let td: i64 = di
      xi = xj; yi = yj; zi = zj; di = dj
      xj = tx; yj = ty; zj = tz; dj = td
  }
  let den: i64 = di - dj
  out[0] = xi + (xj - xi) * di / den
  out[1] = yi + (yj - yi) * di / den
  out[2] = zi + (zj - zi) * di / den
  return 0
36}
37
38// append a triangle, flipping winding so its normal points AWAY from (rx,ry,rz)
39// (a point inside the solid) -> outward-facing.
40func mt_emit(tribuf: *i64, ntri: *i64, cap: i64, p0: *i64, p1: *i64, p2: *i64,
          rx: i64, ry: i64, rz: i64) -> i64 {
  if ntri[0] >= cap { return 0 }
  let ux: i64 = p1[0]-p0[0]; let uy: i64 = p1[1]-p0[1]; let uz: i64 = p1[2]-p0[2]
  let vx: i64 = p2[0]-p0[0]; let vy: i64 = p2[1]-p0[1]; let vz: i64 = p2[2]-p0[2]
  let nx2: i64 = uy*vz - uz*vy
  let ny:  i64 = uz*vx - ux*vz
  let nz:  i64 = ux*vy - uy*vx
  let cxx: i64 = (p0[0]+p1[0]+p2[0])/3
  let cyy: i64 = (p0[1]+p1[1]+p2[1])/3
  let czz: i64 = (p0[2]+p1[2]+p2[2])/3
  let dot: i64 = nx2*(cxx-rx) + ny*(cyy-ry) + nz*(czz-rz)
  let base: i64 = ntri[0]*9
  tribuf[base+0]=p0[0]; tribuf[base+1]=p0[1]; tribuf[base+2]=p0[2]
  if dot >= 0 {
      tribuf[base+3]=p1[0]; tribuf[base+4]=p1[1]; tribuf[base+5]=p1[2]
      tribuf[base+6]=p2[0]; tribuf[base+7]=p2[1]; tribuf[base+8]=p2[2]
  } else {
      tribuf[base+3]=p2[0]; tribuf[base+4]=p2[1]; tribuf[base+5]=p2[2]
      tribuf[base+6]=p1[0]; tribuf[base+7]=p1[1]; tribuf[base+8]=p1[2]
  }
  ntri[0] = ntri[0] + 1
  return 0
63}
64
65// march one tetrahedron (P = 4 corners flat [12], D = 4 field values).  ins/outs/
66// c0..c3 are caller-owned scratch (avoids per-tet mmap storms).
67func mt_tet(tribuf: *i64, ntri: *i64, cap: i64, P: *i64, D: *i64,
         ins: *i64, outs: *i64, c0: *i64, c1: *i64, c2: *i64, c3: *i64) -> i64 {
  var ni: i64 = 0; var no: i64 = 0
  var i: i64 = 0
  while i < 4 {
      if D[i] < 0 { ins[ni] = i; ni = ni + 1 } else { outs[no] = i; no = no + 1 }
      i = i + 1
  }
  if ni == 0 { return 0 }
  if ni == 4 { return 0 }
  if ni == 1 {
      let a: i64 = ins[0]
      mt_cross(c0, P, a, outs[0], D); mt_cross(c1, P, a, outs[1], D); mt_cross(c2, P, a, outs[2], D)
      mt_emit(tribuf, ntri, cap, c0, c1, c2, P[a*3+0], P[a*3+1], P[a*3+2])
  }
  if ni == 3 {
      let d: i64 = outs[0]
      mt_cross(c0, P, d, ins[0], D); mt_cross(c1, P, d, ins[1], D); mt_cross(c2, P, d, ins[2], D)
      let rx: i64 = (P[ins[0]*3+0]+P[ins[1]*3+0]+P[ins[2]*3+0])/3
      let ry: i64 = (P[ins[0]*3+1]+P[ins[1]*3+1]+P[ins[2]*3+1])/3
      let rz: i64 = (P[ins[0]*3+2]+P[ins[1]*3+2]+P[ins[2]*3+2])/3
      mt_emit(tribuf, ntri, cap, c0, c1, c2, rx, ry, rz)
  }
  if ni == 2 {
      let a: i64 = ins[0]; let b: i64 = ins[1]; let cc: i64 = outs[0]; let dd: i64 = outs[1]
      mt_cross(c0, P, a, cc, D)   // a-c
      mt_cross(c1, P, b, cc, D)   // b-c
      mt_cross(c2, P, b, dd, D)   // b-d
      mt_cross(c3, P, a, dd, D)   // a-d
      let rx: i64 = (P[a*3+0]+P[b*3+0])/2
      let ry: i64 = (P[a*3+1]+P[b*3+1])/2
      let rz: i64 = (P[a*3+2]+P[b*3+2])/2
      mt_emit(tribuf, ntri, cap, c0, c1, c2, rx, ry, rz)   // ac, bc, bd
      mt_emit(tribuf, ntri, cap, c0, c2, c3, rx, ry, rz)   // ac, bd, ad
  }
  return 0
103}
104
105// Extract the boolean (a op b) solid over the box [lo,hi] at grid resolution res.
106// k_q14 > 0 rounds the junction (fillet); k_q14 = 0 is a sharp boolean.
107func nx_csg_extract(a: *NxSdfPrim, b: *NxSdfPrim, op: i64,
                  lox: i64, loy: i64, loz: i64,
                  hix: i64, hiy: i64, hiz: i64, res: i64, k_q14: i64) -> *NxMesh {
  let np: i64 = res + 1
  let field: *i64 = (sys_mmap(np * np * np * 8)) as *i64
  var k: i64 = 0
  while k <= res {
      let z: i64 = loz + (hiz - loz) * k / res
      var j: i64 = 0
      while j <= res {
          let y: i64 = loy + (hiy - loy) * j / res
          var ii: i64 = 0
          while ii <= res {
              let x: i64 = lox + (hix - lox) * ii / res
              field[ii + j*np + k*np*np] = nx_csg_field_k(a, b, op, k_q14, x, y, z)
              ii = ii + 1
          }
          j = j + 1
      }
      k = k + 1
  }
128
  // 6-tet (Kuhn) decomposition along cube diagonal 0-7
  let tets: *i64 = (sys_mmap(24 * 8)) as *i64
  tets[0]=0;tets[1]=7;tets[2]=1;tets[3]=3;     tets[4]=0;tets[5]=7;tets[6]=3;tets[7]=2
  tets[8]=0;tets[9]=7;tets[10]=2;tets[11]=6;   tets[12]=0;tets[13]=7;tets[14]=6;tets[15]=4
  tets[16]=0;tets[17]=7;tets[18]=4;tets[19]=5; tets[20]=0;tets[21]=7;tets[22]=5;tets[23]=1
134
  let cap: i64 = 48 * res * res + 4096
  let tribuf: *i64 = (sys_mmap(cap * 9 * 8)) as *i64
  let ntri: *i64 = (sys_mmap(8)) as *i64
  ntri[0] = 0
139
  // scratch (allocate once)
  let cp: *i64 = (sys_mmap(8 * 3 * 8)) as *i64    // 8 corner positions
  let cd: *i64 = (sys_mmap(8 * 8)) as *i64        // 8 corner field values
  let P:  *i64 = (sys_mmap(12 * 8)) as *i64
  let D:  *i64 = (sys_mmap(4 * 8)) as *i64
  let ins:  *i64 = (sys_mmap(4 * 8)) as *i64
  let outs: *i64 = (sys_mmap(4 * 8)) as *i64
  let q0: *i64 = (sys_mmap(24)) as *i64
  let q1: *i64 = (sys_mmap(24)) as *i64
  let q2: *i64 = (sys_mmap(24)) as *i64
  let q3: *i64 = (sys_mmap(24)) as *i64
151
  var ck: i64 = 0
  while ck < res {
      var cj: i64 = 0
      while cj < res {
          var ci: i64 = 0
          while ci < res {
              // gather the cube's 8 corners
              var c: i64 = 0
              while c < 8 {
                  let gi: i64 = ci + (c & 1)
                  let gj: i64 = cj + ((c >> 1) & 1)
                  let gk: i64 = ck + ((c >> 2) & 1)
                  cp[c*3+0] = lox + (hix - lox) * gi / res
                  cp[c*3+1] = loy + (hiy - loy) * gj / res
                  cp[c*3+2] = loz + (hiz - loz) * gk / res
                  cd[c] = field[gi + gj*np + gk*np*np]
                  c = c + 1
              }
              // march the 6 tets
              var t: i64 = 0
              while t < 6 {
                  var v: i64 = 0
                  while v < 4 {
                      let corner: i64 = tets[t*4+v]
                      P[v*3+0] = cp[corner*3+0]
                      P[v*3+1] = cp[corner*3+1]
                      P[v*3+2] = cp[corner*3+2]
                      D[v] = cd[corner]
                      v = v + 1
                  }
                  mt_tet(tribuf, ntri, cap, P, D, ins, outs, q0, q1, q2, q3)
                  t = t + 1
              }
              ci = ci + 1
          }
          cj = cj + 1
      }
      ck = ck + 1
  }
191
  let nt: i64 = ntri[0]
  if nt <= 0 { return 0 as *NxMesh }
  let m: *NxMesh = nx_mesh_alloc(nt * 3, nt, 0)
  if (m as i64) == 0 { return m }
  var t2: i64 = 0
  while t2 < nt {
      let base: i64 = t2 * 9
      nx_mesh_set_vertex(m, t2*3+0, tribuf[base+0], tribuf[base+1], tribuf[base+2], 0)
      nx_mesh_set_vertex(m, t2*3+1, tribuf[base+3], tribuf[base+4], tribuf[base+5], 0)
      nx_mesh_set_vertex(m, t2*3+2, tribuf[base+6], tribuf[base+7], tribuf[base+8], 0)
      nx_mesh_set_triangle(m, t2, t2*3+0, t2*3+1, t2*3+2)
      t2 = t2 + 1
  }
  return m
206}
207
208// convenience: box (half hx,hy,hz at origin) MINUS a Z-cylinder (radius r) through it.
209func nx_csg_box_minus_cyl(hx: i64, hy: i64, hz: i64, r: i64, res: i64) -> *NxMesh {
  let box: *NxSdfPrim = nx_sdf_make(NX_SDF_BOX, 0, 0, 0, hx, hy, hz)
  let cyl: *NxSdfPrim = nx_sdf_make(NX_SDF_CYL, 0, 0, 0, r, hz + 16384, 0)   // taller than box
  let pad: i64 = 16384   // 1mm padding so the surface isn't clipped at the grid edge
  return nx_csg_extract(box, cyl, NX_CSG_DIFF,
                        0 - hx - pad, 0 - hy - pad, 0 - hz - pad,
                        hx + pad, hy + pad, hz + pad, res, 0)
216}
217
218// convenience: a rounded box of overall half-extent H with corner radius r.
219// (Face centres reach +/-H; edges/corners are rounded by r.)
220func nx_csg_roundbox_mesh(H: i64, r: i64, res: i64) -> *NxMesh {
  let rb: *NxSdfPrim = nx_sdf_make4(NX_SDF_ROUNDBOX, 0, 0, 0, H - r, H - r, H - r, r)
  let e: i64 = H + 16384
  return nx_csg_extract(rb, rb, NX_CSG_UNION, 0 - e, 0 - e, 0 - e, e, e, e, res, 0)
224}
225
226// convenience: a single torus (major R, minor r) meshed via marching tetrahedra.
227// (union with itself is idempotent, so we reuse the two-primitive extractor.)
228func nx_csg_torus_mesh(R: i64, r: i64, res: i64) -> *NxMesh {
  let t: *NxSdfPrim = nx_sdf_make(NX_SDF_TORUS, 0, 0, 0, R, r, 0)
  let ext: i64 = R + r + 16384
  let zext: i64 = r + 16384
  return nx_csg_extract(t, t, NX_CSG_UNION, 0 - ext, 0 - ext, 0 - zext, ext, ext, zext, res, 0)
233}
234
235// convenience: smooth (filleted) union of two spheres of radius r, centres at
236// x = +/- sep on the X axis, blended over fillet radius k_q14.  Demonstrates the
237// SDF fillet -- the spheres merge with a smooth neck no B-rep edge-surgery needed.
238func nx_csg_two_sphere_fillet(r: i64, sep: i64, k: i64, res: i64) -> *NxMesh {
  let sa: *NxSdfPrim = nx_sdf_make(NX_SDF_SPHERE, 0 - sep, 0, 0, r, 0, 0)
  let sb: *NxSdfPrim = nx_sdf_make(NX_SDF_SPHERE, sep, 0, 0, r, 0, 0)
  let pad: i64 = k + 16384
  let ex: i64 = sep + r + pad
  let ey: i64 = r + pad
  return nx_csg_extract(sa, sb, NX_CSG_UNION, 0 - ex, 0 - ey, 0 - ey, ex, ey, ey, res, k)
245}
nx_csg.nx source

`nx_csg.nx` source
