nx_recon3d_gate.nx source

1// nx_recon3d_gate.nx -- benchmark the classical multi-view triangulation (cadtwin P1). KAT the fixed-point 2// (projection<->ray consistency, 2-view midpoint recovers a known point), then the CHAMFER benchmark the DTU 3// SOTA uses: synthesize 4 camera views of a KNOWN 9-point cloud (a ~160-unit part), triangulate each point, 4// Chamfer vs ground truth -- CLEAN (near-exact) + NOISY (+/-3px, graceful) + NEG-CONTROL (wrong correspondences 5// => large Chamfer, so the metric discriminates). Reports Chamfer in units + permille of the object diagonal. 6// expect_exit: 0 license_tier: ORIGINAL 7import "nx_recon3d.nx" 8 9func rg_puts(s: *u8) -> i64 { var n: i64 = 0; while s[n] != (0 as u8) { n = n + 1 } sys_write(1, s, n); return 0 } 10func rg_putn(v: i64) -> i64 { 11 let bb: *u8 = sys_mmap(28) 12 var m: i64 = v 13 if m < 0 { rg_puts("-" as *u8); m = 0 - m } 14 let t: *u8 = sys_mmap(28) 15 var k: i64 = 0 16 if m == 0 { t[0] = 48 as u8; k = 1 } 17 while m > 0 { t[k] = (48 + (m % 10)) as u8; m = m / 10; k = k + 1 } 18 var i: i64 = 0 19 while i < k { bb[i] = t[k - 1 - i]; i = i + 1 } 20 sys_write(1, bb, k) 21 return 0 22} 23func rg_tooth(name: *u8, pass: i64, fails: *i64) -> i64 { 24 rg_puts(" " as *u8); rg_puts(name); rg_puts(" -> " as *u8) 25 if pass == 1 { rg_puts("PASS\n" as *u8); return 0 } 26 rg_puts("FAIL\n" as *u8) 27 fails[0] = fails[0] + 1 28 return 0 29} 30// distance from point P to ray (C, d Q14 unit) 31func rg_raydist(px: i64, py: i64, pz: i64, cx: i64, cy: i64, cz: i64, dx: i64, dy: i64, dz: i64) -> i64 { 32 let vx: i64 = px - cx 33 let vy: i64 = py - cy 34 let vz: i64 = pz - cz 35 let proj: i64 = (dx * vx + dy * vy + dz * vz) / R3_Q // coord 36 let qx: i64 = cx + (proj * dx) / R3_Q 37 let qy: i64 = cy + (proj * dy) / R3_Q 38 let qz: i64 = cz + (proj * dz) / R3_Q 39 return r3_dist(px, py, pz, qx, qy, qz) 40} 41 42func main() -> i64 { 43 rg_puts("=== nx_recon3d_gate -- classical multi-view triangulation, Chamfer benchmark (photo->3D core) ===\n" as *u8) 44 let fails: *i64 = sys_mmap(8) as *i64 45 fails[0] = 0 46 47 // ground-truth cloud: 8 box corners (+-800) + center = 9 points. object diagonal ~ 2771 units. 48 let gt: *i64 = sys_mmap(512) as *i64 49 let np: i64 = 9 50 let H: i64 = 800 51 gt[0] = 0 - H; gt[1] = 0 - H; gt[2] = 0 - H 52 gt[3] = H; gt[4] = 0 - H; gt[5] = 0 - H 53 gt[6] = H; gt[7] = H; gt[8] = 0 - H 54 gt[9] = 0 - H; gt[10] = H; gt[11] = 0 - H 55 gt[12] = 0 - H; gt[13] = 0 - H; gt[14] = H 56 gt[15] = H; gt[16] = 0 - H; gt[17] = H 57 gt[18] = H; gt[19] = H; gt[20] = H 58 gt[21] = 0 - H; gt[22] = H; gt[23] = H 59 gt[24] = 0; gt[25] = 0; gt[26] = 0 60 61 // 4 cameras looking at the origin 62 let ncam: i64 = 4 63 let ccx: *i64 = sys_mmap(64) as *i64 64 let ccy: *i64 = sys_mmap(64) as *i64 65 let ccz: *i64 = sys_mmap(64) as *i64 66 ccx[0] = 5000; ccy[0] = 1500; ccz[0] = 3000 67 ccx[1] = 0 - 5000; ccy[1] = 1500; ccz[1] = 3000 68 ccx[2] = 3000; ccy[2] = 1800; ccz[2] = 0 - 5000 69 ccx[3] = 0 - 3000; ccy[3] = 1200; ccz[3] = 0 - 5000 70 let f: i64 = 1200 71 let basis: *i64 = sys_mmap(64 * 4) as *i64 72 var c: i64 = 0 73 while c < ncam { 74 r3_lookat(ccx[c], ccy[c], ccz[c], 0, 0, 0, (basis as i64 + c * 72) as *i64) 75 c = c + 1 76 } 77 78 // KAT1: project gt[6] (a corner) to cam0, back-project the pixel, ray must pass through the corner 79 let uv: *i64 = sys_mmap(16) as *i64 80 let ray: *i64 = sys_mmap(32) as *i64 81 let zc0: i64 = r3_project(ccx[0], ccy[0], ccz[0], (basis as i64 + 0) as *i64, f, gt[6], gt[7], gt[8], uv) 82 r3_ray((basis as i64 + 0) as *i64, f, uv[0], uv[1], ray) 83 let rd: i64 = rg_raydist(gt[6], gt[7], gt[8], ccx[0], ccy[0], ccz[0], ray[0], ray[1], ray[2]) 84 rg_puts(" KAT1 project(corner)->pixel(" as *u8); rg_putn(uv[0]); rg_puts("," as *u8); rg_putn(uv[1]) 85 rg_puts(")->ray ; ray-to-corner dist=" as *u8); rg_putn(rd); rg_puts(" (expect ~0)\n" as *u8) 86 var t1: i64 = 0 87 if zc0 > 0 { if rd < 12 { t1 = 1 } } 88 let ig1: i64 = rg_tooth("KAT1 project<->backproject consistent (ray hits the point)" as *u8, t1, fails) 89 90 // KAT2: two rays (cam0, cam2) toward gt[6] -> midpoint recovers gt[6] 91 let uv2: *i64 = sys_mmap(16) as *i64 92 let ray2: *i64 = sys_mmap(32) as *i64 93 r3_project(ccx[2], ccy[2], ccz[2], (basis as i64 + 2 * 72) as *i64, f, gt[6], gt[7], gt[8], uv2) 94 r3_ray((basis as i64 + 2 * 72) as *i64, f, uv2[0], uv2[1], ray2) 95 let mp: *i64 = sys_mmap(32) as *i64 96 r3_midpoint(ccx[0], ccy[0], ccz[0], ray[0], ray[1], ray[2], ccx[2], ccy[2], ccz[2], ray2[0], ray2[1], ray2[2], mp) 97 let mperr: i64 = r3_dist(mp[0], mp[1], mp[2], gt[6], gt[7], gt[8]) 98 rg_puts(" KAT2 2-view midpoint = (" as *u8); rg_putn(mp[0]); rg_puts("," as *u8); rg_putn(mp[1]); rg_puts("," as *u8); rg_putn(mp[2]) 99 rg_puts(") err=" as *u8); rg_putn(mperr); rg_puts(" (corner=800,800,-800)\n" as *u8) 100 var t2: i64 = 0 101 if mperr < 20 { t2 = 1 } 102 let ig2: i64 = rg_tooth("KAT2 2-view midpoint recovers the known point" as *u8, t2, fails) 103 104 // full N-view triangulation of the whole cloud, CLEAN 105 let recon: *i64 = sys_mmap(512) as *i64 106 let rx: *i64 = sys_mmap(64) as *i64 107 let ry: *i64 = sys_mmap(64) as *i64 108 let rz: *i64 = sys_mmap(64) as *i64 109 let dxa: *i64 = sys_mmap(64) as *i64 110 let dya: *i64 = sys_mmap(64) as *i64 111 let dza: *i64 = sys_mmap(64) as *i64 112 let out3: *i64 = sys_mmap(32) as *i64 113 var noise: i64 = 0 114 var pass: i64 = 0 115 while pass < 3 { 116 // pass0 clean, pass1 +-3px noise, pass2 WRONG correspondences (neg-control) 117 var pi: i64 = 0 118 while pi < np { 119 var cc: i64 = 0 120 while cc < ncam { 121 let bb: i64 = (basis as i64 + cc * 72) 122 let uvp: *i64 = sys_mmap(16) as *i64 123 r3_project(ccx[cc], ccy[cc], ccz[cc], bb as *i64, f, gt[pi * 3], gt[pi * 3 + 1], gt[pi * 3 + 2], uvp) 124 var uu: i64 = uvp[0] 125 var vv: i64 = uvp[1] 126 if pass == 1 { 127 // deterministic +-3px pseudo-noise 128 uu = uu + ((pi * 7 + cc * 13) % 7) - 3 129 vv = vv + ((pi * 5 + cc * 11) % 7) - 3 130 } 131 if pass == 2 { 132 // NEG-CONTROL: use a DIFFERENT point's projection for some cams = wrong correspondence 133 if cc >= 2 { 134 let wrong: i64 = (pi + 3) % np 135 r3_project(ccx[cc], ccy[cc], ccz[cc], bb as *i64, f, gt[wrong * 3], gt[wrong * 3 + 1], gt[wrong * 3 + 2], uvp) 136 uu = uvp[0]; vv = uvp[1] 137 } 138 } 139 let rr: *i64 = sys_mmap(32) as *i64 140 r3_ray(bb as *i64, f, uu, vv, rr) 141 rx[cc] = ccx[cc]; ry[cc] = ccy[cc]; rz[cc] = ccz[cc] 142 dxa[cc] = rr[0]; dya[cc] = rr[1]; dza[cc] = rr[2] 143 cc = cc + 1 144 } 145 r3_triangulate_n(rx, ry, rz, dxa, dya, dza, ncam, out3) 146 recon[pi * 3] = out3[0]; recon[pi * 3 + 1] = out3[1]; recon[pi * 3 + 2] = out3[2] 147 pi = pi + 1 148 } 149 let cham: i64 = r3_chamfer(recon, np, gt, np) 150 // object diagonal ~ 2771; report permille of diagonal 151 let permille: i64 = (cham * 1000) / 2771 152 if pass == 0 { 153 rg_puts(" CLEAN Chamfer=" as *u8); rg_putn(cham); rg_puts(" units = " as *u8); rg_putn(permille); rg_puts(" permille of object diagonal\n" as *u8) 154 var t3: i64 = 0 155 if cham < 20 { t3 = 1 } 156 let ig3: i64 = rg_tooth("T3 CLEAN reconstruction near-EXACT (Chamfer < 20 units)" as *u8, t3, fails) 157 } 158 if pass == 1 { 159 rg_puts(" NOISY Chamfer=" as *u8); rg_putn(cham); rg_puts(" units = " as *u8); rg_putn(permille); rg_puts(" permille (+/-3px)\n" as *u8) 160 var t4: i64 = 0 161 if cham < 120 { t4 = 1 } // graceful: noise degrades but no blowup 162 let ig4: i64 = rg_tooth("T4 NOISY graceful (Chamfer < 120 units, no blowup)" as *u8, t4, fails) 163 noise = cham 164 } 165 if pass == 2 { 166 rg_puts(" WRONG Chamfer=" as *u8); rg_putn(cham); rg_puts(" units (neg-control: bad correspondences)\n" as *u8) 167 var t5: i64 = 0 168 if cham > 150 { t5 = 1 } // metric discriminates wrong data 169 let ig5: i64 = rg_tooth("T5 NEG-CONTROL wrong correspondences -> large Chamfer (metric works)" as *u8, t5, fails) 170 } 171 pass = pass + 1 172 } 173 174 rg_puts("\nSOTA CONTEXT: Chamfer is the DTU reconstruction metric (QGS ~0.54mm). Our CLASSICAL triangulation is\n" as *u8) 175 rg_puts("near-EXACT on clean correspondences; the 2025-26 feed-forward SOTA (DUSt3R/MASt3R/VGGT) is NEURAL =\n" as *u8) 176 rg_puts("trained-model-bound (our named gap). REMAINING sovereign rungs: feature match (nx_features) + pose est.\n" as *u8) 177 rg_puts("\nfails=" as *u8); rg_putn(fails[0]); rg_puts("\n" as *u8) 178 if fails[0] == 0 { rg_puts("verdict=GREEN -- multi-view triangulation 5/5, Chamfer-benchmarked (photo->3D classical core)\n" as *u8); return 0 } 179 rg_puts("RED\n" as *u8) 180 return 1 181}