nx_isosurf.nx source
↩ module page · 145 lines · 7105 B
1// nx_isosurf.nx -- ★F0 UNIFY implicit + explicit (from the graphics foundation census): a SURFACE-NETS polygonizer that
2// turns a signed-distance FIELD into a real TRIANGLE MESH. This is the bridge that lets us MESH THE BEING (or any SDF):
3// sample the SDF on a grid, place one vertex per surface-straddling cell (averaged edge crossings), connect the 4 cells
4// around every sign-changing grid edge into a quad. Per-vertex normals come from the SDF GRADIENT (so smooth Gouraud
5// shading is winding-independent). Fills the nx_trimesh buffers -> render with the existing rasterizer. license_tier: ORIGINAL
6import "nx_syscalls.nx"
7import "nx_trimesh.nx"
8import "nx_vecmath.nx"
9const K_MAGIC_1024: i64 = 1024
10const K_MAGIC_4096: i64 = 4096
11
12func iso_isqrt(v: i64) -> i64 { return vm_isqrt(v) }
13
14// ★polygonize a corner-SDF grid (dims (GX+1)x(GY+1)x(GZ+1)) at iso-level `iso`. Cells GX x GY x GZ. Origin ox,oy,oz;
15// corner (i,j,k) is at (ox+i*cell, oy+j*cell, oz+k*cell). grid[(i*(GY+1)+j)*(GZ+1)+k] = SDF value. Emits into trimesh.
16func surface_nets(grid: *i64, GX: i64, GY: i64, GZ: i64, ox: i64, oy: i64, oz: i64, cell: i64, iso: i64, col: i64) -> i64 {
17 let GY1: i64 = GY+1; let GZ1: i64 = GZ+1
18 let cv: *i64 = sys_mmap(GX*GY*GZ*8) as *i64
19 let cval: *i64 = sys_mmap(8*8) as *i64 // scratch, reused per cell (NOT allocated in the hot loop)
20 let ea: *i64 = sys_mmap(12*8) as *i64
21 let eb: *i64 = sys_mmap(12*8) as *i64
22 ea[0]=0;eb[0]=1; ea[1]=2;eb[1]=3; ea[2]=4;eb[2]=5; ea[3]=6;eb[3]=7
23 ea[4]=0;eb[4]=2; ea[5]=1;eb[5]=3; ea[6]=4;eb[6]=6; ea[7]=5;eb[7]=7
24 ea[8]=0;eb[8]=4; ea[9]=1;eb[9]=5; ea[10]=2;eb[10]=6; ea[11]=3;eb[11]=7
25 // ---- pass 1: one vertex per surface-straddling cell ----
26 var ci: i64 = 0
27 while ci < GX {
28 var cj: i64 = 0
29 while cj < GY {
30 var ck: i64 = 0
31 while ck < GZ {
32 cv[(ci*GY+cj)*GZ+ck] = 0-1
33 var bit: i64 = 0
34 while bit < 8 {
35 let dx: i64 = bit&1; let dy: i64 = (bit>>1)&1; let dz: i64 = (bit>>2)&1
36 cval[bit] = grid[((ci+dx)*GY1 + (cj+dy))*GZ1 + (ck+dz)]
37 bit = bit + 1
38 }
39 var sx: i64 = 0; var sy: i64 = 0; var sz: i64 = 0; var cnt: i64 = 0
40 var e: i64 = 0
41 while e < 12 {
42 let a: i64 = ea[e]; let b: i64 = eb[e]
43 var ina: i64 = 0; if cval[a] <= iso { ina = 1 } // inside = val<=iso (handles exact-0 corners)
44 var inb: i64 = 0; if cval[b] <= iso { inb = 1 }
45 if ina != inb {
46 let den: i64 = cval[b]-cval[a]
47 var t: i64 = 0
48 if den != 0 { t = (iso-cval[a])*K_MAGIC_1024/den }
49 let ax: i64 = ox+(ci+(a&1))*cell; let ay: i64 = oy+(cj+((a>>1)&1))*cell; let az: i64 = oz+(ck+((a>>2)&1))*cell
50 let bx: i64 = ox+(ci+(b&1))*cell; let by: i64 = oy+(cj+((b>>1)&1))*cell; let bz: i64 = oz+(ck+((b>>2)&1))*cell
51 sx = sx + ax + (bx-ax)*t/K_MAGIC_1024
52 sy = sy + ay + (by-ay)*t/K_MAGIC_1024
53 sz = sz + az + (bz-az)*t/K_MAGIC_1024
54 cnt = cnt + 1
55 }
56 e = e + 1
57 }
58 if cnt > 0 {
59 let vx: i64 = sx/cnt; let vy: i64 = sy/cnt; let vz: i64 = sz/cnt
60 // normal = SDF gradient at the base corner (central differences, clamped)
61 var ai: i64 = ci+1; if ai>GX {ai=GX}
62 var bi: i64 = ci-1; if bi<0 {bi=0}
63 var aj: i64 = cj+1; if aj>GY {aj=GY}
64 var bj: i64 = cj-1; if bj<0 {bj=0}
65 var ak: i64 = ck+1; if ak>GZ {ak=GZ}
66 var bk: i64 = ck-1; if bk<0 {bk=0}
67 let gx: i64 = grid[(ai*GY1+cj)*GZ1+ck] - grid[(bi*GY1+cj)*GZ1+ck]
68 let gy: i64 = grid[(ci*GY1+aj)*GZ1+ck] - grid[(ci*GY1+bj)*GZ1+ck]
69 let gz: i64 = grid[(ci*GY1+cj)*GZ1+ak] - grid[(ci*GY1+cj)*GZ1+bk]
70 let gl: i64 = iso_isqrt(gx*gx+gy*gy+gz*gz) + 1
71 let idx: i64 = tm_vert_n(vx, vy, vz, gx*K_MAGIC_4096/gl, gy*K_MAGIC_4096/gl, gz*K_MAGIC_4096/gl)
72 cv[(ci*GY+cj)*GZ+ck] = idx
73 }
74 ck = ck + 1
75 }
76 cj = cj + 1
77 }
78 ci = ci + 1
79 }
80 // ---- pass 2: quad per sign-changing INTERIOR corner-grid edge (4 surrounding cells) ----
81 // x-edges
82 var i: i64 = 0
83 while i < GX {
84 var j: i64 = 1
85 while j < GY {
86 var k: i64 = 1
87 while k < GZ {
88 var q0: i64 = 0; if grid[(i*GY1+j)*GZ1+k] <= iso { q0 = 1 }
89 var q1: i64 = 0; if grid[((i+1)*GY1+j)*GZ1+k] <= iso { q1 = 1 }
90 if q0 != q1 {
91 let m: i64 = cv[(i*GY+(j-1))*GZ+(k-1)]; let n: i64 = cv[(i*GY+j)*GZ+(k-1)]; let o: i64 = cv[(i*GY+j)*GZ+k]; let p: i64 = cv[(i*GY+(j-1))*GZ+k]
92 if m>=0 { if n>=0 { if o>=0 { if p>=0 {
93 if q0 == 1 { tm_quad(m,n,o,p,col) } else { tm_quad(m,p,o,n,col) } // sign-consistent outward winding
94 } } } }
95 }
96 k = k + 1
97 }
98 j = j + 1
99 }
100 i = i + 1
101 }
102 // y-edges
103 i = 1
104 while i < GX {
105 var j: i64 = 0
106 while j < GY {
107 var k: i64 = 1
108 while k < GZ {
109 var q0: i64 = 0; if grid[(i*GY1+j)*GZ1+k] <= iso { q0 = 1 }
110 var q1: i64 = 0; if grid[(i*GY1+(j+1))*GZ1+k] <= iso { q1 = 1 }
111 if q0 != q1 {
112 let m: i64 = cv[((i-1)*GY+j)*GZ+(k-1)]; let n: i64 = cv[(i*GY+j)*GZ+(k-1)]; let o: i64 = cv[(i*GY+j)*GZ+k]; let p: i64 = cv[((i-1)*GY+j)*GZ+k]
113 if m>=0 { if n>=0 { if o>=0 { if p>=0 {
114 if q0 == 1 { tm_quad(m,n,o,p,col) } else { tm_quad(m,p,o,n,col) } // sign-consistent outward winding
115 } } } }
116 }
117 k = k + 1
118 }
119 j = j + 1
120 }
121 i = i + 1
122 }
123 // z-edges
124 i = 1
125 while i < GX {
126 var j: i64 = 1
127 while j < GY {
128 var k: i64 = 0
129 while k < GZ {
130 var q0: i64 = 0; if grid[(i*GY1+j)*GZ1+k] <= iso { q0 = 1 }
131 var q1: i64 = 0; if grid[(i*GY1+j)*GZ1+(k+1)] <= iso { q1 = 1 }
132 if q0 != q1 {
133 let m: i64 = cv[((i-1)*GY+(j-1))*GZ+k]; let n: i64 = cv[(i*GY+(j-1))*GZ+k]; let o: i64 = cv[(i*GY+j)*GZ+k]; let p: i64 = cv[((i-1)*GY+j)*GZ+k]
134 if m>=0 { if n>=0 { if o>=0 { if p>=0 {
135 if q0 == 1 { tm_quad(m,n,o,p,col) } else { tm_quad(m,p,o,n,col) } // sign-consistent outward winding
136 } } } }
137 }
138 k = k + 1
139 }
140 j = j + 1
141 }
142 i = i + 1
143 }
144 return 0
145}