nx_planesweep.nx source
↩ module page · 195 lines · 8892 B
1// nx_planesweep.nx -- Reconstructs depth maps from two stereo views using plane sweep and photo-consistency scoring.
2// nx_planesweep.nx -- DENSE stereo by PLANE SWEEP: given two posed views, recover a DEPTH per pixel by
3// photo-consistency. This is the classical route to the dense surface that the cadtwin census records as
4// P1's remaining half. That half was assumed float-autograd-bound because the 2025-26 SOTA reaches it with
5// differentiable splatting -- but dense MVS does NOT require autograd: sweeping a depth hypothesis and
6// scoring patch agreement is pure search, and the census's own reframe notes classical COLMAP MVS
7// (0.532mm DTU) is accuracy-competitive with neural VGGT. All-integer, so the depth map is reproducible.
8// license_tier: ORIGINAL
9import "nx_recon3d.nx"
10
11const PS_PATCH: i64 = 4 // patch half-width -> 9x9 support; a 5x5 patch spanned only ~2 texture
12 // blocks and under-constrained the match
13const PS_BAD: i64 = 999999999 // score for an unusable hypothesis (projects outside the other view)
14
15// deterministic high-frequency texture on the world surface. Matching is only meaningful where the surface
16// actually carries detail, so the synthetic scene must be textured, not flat colour.
17// âš BLOCK-QUANTISED, and non-negative. Hashing raw world coordinates gives texture that changes every world
18// unit -- far finer than a pixel -- which ALIASES: the patch sees noise rather than structure and the match
19// score stops being meaningful. Blocks must be a few pixels across at the working depth.
20// âš MULTI-SCALE, and non-negative. Two failures live here. (a) Hashing raw world coordinates changes every
21// world unit -- far finer than a pixel -- which ALIASES: the patch sees noise, not structure. (b) A SINGLE
22// block scale gives only ~2 blocks per patch, which under-constrains the match: many depth hypotheses score
23// almost the same and the winner is close to arbitrary (measured: 230 units of error against a 54-unit
24// sweep step). A fine octave for per-pixel detail plus a coarse octave for REGIONAL uniqueness fixes both.
25func ps_tex(x: i64, y: i64) -> i64 {
26 var bx: i64 = x / 16
27 var by: i64 = y / 16
28 if x < 0 { bx = bx - 1 }
29 if y < 0 { by = by - 1 }
30 var h: i64 = (bx * 7919 + by * 104729) % 1021
31 if h < 0 { h = 0 - h }
32 var cx: i64 = x / 72
33 var cy: i64 = y / 72
34 if x < 0 { cx = cx - 1 }
35 if y < 0 { cy = cy - 1 }
36 var g: i64 = (cx * 40503 + cy * 65537) % 733
37 if g < 0 { g = 0 - g }
38 return (h % 150) + (g % 90)
39}
40// sample an image with bounds check; -1 = outside
41func ps_at(img: *i64, w: i64, h: i64, px: i64, py: i64) -> i64 {
42 if px < 0 { return 0 - 1 }
43 if py < 0 { return 0 - 1 }
44 if px >= w { return 0 - 1 }
45 if py >= h { return 0 - 1 }
46 return img[py * w + px]
47}
48// render the synthetic fronto-parallel textured plane at world z = zp, as seen by one camera
49func ps_render(cx: i64, cy: i64, cz: i64, basis: *i64, f: i64, zp: i64, w: i64, h: i64, img: *i64) -> i64 {
50 let d3: *i64 = sys_mmap(32) as *i64
51 var py: i64 = 0
52 while py < h {
53 var px: i64 = 0
54 while px < w {
55 let u: i64 = px - w / 2
56 let v: i64 = py - h / 2
57 r3_ray(basis, f, u, v, d3)
58 var val: i64 = 0
59 if d3[2] != 0 {
60 let s: i64 = ((zp - cz) * R3_Q) / d3[2]
61 if s > 0 {
62 let wx: i64 = cx + (d3[0] * s) / R3_Q
63 let wy: i64 = cy + (d3[1] * s) / R3_Q
64 val = ps_tex(wx, wy)
65 }
66 }
67 img[py * w + px] = val
68 px = px + 1
69 }
70 py = py + 1
71 }
72 return 0
73}
74// sum of absolute differences over a patch: A centred at (ax,ay), B centred at (bx,by). -1 = unusable.
75func ps_sad(a: *i64, b: *i64, w: i64, h: i64, ax: i64, ay: i64, bx: i64, by: i64) -> i64 {
76 var s: i64 = 0
77 var dy: i64 = 0 - PS_PATCH
78 while dy <= PS_PATCH {
79 var dx: i64 = 0 - PS_PATCH
80 while dx <= PS_PATCH {
81 let va: i64 = ps_at(a, w, h, ax + dx, ay + dy)
82 let vb: i64 = ps_at(b, w, h, bx + dx, by + dy)
83 if va < 0 { return 0 - 1 }
84 if vb < 0 { return 0 - 1 }
85 var d: i64 = va - vb
86 if d < 0 { d = 0 - d }
87 s = s + d
88 dx = dx + 1
89 }
90 dy = dy + 1
91 }
92 return s
93}
94// MEAN-NORMALISED patch difference: subtract each patch's own mean before differencing, which makes the
95// score invariant to a constant brightness offset between the two views.
96// ⚠REAL CAMERAS DO NOT SHARE AN EXPOSURE. Plain SAD compares absolute intensities, so two photographs of the
97// same surface taken at different exposures score as a mismatch everywhere -- the matcher does not degrade,
98// it fails. This is the difference between "works on rendered images" and "works on photographs".
99func ps_sad_norm(a: *i64, b: *i64, w: i64, h: i64, ax: i64, ay: i64, bx: i64, by: i64) -> i64 {
100 var sa: i64 = 0
101 var sb: i64 = 0
102 var cnt: i64 = 0
103 var dy: i64 = 0 - PS_PATCH
104 while dy <= PS_PATCH {
105 var dx: i64 = 0 - PS_PATCH
106 while dx <= PS_PATCH {
107 let va: i64 = ps_at(a, w, h, ax + dx, ay + dy)
108 let vb: i64 = ps_at(b, w, h, bx + dx, by + dy)
109 if va < 0 { return 0 - 1 }
110 if vb < 0 { return 0 - 1 }
111 sa = sa + va
112 sb = sb + vb
113 cnt = cnt + 1
114 dx = dx + 1
115 }
116 dy = dy + 1
117 }
118 if cnt == 0 { return 0 - 1 }
119 let ma: i64 = sa / cnt
120 let mb: i64 = sb / cnt
121 var s: i64 = 0
122 dy = 0 - PS_PATCH
123 while dy <= PS_PATCH {
124 var dx2: i64 = 0 - PS_PATCH
125 while dx2 <= PS_PATCH {
126 let va: i64 = ps_at(a, w, h, ax + dx2, ay + dy) - ma
127 let vb: i64 = ps_at(b, w, h, bx + dx2, by + dy) - mb
128 var d: i64 = va - vb
129 if d < 0 { d = 0 - d }
130 s = s + d
131 dx2 = dx2 + 1
132 }
133 dy = dy + 1
134 }
135 return s
136}
137// THE SWEEP. For every pixel of view A, walk depth hypotheses s0..s1, project the hypothesised world point
138// into view B, and keep the depth whose patch agrees best. Writes the winning s per pixel into depth[].
139// Pixels whose best hypothesis is unusable (always projecting outside B) are written as -1 -- a MISS is a
140// miss, never a silently plausible depth.
141func ps_sweep(ca: *i64, ba: *i64, cb: *i64, bb: *i64, f: i64, ia: *i64, ib: *i64, w: i64, h: i64, s0: i64, s1: i64, ns: i64, depth: *i64) -> i64 {
142 // ⚠existing call sites pass 13 args; nx_cc does NOT check arity, so the norm-aware form is a SEPARATE
143 // function and this stays a wrapper rather than growing a 14th parameter under old callers.
144 return ps_sweep_n(ca, ba, cb, bb, f, ia, ib, w, h, s0, s1, ns, depth, 0)
145}
146func ps_sweep_n(ca: *i64, ba: *i64, cb: *i64, bb: *i64, f: i64, ia: *i64, ib: *i64, w: i64, h: i64, s0: i64, s1: i64, ns: i64, depth: *i64, norm: i64) -> i64 {
147 let d3: *i64 = sys_mmap(32) as *i64
148 let o2: *i64 = sys_mmap(32) as *i64
149 var found: i64 = 0
150 var py: i64 = 0
151 while py < h {
152 var px: i64 = 0
153 while px < w {
154 let u: i64 = px - w / 2
155 let v: i64 = py - h / 2
156 r3_ray(ba, f, u, v, d3)
157 var best: i64 = PS_BAD
158 var bests: i64 = 0 - 1
159 var k: i64 = 0
160 while k < ns {
161 let s: i64 = s0 + ((s1 - s0) * k) / (ns - 1)
162 let wx: i64 = ca[0] + (d3[0] * s) / R3_Q
163 let wy: i64 = ca[1] + (d3[1] * s) / R3_Q
164 let wz: i64 = ca[2] + (d3[2] * s) / R3_Q
165 let zc: i64 = r3_project(cb[0], cb[1], cb[2], bb, f, wx, wy, wz, o2)
166 if zc > 0 {
167 let bx: i64 = o2[0] + w / 2
168 let by: i64 = o2[1] + h / 2
169 var sc: i64 = ps_sad(ia, ib, w, h, px, py, bx, by)
170 if norm == 1 { sc = ps_sad_norm(ia, ib, w, h, px, py, bx, by) }
171 if sc >= 0 {
172 if sc < best {
173 best = sc
174 bests = s
175 }
176 }
177 }
178 k = k + 1
179 }
180 depth[py * w + px] = bests
181 if bests >= 0 { found = found + 1 }
182 px = px + 1
183 }
184 py = py + 1
185 }
186 return found
187}
188// ground-truth ray distance to the plane z = zp for pixel (px,py) of a camera -- the exact answer the sweep
189// is trying to find, so the benchmark needs no dataset
190func ps_true_s(cx: i64, cy: i64, cz: i64, basis: *i64, f: i64, zp: i64, w: i64, h: i64, px: i64, py: i64) -> i64 {
191 let d3: *i64 = sys_mmap(32) as *i64
192 r3_ray(basis, f, px - w / 2, py - h / 2, d3)
193 if d3[2] == 0 { return 0 - 1 }
194 return ((zp - cz) * R3_Q) / d3[2]
195}