nx_loop_subdiv.nx source
↩ module page · 147 lines · 5694 B
1// nx_loop_subdiv.nx -- LOOP SUBDIVISION (Charles Loop, 1987): the triangle-mesh SMOOTH-LIMIT subdivision scheme.
2// Unlike linear midpoint refinement (nx_subdiv sd_subdivide, which keeps original vertices FIXED = same silhouette,
3// only more triangles), Loop REPOSITIONS even (original) vertices toward the neighbour centroid, and weights each
4// new odd (edge) vertex by the two OPPOSITE vertices (3/8(a+b) + 1/8(c+d)) -> the mesh converges to a smooth limit
5// surface (a coarse octahedron -> a round blob). Integer fx256, deterministic. This is the census GAP
6// "no Catmull-Clark smooth-limit subdivision" answered for triangle meshes. license_tier: ORIGINAL
7import "nx_subdiv.nx"
8
9// distinct neighbours of vertex vi in meshIn -> valence returned; neighbour position SUM into sum[3].
10// seen[] = caller scratch (>= nverts flags, all 0 on entry); left all-0 on exit (reusable).
11func lsd_neighbours(meshIn: i64, vi: i64, sum: *i64, seen: *i64) -> i64 {
12 let h: *i64 = m3_hdr(meshIn)
13 var n: i64 = 0
14 sum[0] = 0; sum[1] = 0; sum[2] = 0
15 var t: i64 = 0
16 while t < h[1] {
17 let tr: *i64 = m3_tri(meshIn, t)
18 var has: i64 = 0
19 if tr[0] == vi { has = 1 }
20 if tr[1] == vi { has = 1 }
21 if tr[2] == vi { has = 1 }
22 if has == 1 {
23 var k: i64 = 0
24 while k < 3 {
25 let nb: i64 = tr[k]
26 if nb != vi {
27 if seen[nb] == 0 {
28 seen[nb] = 1
29 let vp: *i64 = m3_vert(meshIn, nb)
30 sum[0] = sum[0] + vp[0]; sum[1] = sum[1] + vp[1]; sum[2] = sum[2] + vp[2]
31 n = n + 1
32 }
33 }
34 k = k + 1
35 }
36 }
37 t = t + 1
38 }
39 // clear the flags we set (re-scan the incident tris) so seen[] is reusable
40 t = 0
41 while t < h[1] {
42 let tr: *i64 = m3_tri(meshIn, t)
43 var has: i64 = 0
44 if tr[0] == vi { has = 1 }
45 if tr[1] == vi { has = 1 }
46 if tr[2] == vi { has = 1 }
47 if has == 1 {
48 var k: i64 = 0
49 while k < 3 { seen[tr[k]] = 0; k = k + 1 }
50 }
51 t = t + 1
52 }
53 return n
54}
55
56// opposite vertex of edge (a,b) in the triangle that is NOT curT; returns vertex id, or -1 for a boundary edge.
57func lsd_opposite(meshIn: i64, a: i64, b: i64, curT: i64) -> i64 {
58 let h: *i64 = m3_hdr(meshIn)
59 var t: i64 = 0
60 while t < h[1] {
61 if t != curT {
62 let tr: *i64 = m3_tri(meshIn, t)
63 var ha: i64 = 0
64 var hb: i64 = 0
65 if tr[0] == a { ha = 1 }
66 if tr[1] == a { ha = 1 }
67 if tr[2] == a { ha = 1 }
68 if tr[0] == b { hb = 1 }
69 if tr[1] == b { hb = 1 }
70 if tr[2] == b { hb = 1 }
71 if ha == 1 {
72 if hb == 1 {
73 var opp: i64 = 0 - 1
74 var k: i64 = 0
75 while k < 3 { if tr[k] != a { if tr[k] != b { opp = tr[k] } } k = k + 1 }
76 return opp
77 }
78 }
79 }
80 t = t + 1
81 }
82 return 0 - 1
83}
84
85// odd (edge) vertex for edge (a,b) of triangle curT whose third vertex is c. Dedup by exact position (sd_getvert):
86// both triangles sharing edge (a,b) compute the IDENTICAL position, so the edge vertex is shared (watertight).
87func lsd_oddvert(meshIn: i64, meshOut: i64, a: i64, b: i64, c: i64, curT: i64) -> i64 {
88 let pa: *i64 = m3_vert(meshIn, a)
89 let pb: *i64 = m3_vert(meshIn, b)
90 let d: i64 = lsd_opposite(meshIn, a, b, curT)
91 var ox: i64 = (pa[0] + pb[0]) / 2
92 var oy: i64 = (pa[1] + pb[1]) / 2
93 var oz: i64 = (pa[2] + pb[2]) / 2
94 if d >= 0 {
95 let pc: *i64 = m3_vert(meshIn, c)
96 let pd: *i64 = m3_vert(meshIn, d)
97 ox = (3 * (pa[0] + pb[0]) + (pc[0] + pd[0])) / 8
98 oy = (3 * (pa[1] + pb[1]) + (pc[1] + pd[1])) / 8
99 oz = (3 * (pa[2] + pb[2]) + (pc[2] + pd[2])) / 8
100 }
101 return sd_getvert(meshOut, ox, oy, oz)
102}
103
104// one Loop subdivision level: meshIn -> meshOut. seen/sum = caller scratch.
105func loop_subdivide(meshIn: i64, meshOut: i64, seen: *i64, sum: *i64) -> i64 {
106 m3_init(meshOut)
107 let h: *i64 = m3_hdr(meshIn)
108 let nv: i64 = h[0]
109 // Step 1: repositioned EVEN vertices (original indices 0..nv-1 preserved in meshOut)
110 var i: i64 = 0
111 while i < nv {
112 let n: i64 = lsd_neighbours(meshIn, i, sum, seen)
113 let v: *i64 = m3_vert(meshIn, i)
114 var ox: i64 = v[0]
115 var oy: i64 = v[1]
116 var oz: i64 = v[2]
117 if n == 3 {
118 ox = v[0] + 3 * (sum[0] - 3 * v[0]) / 16
119 oy = v[1] + 3 * (sum[1] - 3 * v[1]) / 16
120 oz = v[2] + 3 * (sum[2] - 3 * v[2]) / 16
121 }
122 if n >= 4 {
123 ox = v[0] + 3 * (sum[0] - n * v[0]) / (8 * n)
124 oy = v[1] + 3 * (sum[1] - n * v[1]) / (8 * n)
125 oz = v[2] + 3 * (sum[2] - n * v[2]) / (8 * n)
126 }
127 m3_add_vert(meshOut, ox, oy, oz)
128 i = i + 1
129 }
130 // Steps 2 + 3: odd (edge) vertices (position-deduped) + retriangulate each face into 4
131 var t: i64 = 0
132 while t < h[1] {
133 let tr: *i64 = m3_tri(meshIn, t)
134 let a: i64 = tr[0]
135 let b: i64 = tr[1]
136 let c: i64 = tr[2]
137 let iab: i64 = lsd_oddvert(meshIn, meshOut, a, b, c, t)
138 let ibc: i64 = lsd_oddvert(meshIn, meshOut, b, c, a, t)
139 let ica: i64 = lsd_oddvert(meshIn, meshOut, c, a, b, t)
140 m3_add_tri(meshOut, a, iab, ica)
141 m3_add_tri(meshOut, iab, b, ibc)
142 m3_add_tri(meshOut, ica, ibc, c)
143 m3_add_tri(meshOut, iab, ibc, ica)
144 t = t + 1
145 }
146 return h[1] * 4
147}