nx_visual_hull.nx source
↩ module page · 168 lines · 7349 B
1// nx_visual_hull.nx -- SILHOUETTE -> 3D FORM. Builds a signed field by intersecting the back-
2// projections of N silhouettes (the visual hull), which surface_nets then polygonizes into a mesh.
3//
4// This is the missing rung between a photograph and geometry. The ecosystem could decode JPEG
5// (nx_jpeg_decode_rgb), polygonize a field (nx_isosurf), render a mesh (nx_trimesh) and import OBJ
6// (nx_obj_import) -- but had no way to turn images INTO a field. Verified absent before writing:
7// every `carve` in the corpus is terrain/geology.
8//
9// ★VIEWS MUST DIFFER IN ANGLE. An anterior and a posterior photograph project along the SAME axis
10// and therefore yield the SAME silhouette outline (mirrored) -- that is ONE shape constraint, not
11// two, however many images you have. Front + LATERAL is what constrains a form. Anterior/posterior
12// pairs are still valuable, but for per-surface ALBEDO, not for shape.
13//
14// ★CSG: the visual hull is an INTERSECTION of back-projected silhouette prisms, and intersection of
15// signed distance fields is max(). Each view is folded in with max, so views can be added
16// incrementally and in any order.
17//
18// Caller supplies cos/sin of the view yaw in Q10 (1024 = 1.0) so this organ needs no trig table and
19// stays exact for the axis-aligned cases that matter (front 1024/0, lateral 0/1024).
20// license_tier: ORIGINAL
21import "nx_syscalls.nx"
22
23const VH_BIG: i64 = 1000000
24const VH_Q: i64 = 1024
25
26func vh_abs(v: i64) -> i64 { if v < 0 { return 0 - v } return v }
27func vh_min(a: i64, b: i64) -> i64 { if a < b { return a } return b }
28func vh_max(a: i64, b: i64) -> i64 { if a > b { return a } return b }
29
30// Background segmentation for a STUDIO plate: the backdrop colour is estimated from the four
31// corners (median-of-4 per channel via mid-two average), and a pixel is FOREGROUND when its
32// Manhattan distance from that colour exceeds `tol`. mask[] gets 1 = subject, 0 = background.
33func vh_bg_mask(rgb: *u8, w: i64, h: i64, tol: i64, mask: *u8) -> i64 {
34 if w <= 0 { return 0 - 1 }
35 if h <= 0 { return 0 - 1 }
36 let cx: *i64 = sys_mmap(4 * 8) as *i64
37 let cy: *i64 = sys_mmap(4 * 8) as *i64
38 cx[0] = 0; cy[0] = 0
39 cx[1] = w - 1; cy[1] = 0
40 cx[2] = 0; cy[2] = h - 1
41 cx[3] = w - 1; cy[3] = h - 1
42 var br: i64 = 0
43 var bg: i64 = 0
44 var bb: i64 = 0
45 var i: i64 = 0
46 while i < 4 {
47 let o: i64 = (cy[i] * w + cx[i]) * 3
48 br = br + (rgb[o] as i64)
49 bg = bg + (rgb[o+1] as i64)
50 bb = bb + (rgb[o+2] as i64)
51 i = i + 1
52 }
53 br = br / 4; bg = bg / 4; bb = bb / 4
54 var y: i64 = 0
55 var fg: i64 = 0
56 while y < h {
57 var x: i64 = 0
58 while x < w {
59 let o: i64 = (y * w + x) * 3
60 let d: i64 = vh_abs((rgb[o] as i64) - br) + vh_abs((rgb[o+1] as i64) - bg) + vh_abs((rgb[o+2] as i64) - bb)
61 if d > tol { mask[y*w + x] = 1 as u8; fg = fg + 1 } else { mask[y*w + x] = 0 as u8 }
62 x = x + 1
63 }
64 y = y + 1
65 }
66 sys_munmap(cx as *u8, 32)
67 sys_munmap(cy as *u8, 32)
68 return fg
69}
70
71// Signed 2D distance from a binary mask by two-pass chamfer (3-4 metric, scaled by 3 so the
72// diagonal step is an integer). NEGATIVE inside the subject, POSITIVE outside -- the SDF convention
73// surface_nets expects. out[] is w*h.
74func vh_sdf2d(mask: *u8, w: i64, h: i64, out: *i64) -> i64 {
75 let din: *i64 = sys_mmap(w * h * 8) as *i64
76 let dout: *i64 = sys_mmap(w * h * 8) as *i64
77 var i: i64 = 0
78 while i < w * h {
79 if mask[i] == (1 as u8) { din[i] = VH_BIG; dout[i] = 0 } else { din[i] = 0; dout[i] = VH_BIG }
80 i = i + 1
81 }
82 // forward then backward chamfer over BOTH fields
83 var y: i64 = 0
84 while y < h {
85 var x: i64 = 0
86 while x < w {
87 let k: i64 = y * w + x
88 if x > 0 { din[k] = vh_min(din[k], din[k-1] + 3); dout[k] = vh_min(dout[k], dout[k-1] + 3) }
89 if y > 0 { din[k] = vh_min(din[k], din[k-w] + 3); dout[k] = vh_min(dout[k], dout[k-w] + 3) }
90 if y > 0 { if x > 0 { din[k] = vh_min(din[k], din[k-w-1] + 4); dout[k] = vh_min(dout[k], dout[k-w-1] + 4) } }
91 if y > 0 { if x < w-1 { din[k] = vh_min(din[k], din[k-w+1] + 4); dout[k] = vh_min(dout[k], dout[k-w+1] + 4) } }
92 x = x + 1
93 }
94 y = y + 1
95 }
96 y = h - 1
97 while y >= 0 {
98 var x: i64 = w - 1
99 while x >= 0 {
100 let k: i64 = y * w + x
101 if x < w-1 { din[k] = vh_min(din[k], din[k+1] + 3); dout[k] = vh_min(dout[k], dout[k+1] + 3) }
102 if y < h-1 { din[k] = vh_min(din[k], din[k+w] + 3); dout[k] = vh_min(dout[k], dout[k+w] + 3) }
103 if y < h-1 { if x < w-1 { din[k] = vh_min(din[k], din[k+w+1] + 4); dout[k] = vh_min(dout[k], dout[k+w+1] + 4) } }
104 if y < h-1 { if x > 0 { din[k] = vh_min(din[k], din[k+w-1] + 4); dout[k] = vh_min(dout[k], dout[k+w-1] + 4) } }
105 x = x - 1
106 }
107 y = y - 1
108 }
109 i = 0
110 while i < w * h {
111 // inside -> negative. /3 converts chamfer units back to pixels.
112 if mask[i] == (1 as u8) { out[i] = 0 - (din[i] / 3) } else { out[i] = dout[i] / 3 }
113 i = i + 1
114 }
115 sys_munmap(din as *u8, w * h * 8)
116 sys_munmap(dout as *u8, w * h * 8)
117 return 0
118}
119
120// sample a 2D field with clamping
121func vh_samp(f: *i64, w: i64, h: i64, x: i64, y: i64) -> i64 {
122 var xx: i64 = x
123 var yy: i64 = y
124 if xx < 0 { xx = 0 }
125 if yy < 0 { yy = 0 }
126 if xx > w - 1 { xx = w - 1 }
127 if yy > h - 1 { yy = h - 1 }
128 return f[yy * w + xx]
129}
130
131// Fold ONE view into the 3D field by back-projection. grid dims are (GX+1)(GY+1)(GZ+1) corners,
132// origin (ox,oy,oz), spacing `cell`. The view looks along a yaw given by (cyaw,syaw) in Q10.
133// `first` = 1 initialises the grid; 0 intersects (max) into it.
134// px_per_unit is Q10: image pixels per world unit. (icx,icy) is the image-space centre of the
135// subject, (wcx,wcy) the world-space centre it corresponds to.
136func vh_carve_view(f2d: *i64, w: i64, h: i64, cyaw: i64, syaw: i64,
137 grid: *i64, GX: i64, GY: i64, GZ: i64,
138 ox: i64, oy: i64, oz: i64, cell: i64,
139 icx: i64, icy: i64, wcx: i64, wcy: i64, px_per_unit: i64, first: i64) -> i64 {
140 if px_per_unit <= 0 { return 0 - 1 }
141 let GY1: i64 = GY + 1
142 let GZ1: i64 = GZ + 1
143 var i: i64 = 0
144 while i <= GX {
145 let wx: i64 = ox + i * cell
146 var j: i64 = 0
147 while j <= GY {
148 let wy: i64 = oy + j * cell
149 // vertical maps straight through for any yaw about the Y axis
150 let py: i64 = icy + ((wy - wcy) * px_per_unit) / VH_Q
151 var k: i64 = 0
152 while k <= GZ {
153 let wz: i64 = oz + k * cell
154 // horizontal image axis = the world axis perpendicular to the view direction
155 let u: i64 = (wx * cyaw + wz * syaw) / VH_Q
156 let px: i64 = icx + ((u - wcx) * px_per_unit) / VH_Q
157 // image y grows DOWN, world y grows UP
158 let d: i64 = vh_samp(f2d, w, h, px, (h - 1) - (py - 0))
159 let idx: i64 = (i * GY1 + j) * GZ1 + k
160 if first == 1 { grid[idx] = d } else { grid[idx] = vh_max(grid[idx], d) }
161 k = k + 1
162 }
163 j = j + 1
164 }
165 i = i + 1
166 }
167 return 0
168}