nx_photogram_front_diag.nx source

1// nx_photogram_front_diag.nx -- measurement probe, not a gate. T5 failed with the noise control scoring 2// 18/20, which says the epipolar residual is not separating true correspondences from garbage. Before 3// touching a tolerance constant I want the actual numbers: what does e3_residual return for a TRUE 4// correspondence, for a planted outlier, and for pure noise -- both under an E built from the known pose 5// (ep_build_E) and under an E estimated from the data (e3_estimate_norm). Tuning a threshold without 6// looking at the distribution is how a gate gets faked green. 7// license_tier: ORIGINAL expect_exit: 0 8import "nx_photogram_front.nx" 9import "nx_epose.nx" 10 11const D_N: i64 = 20 12const D_OUT: i64 = 6 13 14func dw(s: *u8) -> i64 { 15 var n: i64 = 0 16 while s[n] != (0 as u8) { n = n + 1 } 17 sys_write(1, s, n) 18 return 0 19} 20 21func dn(v: i64) -> i64 { 22 let bb: *u8 = sys_mmap(32) 23 let t: *u8 = sys_mmap(32) 24 var m: i64 = v 25 var neg: i64 = 0 26 if m < 0 { neg = 1; m = 0 - m } 27 var k: i64 = 0 28 if m == 0 { t[0] = 48 as u8; k = 1 } 29 while m > 0 { t[k] = (48 + (m % 10)) as u8; m = m / 10; k = k + 1 } 30 var i: i64 = 0 31 if neg == 1 { bb[0] = 45 as u8; i = 1 } 32 var j: i64 = 0 33 while j < k { bb[i + j] = t[k - 1 - j]; j = j + 1 } 34 sys_write(1, bb, i + k) 35 return 0 36} 37 38func main() -> i64 { 39 pf_init() 40 dw("=== residual distribution probe ===\n" as *u8) 41 let Rt: *i64 = sys_mmap(9 * 8 + 64) as *i64 42 Rt[0] = 15137; Rt[1] = 0; Rt[2] = 6270 43 Rt[3] = 0; Rt[4] = 16384; Rt[5] = 0 44 Rt[6] = 0 - 6270; Rt[7] = 0; Rt[8] = 15137 45 let tt: *i64 = sys_mmap(3 * 8 + 64) as *i64 46 tt[0] = 1400 47 tt[1] = 250 48 tt[2] = 320 49 let c0: *i64 = sys_mmap(D_N * 2 * 8 + 64) as *i64 50 let c1: *i64 = sys_mmap(D_N * 2 * 8 + 64) as *i64 51 var st: i64 = 987654321 52 var n: i64 = 0 53 while n < D_N { 54 st = (st * PF_LCG_MUL + PF_LCG_ADD) & PF_LCG_MASK 55 let X: i64 = (st % 4000) - 2000 56 st = (st * PF_LCG_MUL + PF_LCG_ADD) & PF_LCG_MASK 57 let Y: i64 = (st % 4000) - 2000 58 st = (st * PF_LCG_MUL + PF_LCG_ADD) & PF_LCG_MASK 59 let Z: i64 = 4000 + (st % 4000) 60 c0[n * 2] = (X * PF_Q) / Z 61 c0[n * 2 + 1] = (Y * PF_Q) / Z 62 let p2x: i64 = (Rt[0] * X + Rt[1] * Y + Rt[2] * Z) / PF_Q + tt[0] 63 let p2y: i64 = (Rt[3] * X + Rt[4] * Y + Rt[5] * Z) / PF_Q + tt[1] 64 let p2z: i64 = (Rt[6] * X + Rt[7] * Y + Rt[8] * Z) / PF_Q + tt[2] 65 c1[n * 2] = (p2x * PF_Q) / p2z 66 c1[n * 2 + 1] = (p2y * PF_Q) / p2z 67 n = n + 1 68 } 69 // E from the KNOWN pose. NOTE the convention: for x2^T E x1 = 0 with P2 = [R | t], E = [t]_x R. 70 let Etrue: *i64 = sys_mmap(9 * 8 + 64) as *i64 71 ep_build_E(tt[0], tt[1], tt[2], Rt, Etrue) 72 dw("\n-- residual under E built from the TRUE pose (should be ~0 for all 20) --\n" as *u8) 73 var i: i64 = 0 74 while i < D_N { 75 let r: i64 = e3_residual(Etrue, c0[i * 2], c0[i * 2 + 1], c1[i * 2], c1[i * 2 + 1]) 76 dw(" pt " as *u8) 77 dn(i) 78 dw(" resid=" as *u8) 79 dn(r) 80 dw("\n" as *u8) 81 i = i + 1 82 } 83 // E estimated from all 20 clean correspondences 84 let Eest: *i64 = sys_mmap(9 * 8 + 64) as *i64 85 e3_estimate_norm(c0, c1, D_N, Eest) 86 dw("\n-- E entries: TRUE vs ESTIMATED (both max-entry normalized to Q14) --\n" as *u8) 87 var e: i64 = 0 88 while e < 9 { 89 dw(" E[" as *u8) 90 dn(e) 91 dw("] true=" as *u8) 92 dn(Etrue[e]) 93 dw(" est=" as *u8) 94 dn(Eest[e]) 95 dw("\n" as *u8) 96 e = e + 1 97 } 98 dw("\n-- residual under the ESTIMATED E (clean data) --\n" as *u8) 99 var mx: i64 = 0 100 i = 0 101 while i < D_N { 102 let r2: i64 = e3_residual(Eest, c0[i * 2], c0[i * 2 + 1], c1[i * 2], c1[i * 2 + 1]) 103 if r2 > mx { mx = r2 } 104 dw(" pt " as *u8) 105 dn(i) 106 dw(" resid=" as *u8) 107 dn(r2) 108 dw("\n" as *u8) 109 i = i + 1 110 } 111 dw(" max clean residual = " as *u8) 112 dn(mx) 113 dw("\n" as *u8) 114 // now: what does a GARBAGE correspondence score under the same estimated E? 115 dw("\n-- residual for RANDOM (garbage) pairs under the estimated E --\n" as *u8) 116 var g: i64 = 0 117 var gmin: i64 = 0 - 1 118 while g < 12 { 119 st = (st * PF_LCG_MUL + PF_LCG_ADD) & PF_LCG_MASK 120 let ax: i64 = (st % 8000) - 4000 121 st = (st * PF_LCG_MUL + PF_LCG_ADD) & PF_LCG_MASK 122 let ay: i64 = (st % 8000) - 4000 123 st = (st * PF_LCG_MUL + PF_LCG_ADD) & PF_LCG_MASK 124 let bx: i64 = (st % 8000) - 4000 125 st = (st * PF_LCG_MUL + PF_LCG_ADD) & PF_LCG_MASK 126 let by: i64 = (st % 8000) - 4000 127 let r3: i64 = e3_residual(Eest, ax, ay, bx, by) 128 if gmin < 0 { gmin = r3 } else { if r3 < gmin { gmin = r3 } } 129 dw(" rnd " as *u8) 130 dn(g) 131 dw(" resid=" as *u8) 132 dn(r3) 133 dw("\n" as *u8) 134 g = g + 1 135 } 136 dw(" min garbage residual = " as *u8) 137 dn(gmin) 138 dw("\n\nSEPARATION (e3_estimate_norm): clean max=" as *u8) 139 dn(mx) 140 dw(" garbage min=" as *u8) 141 dn(gmin) 142 dw("\n" as *u8) 143 144 // ---- A/B: the same data through pf_estimate_e (repeated squaring) ---- 145 dw("\n=== A/B: pf_estimate_e (repeated squaring) on the SAME data ===\n" as *u8) 146 let Enew: *i64 = sys_mmap(9 * 8 + 64) as *i64 147 pf_estimate_e(c0, c1, D_N, Enew) 148 var e2: i64 = 0 149 while e2 < 9 { 150 dw(" E[" as *u8) 151 dn(e2) 152 dw("] true=" as *u8) 153 dn(Etrue[e2]) 154 dw(" pf_est=" as *u8) 155 dn(Enew[e2]) 156 dw("\n" as *u8) 157 e2 = e2 + 1 158 } 159 var mx2: i64 = 0 160 i = 0 161 while i < D_N { 162 let r4: i64 = e3_residual(Enew, c0[i * 2], c0[i * 2 + 1], c1[i * 2], c1[i * 2 + 1]) 163 if r4 > mx2 { mx2 = r4 } 164 i = i + 1 165 } 166 var gmin2: i64 = 0 - 1 167 var g2: i64 = 0 168 var st2: i64 = 55555555 169 while g2 < 200 { 170 st2 = (st2 * PF_LCG_MUL + PF_LCG_ADD) & PF_LCG_MASK 171 let ax2: i64 = (st2 % 8000) - 4000 172 st2 = (st2 * PF_LCG_MUL + PF_LCG_ADD) & PF_LCG_MASK 173 let ay2: i64 = (st2 % 8000) - 4000 174 st2 = (st2 * PF_LCG_MUL + PF_LCG_ADD) & PF_LCG_MASK 175 let bx2: i64 = (st2 % 8000) - 4000 176 st2 = (st2 * PF_LCG_MUL + PF_LCG_ADD) & PF_LCG_MASK 177 let by2: i64 = (st2 % 8000) - 4000 178 let r5: i64 = e3_residual(Enew, ax2, ay2, bx2, by2) 179 if gmin2 < 0 { gmin2 = r5 } else { if r5 < gmin2 { gmin2 = r5 } } 180 g2 = g2 + 1 181 } 182 dw("\nSEPARATION (pf_estimate_e): clean max=" as *u8) 183 dn(mx2) 184 dw(" garbage min over 200 draws=" as *u8) 185 dn(gmin2) 186 dw("\n" as *u8) 187 dw("VERDICT: " as *u8) 188 if mx2 < mx { dw("pf_estimate_e is MORE accurate on clean data" as *u8) } else { dw("NO improvement on clean data" as *u8) } 189 dw("\n" as *u8) 190 return 0 191}