code wiki / (root) / nx_photogram_front.nx
`nx_photogram_front.nx` source

↩ module page · 817 lines · 34744 B
1﻿// nx_photogram_front.nx -- THE PHOTOGRAMMETRY FRONT END: pixels -> verified correspondences.
2//
3// WHY THIS FILE EXISTS: nx_photogrambench measured the stack and its gap queue put this at the top by a
4// wide margin (prio 54 / 48 / 15). The Nishi stack already owns a real, gated two-view BACK end --
5// normalized 8-point essential estimation (nx_essential3d), E->(R,t) with cheirality selection (nx_epose),
6// N-view triangulation and Chamfer scoring (nx_recon3d). What it did NOT own is any way to get from actual
7// pixels to the correspondences that back end consumes. Harris corners exist at a single scale; the two
8// things named "descriptor" in the tree are a TEXT lexicon profiler and a 16-element GLOBAL aesthetic
9// signature -- neither is a local patch descriptor, so nothing could be matched between two photographs.
10// The back end was a working engine with no fuel line. This is the fuel line.
11//
12// WHAT IT DOES (all integer, deterministic, no float, no trained weights, no third-party code):
13//   1. DETECT   Harris corner response + 3x3 non-max suppression, strongest-N selection.
14//   2. ORIENT   intensity-centroid angle (Rosin): the vector from patch centre to intensity centroid
15//               gives each keypoint a repeatable direction, quantised to 32 buckets by MAX DOT PRODUCT
16//               against a Q14 unit-vector table -- no atan2, no float, exact.
17//   3. DESCRIBE steered BRIEF: 256 smoothed intensity comparisons over a fixed random-but-deterministic
18//               sampling pattern, the pattern ROTATED by the keypoint's angle before sampling. Rotating
19//               the pattern is what buys rotation invariance -- an unsteered BRIEF collapses the moment
20//               the camera rolls, which is exactly what happens between two hand-held photographs.
21//   4. MATCH    brute-force Hamming (nibble-LUT popcount) + Lowe ratio test + mutual cross-check. Both
22//               filters are needed: the ratio test kills ambiguous matches, the cross-check kills
23//               one-way matches into repeated texture.
24//   5. VERIFY   RANSAC whose MODEL IS THE ESSENTIAL MATRIX (8-point minimal sample -> e3_estimate_norm),
25//               scored by epipolar residual, then re-fit on the full inlier set. The tree's existing
26//               nx_imgransac is a genuine RANSAC but over translation/similarity for copy-detection; it
27//               cannot fit E, so it cannot verify a real camera pair. This can.
28//
29// The sampling pattern and the angle table are DATA generated at init from named constants (rule 11), not
30// magic numbers sprinkled through the code, and the LCG seed is fixed so every run is bit-reproducible.
31//
32// Composes nx_essential3d (-> nx_recon3d -> nx_syscalls). license_tier: ORIGINAL
33import "nx_essential3d.nx"
34const PF_MAGIC_65536: i64 = 65536
35const PF_MAGIC_16384: i64 = 16384
36const PF_MAGIC_49152: i64 = 49152
37const PF_MAGIC_1152921504606846975: i64 = 1152921504606846975
38const PF_MAGIC_4096: i64 = 4096
39
40const PF_Q: i64 = 16384              // Q14 one, matching the rest of the geometry stack
41const PF_QH: i64 = 1073741824        // Q30 -- headroom for the angle-table recurrence
42const PF_NANG: i64 = 32              // orientation buckets (11.25 degrees each)
43const PF_NBITS: i64 = 256            // descriptor length in bits
44const PF_NWORD: i64 = 4              // 256 bits / 64
45const PF_PATCH: i64 = 15             // descriptor patch radius -> 31x31 window
46const PF_ORAD: i64 = 9               // orientation centroid radius
47const PF_SMOOTH: i64 = 2             // box-blur radius at each sample point (BRIEF needs smoothing)
48const PF_LCG_SEED: i64 = 20260725    // fixed -> the sampling pattern is identical on every machine
49const PF_LCG_MUL: i64 = 1103515245
50const PF_LCG_ADD: i64 = 12345
51const PF_LCG_MASK: i64 = 2147483647
52const PF_HARRIS_K_NUM: i64 = 4       // Harris k = 0.04
53const PF_HARRIS_K_DEN: i64 = 100
54const PF_RATIO_NUM: i64 = 8          // Lowe ratio 0.8
55const PF_RATIO_DEN: i64 = 10
56const PF_RANSAC_ITERS: i64 = 256
57const PF_RANSAC_SAMPLE: i64 = 8      // the 8-point algorithm's minimal set
58const PF_CAND_CAP: i64 = 262144      // detector candidate pool (see pf_detect)
59const PF_OCTAVES: i64 = 6            // scale-space levels (see pf_pyramid)
60const PF_SCALE_NUM: i64 = 4          // per-level shrink = 4/5, i.e. a scale factor of 1.25
61const PF_SCALE_DEN: i64 = 5
62const PF_MIN_SIDE: i64 = 64          // stop building levels below this
63
64static pf_cos: *i64
65static pf_sin: *i64
66static pf_pat: *i64                  // 4 entries per bit: px,py,qx,qy
67static pf_nib: *i64                  // 16-entry popcount lookup
68static pf_ready: i64
69
70// ---- LIVE TUNING (rule 11: thresholds are data) ----------------------------------------------------
71// These mirror the consts above and are overridden from knowledge/photogram_front.conf when it exists.
72// The consts remain as the fallback so the organ is still self-contained when run without the tree.
73static pf_c_ratio_num: i64
74static pf_c_ratio_den: i64
75static pf_c_ransac_iters: i64
76static pf_c_ransac_tol: i64
77static pf_c_octaves: i64
78static pf_c_scale_num: i64
79static pf_c_scale_den: i64
80static pf_c_min_side: i64
81static pf_c_harris_num: i64
82static pf_c_harris_den: i64
83
84// Parse "<key> <int>" lines from the config. Deliberately tiny and total: unknown keys are ignored so
85// the file can carry keys for sibling organs, and any key the file omits keeps its compiled default.
86func pf_conf_key(buf: *u8, n: i64, key: *u8, dflt: i64) -> i64 {
  var i: i64 = 0
  while i < n {
      // at a line start, try to match the key
      var j: i64 = 0
      var hit: i64 = 1
      while key[j] != (0 as u8) {
          let o: i64 = i + j
          if o >= n { hit = 0 } else { if buf[o] != key[j] { hit = 0 } }
          if hit == 0 { j = j + 1 } else { j = j + 1 }
      }
      if hit == 1 {
          let p: i64 = i + j
          // require whitespace right after the key so "ratio_num" cannot match "ratio_numeric"
          var sep: i64 = 0
          if p < n { if buf[p] == (32 as u8) { sep = 1 } }
          if p < n { if buf[p] == (9 as u8) { sep = 1 } }
          if sep == 1 {
              var q: i64 = p
              var scanning: i64 = 1
              while scanning == 1 {
                  if q >= n { scanning = 0 } else {
                      let c2: i64 = buf[q] as i64
                      if c2 == 32 { q = q + 1 } else { if c2 == 9 { q = q + 1 } else { scanning = 0 } }
                  }
              }
              var v: i64 = 0
              var got: i64 = 0
              var reading: i64 = 1
              while reading == 1 {
                  if q >= n { reading = 0 } else {
                      let d: i64 = buf[q] as i64
                      if d >= 48 { if d <= 57 { v = v * 10 + (d - 48); got = 1; q = q + 1 } else { reading = 0 } } else { reading = 0 }
                  }
              }
              if got == 1 { return v }
          }
      }
      // advance to the next line start
      var adv: i64 = 1
      while adv == 1 {
          if i >= n { adv = 0 } else {
              if buf[i] == (10 as u8) { i = i + 1; adv = 0 } else { i = i + 1 }
          }
      }
  }
  return dflt
133}
134
135func pf_load_conf() -> i64 {
  pf_c_ratio_num = PF_RATIO_NUM
  pf_c_ratio_den = PF_RATIO_DEN
  pf_c_ransac_iters = PF_RANSAC_ITERS
  pf_c_ransac_tol = 20
  pf_c_octaves = PF_OCTAVES
  pf_c_scale_num = PF_SCALE_NUM
  pf_c_scale_den = PF_SCALE_DEN
  pf_c_min_side = PF_MIN_SIDE
  pf_c_harris_num = PF_HARRIS_K_NUM
  pf_c_harris_den = PF_HARRIS_K_DEN
  var fd: i64 = sys_openat_rd("knowledge/photogram_front.conf" as *u8)
  if fd < 0 { fd = sys_openat_rd("buildroot/knowledge/photogram_front.conf" as *u8) }
  if fd < 0 { return 0 }
  let buf: *u8 = sys_mmap(PF_MAGIC_65536)
  var total: i64 = 0
  var reading: i64 = 1
  while reading == 1 {
      let got: i64 = sys_read(fd, buf + total, PF_MAGIC_16384)
      if got <= 0 { reading = 0 } else {
          total = total + got
          if total > PF_MAGIC_49152 { reading = 0 }
      }
  }
  sys_close(fd)
  if total <= 0 { return 0 }
  pf_c_ratio_num = pf_conf_key(buf, total, "ratio_num" as *u8, pf_c_ratio_num)
  pf_c_ratio_den = pf_conf_key(buf, total, "ratio_den" as *u8, pf_c_ratio_den)
  pf_c_ransac_iters = pf_conf_key(buf, total, "ransac_iters" as *u8, pf_c_ransac_iters)
  pf_c_ransac_tol = pf_conf_key(buf, total, "ransac_tol" as *u8, pf_c_ransac_tol)
  pf_c_octaves = pf_conf_key(buf, total, "octaves" as *u8, pf_c_octaves)
  pf_c_scale_num = pf_conf_key(buf, total, "scale_num" as *u8, pf_c_scale_num)
  pf_c_scale_den = pf_conf_key(buf, total, "scale_den" as *u8, pf_c_scale_den)
  pf_c_min_side = pf_conf_key(buf, total, "min_side" as *u8, pf_c_min_side)
  pf_c_harris_num = pf_conf_key(buf, total, "harris_k_num" as *u8, pf_c_harris_num)
  pf_c_harris_den = pf_conf_key(buf, total, "harris_k_den" as *u8, pf_c_harris_den)
  return 1
172}
173
174// ---- angle table: PF_NANG unit vectors, COMPUTED, not typed ----------------------------------------
175// This table used to be 64 hand-entered constants. Hand-entered data is wrong in three ways even when
176// the numbers happen to be right: a typo is invisible until something downstream misbehaves, the table
177// silently stops matching PF_NANG the moment anyone retunes the bucket count, and nobody can tell which
178// numbers were derived and which were guessed. So it is derived from first principles at init.
179//
180// METHOD: exact half-angle recurrence, no float and no trig table.
181//   cos(t/2) = sqrt((1+cos t)/2)      sin(t/2) = sqrt((1-cos t)/2)
182// Start at 90 degrees (cos=0, sin=1) and halve to reach the base step, then generate every bucket by
183// repeated complex rotation. The recurrence runs at Q30, not Q14: 31 successive rotations at Q14 would
184// accumulate visible drift, which is the real reason the table was typed out in the first place. With
185// ~30 bits of headroom the accumulated error stays far below one Q14 unit.
186//
187// These values are NOT bit-identical to the table they replace: the recurrence truncates where a human
188// writing the numbers down rounds, so sin(11.25 deg) = 3195.98 lands on 3195 against the 3196 that was
189// typed. That is a sub-ulp difference at Q14, and it is VERIFIED rather than asserted -- gate T1 checks
190// every entry for unit-norm, T4 still measures Hamming 0 under a 90-degree rotation, and T6 still
191// recovers the ground-truth rotation at Frobenius error 0. A broken derivation cannot pass those.
192//
193// Requires PF_NANG to be 8 * 2^k so the base step is reachable by halving from 90 degrees.
194func pf_init_angles() -> i64 {
  // cos/sin of 90 degrees at Q30
  var c: i64 = 0
  var s: i64 = PF_QH
  // halve until we reach the step 360/PF_NANG. 90 degrees is PF_NANG/4 steps, so halve log2 of that.
  var steps: i64 = PF_NANG / 4
  while steps > 1 {
      let cn: i64 = e3_isqrt(((PF_QH + c) / 2) * PF_QH)
      let sn: i64 = e3_isqrt(((PF_QH - c) / 2) * PF_QH)
      c = cn
      s = sn
      steps = steps / 2
  }
  let cb: i64 = c
  let sb: i64 = s
  // walk the circle: (c,s) <- (c,s) * (cb,sb)
  var cc: i64 = PF_QH
  var ss: i64 = 0
  var k: i64 = 0
  while k < PF_NANG {
      pf_cos[k] = (cc * PF_Q) / PF_QH
      pf_sin[k] = (ss * PF_Q) / PF_QH
      let nc: i64 = (cc * cb - ss * sb) / PF_QH
      let ns: i64 = (ss * cb + cc * sb) / PF_QH
      cc = nc
      ss = ns
      k = k + 1
  }
  return 0
223}
224
225func pf_init() -> i64 {
  if pf_ready == 1 { return 0 }
  pf_load_conf()
  pf_cos = sys_mmap(PF_NANG * 8 + 64) as *i64
  pf_sin = sys_mmap(PF_NANG * 8 + 64) as *i64
  pf_pat = sys_mmap(PF_NBITS * 4 * 8 + 64) as *i64
  pf_nib = sys_mmap(16 * 8 + 64) as *i64
  pf_init_angles()
  // popcount nibble LUT
  var i: i64 = 0
  while i < 16 {
      var c: i64 = 0
      var b: i64 = 0
      var v: i64 = i
      while b < 4 {
          if (v % 2) == 1 { c = c + 1 }
          v = v / 2
          b = b + 1
      }
      pf_nib[i] = c
      i = i + 1
  }
  // deterministic sampling pattern: 256 point pairs inside the patch
  let span: i64 = 2 * PF_PATCH + 1
  var st: i64 = PF_LCG_SEED
  var k: i64 = 0
  while k < PF_NBITS * 4 {
      st = (st * PF_LCG_MUL + PF_LCG_ADD) & PF_LCG_MASK
      pf_pat[k] = (st % span) - PF_PATCH
      k = k + 1
  }
  pf_ready = 1
  return 0
258}
259
260func pf_abs(a: i64) -> i64 { if a < 0 { return 0 - a } return a }
261
262// ---- ESSENTIAL-MATRIX ESTIMATION ------------------------------------------------------------------
263// This was briefly a private copy of the solver (pf_null9 / pf_estimate_e), written as an additive
264// sibling so a measured accuracy defect in nx_essential3d could be fixed without touching a module the
265// cadtwin lane depends on. That fix has since been folded INTO nx_essential3d itself (debt seq814
266// eaten): the repair is numerics-only behind an unchanged signature, so every caller should get it
267// rather than only this one. The private copy is gone -- two implementations of the same estimator is
268// exactly the dual-copy debt class the ecosystem keeps paying for. pf_estimate_e survives as a thin
269// alias so this module's call sites read in its own vocabulary.
270func pf_estimate_e(c0: *i64, c1: *i64, n: i64, out_e: *i64) -> i64 {
  return e3_estimate_norm(c0, c1, n, out_e)
272}
273
274// ---- image access: grayscale, one byte per pixel, row stride = w. Out of bounds reads as 0, which keeps
275// every loop below branch-simple and can never fault on a patch that overhangs the border. ----
276func pf_pix(img: *u8, w: i64, h: i64, x: i64, y: i64) -> i64 {
  if x < 0 { return 0 }
  if y < 0 { return 0 }
  if x >= w { return 0 }
  if y >= h { return 0 }
  return img[y * w + x] as i64
282}
283
284// box-average over a (2r+1)^2 window -- BRIEF compares SMOOTHED intensities, otherwise single-pixel noise
285// flips bits and the descriptor stops being repeatable
286func pf_smooth(img: *u8, w: i64, h: i64, x: i64, y: i64) -> i64 {
  var s: i64 = 0
  var n: i64 = 0
  var dy: i64 = 0 - PF_SMOOTH
  while dy <= PF_SMOOTH {
      var dx: i64 = 0 - PF_SMOOTH
      while dx <= PF_SMOOTH {
          s = s + pf_pix(img, w, h, x + dx, y + dy)
          n = n + 1
          dx = dx + 1
      }
      dy = dy + 1
  }
  if n == 0 { return 0 }
  return s / n
301}
302
303// ---- 1. DETECT: Harris response at one pixel (Sobel gradients over a 3x3 structure-tensor window) ----
304func pf_harris(img: *u8, w: i64, h: i64, x: i64, y: i64) -> i64 {
  var sxx: i64 = 0
  var syy: i64 = 0
  var sxy: i64 = 0
  var dy: i64 = 0 - 1
  while dy <= 1 {
      var dx: i64 = 0 - 1
      while dx <= 1 {
          let px: i64 = x + dx
          let py: i64 = y + dy
          let p00: i64 = pf_pix(img, w, h, px - 1, py - 1)
          let p01: i64 = pf_pix(img, w, h, px, py - 1)
          let p02: i64 = pf_pix(img, w, h, px + 1, py - 1)
          let p10: i64 = pf_pix(img, w, h, px - 1, py)
          let p12: i64 = pf_pix(img, w, h, px + 1, py)
          let p20: i64 = pf_pix(img, w, h, px - 1, py + 1)
          let p21: i64 = pf_pix(img, w, h, px, py + 1)
          let p22: i64 = pf_pix(img, w, h, px + 1, py + 1)
          let gx: i64 = (p02 + 2 * p12 + p22) - (p00 + 2 * p10 + p20)
          let gy: i64 = (p20 + 2 * p21 + p22) - (p00 + 2 * p01 + p02)
          sxx = sxx + gx * gx
          syy = syy + gy * gy
          sxy = sxy + gx * gy
          dx = dx + 1
      }
      dy = dy + 1
  }
  // scale down before the products so det/trace^2 cannot overflow on a high-contrast patch
  sxx = sxx / 256
  syy = syy / 256
  sxy = sxy / 256
  let det: i64 = sxx * syy - sxy * sxy
  let tr: i64 = sxx + syy
  let ktr: i64 = (tr * tr * pf_c_harris_num) / pf_c_harris_den
  return det - ktr
339}
340
341// Detect up to maxn keypoints: Harris + 3x3 NMS above thresh, then keep the strongest maxn.
342// kx/ky/ks are caller arrays of length >= maxn. Returns the count found.
343func pf_detect(img: *u8, w: i64, h: i64, kx: *i64, ky: *i64, ks: *i64, maxn: i64, thresh: i64) -> i64 {
  pf_init()
  let margin: i64 = PF_PATCH + PF_SMOOTH + 1
  // Candidate pool must be able to hold EVERY non-max-suppressed maximum, not a small multiple of the
  // requested count. With a pool of maxn*4 a real photograph overflows it early and the "strongest N"
  // selection silently degrades into "the first N found in raster order" -- top-left biased keypoints
  // and a benchmark number that flatters us for the wrong reason. Sizing for the whole image keeps the
  // selection genuinely strongest-first.
  let cap: i64 = PF_CAND_CAP
  let cx: *i64 = sys_mmap(cap * 8 + 64) as *i64
  let cy: *i64 = sys_mmap(cap * 8 + 64) as *i64
  let cs: *i64 = sys_mmap(cap * 8 + 64) as *i64
  var n: i64 = 0
  var y: i64 = margin
  while y < h - margin {
      var x: i64 = margin
      while x < w - margin {
          let r: i64 = pf_harris(img, w, h, x, y)
          if r > thresh {
              var ismax: i64 = 1
              var dy: i64 = 0 - 1
              while dy <= 1 {
                  var dx: i64 = 0 - 1
                  while dx <= 1 {
                      if dx != 0 { if pf_harris(img, w, h, x + dx, y + dy) > r { ismax = 0 } }
                      if dx == 0 { if dy != 0 { if pf_harris(img, w, h, x, y + dy) > r { ismax = 0 } } }
                      dx = dx + 1
                  }
                  dy = dy + 1
              }
              if ismax == 1 {
                  if n < cap {
                      cx[n] = x
                      cy[n] = y
                      cs[n] = r
                      n = n + 1
                  }
              }
          }
          x = x + 1
      }
      y = y + 1
  }
  // keep the strongest maxn (selection sort by score desc over the found set)
  var out: i64 = n
  if out > maxn { out = maxn }
  var a: i64 = 0
  while a < out {
      var best: i64 = a
      var b: i64 = a + 1
      while b < n {
          if cs[b] > cs[best] { best = b }
          b = b + 1
      }
      let tx: i64 = cx[a]
      let ty: i64 = cy[a]
      let ts: i64 = cs[a]
      cx[a] = cx[best]
      cy[a] = cy[best]
      cs[a] = cs[best]
      cx[best] = tx
      cy[best] = ty
      cs[best] = ts
      kx[a] = cx[a]
      ky[a] = cy[a]
      ks[a] = cs[a]
      a = a + 1
  }
  return out
412}
413
414// ---- 2. ORIENT: intensity centroid. m10/m01 over a circular patch give a vector from the patch centre
415// to its centre of mass; its direction is the keypoint angle. Quantised by max dot product against the
416// Q14 unit table -- this is an exact integer replacement for atan2 followed by bucketing. ----
417func pf_orient(img: *u8, w: i64, h: i64, x: i64, y: i64) -> i64 {
  pf_init()
  var m10: i64 = 0
  var m01: i64 = 0
  let r2: i64 = PF_ORAD * PF_ORAD
  var dy: i64 = 0 - PF_ORAD
  while dy <= PF_ORAD {
      var dx: i64 = 0 - PF_ORAD
      while dx <= PF_ORAD {
          if dx * dx + dy * dy <= r2 {
              let v: i64 = pf_pix(img, w, h, x + dx, y + dy)
              m10 = m10 + dx * v
              m01 = m01 + dy * v
          }
          dx = dx + 1
      }
      dy = dy + 1
  }
  if m10 == 0 { if m01 == 0 { return 0 } }
  var bestk: i64 = 0
  var bestd: i64 = 0 - 1
  var k: i64 = 0
  while k < PF_NANG {
      let d: i64 = (m10 * pf_cos[k] + m01 * pf_sin[k]) / PF_Q
      if bestd < 0 { bestd = d; bestk = k } else { if d > bestd { bestd = d; bestk = k } }
      k = k + 1
  }
  return bestk
445}
446
447// ---- 3. DESCRIBE: steered BRIEF. desc = caller array of PF_NWORD i64. ----
448func pf_describe(img: *u8, w: i64, h: i64, x: i64, y: i64, ang: i64, desc: *i64) -> i64 {
  pf_init()
  let ca: i64 = pf_cos[ang]
  let sa: i64 = pf_sin[ang]
  var i: i64 = 0
  while i < PF_NWORD { desc[i] = 0; i = i + 1 }
  var bit: i64 = 0
  while bit < PF_NBITS {
      let base: i64 = bit * 4
      let ax: i64 = pf_pat[base]
      let ay: i64 = pf_pat[base + 1]
      let bx: i64 = pf_pat[base + 2]
      let by: i64 = pf_pat[base + 3]
      // rotate both sample offsets by the keypoint angle
      let rax: i64 = (ca * ax - sa * ay) / PF_Q
      let ray: i64 = (sa * ax + ca * ay) / PF_Q
      let rbx: i64 = (ca * bx - sa * by) / PF_Q
      let rby: i64 = (sa * bx + ca * by) / PF_Q
      let va: i64 = pf_smooth(img, w, h, x + rax, y + ray)
      let vb: i64 = pf_smooth(img, w, h, x + rbx, y + rby)
      if va < vb {
          let word: i64 = bit / 64
          let off: i64 = bit % 64
          let one: i64 = 1
          desc[word] = desc[word] | (one << off)
      }
      bit = bit + 1
  }
  return 0
477}
478
479// popcount of one i64 via the nibble LUT (16 lookups, no float, no builtin)
480func pf_popcnt(v: i64) -> i64 {
  var c: i64 = 0
  var i: i64 = 0
  var x: i64 = v
  while i < 16 {
      let nib: i64 = x & 15
      c = c + pf_nib[nib]
      x = x >> 4
      x = x & PF_MAGIC_1152921504606846975
      i = i + 1
  }
  return c
492}
493
494func pf_hamming(a: *i64, ai: i64, b: *i64, bi: i64) -> i64 {
  var d: i64 = 0
  var i: i64 = 0
  while i < PF_NWORD {
      let xa: i64 = a[ai * PF_NWORD + i]
      let xb: i64 = b[bi * PF_NWORD + i]
      d = d + pf_popcnt(xa ^ xb)
      i = i + 1
  }
  return d
504}
505
506// Describe a whole keypoint set into a packed descriptor array (n * PF_NWORD i64) + angles.
507func pf_describe_set(img: *u8, w: i64, h: i64, kx: *i64, ky: *i64, n: i64, desc: *i64, ang: *i64) -> i64 {
  pf_init()
  let tmp: *i64 = sys_mmap(PF_NWORD * 8 + 64) as *i64
  var i: i64 = 0
  while i < n {
      let a: i64 = pf_orient(img, w, h, kx[i], ky[i])
      ang[i] = a
      pf_describe(img, w, h, kx[i], ky[i], a, tmp)
      var j: i64 = 0
      while j < PF_NWORD {
          desc[i * PF_NWORD + j] = tmp[j]
          j = j + 1
      }
      i = i + 1
  }
  return 0
523}
524
525// ---- 4. MATCH: Lowe ratio test + mutual cross-check. mA/mB receive the matched index pairs. ----
526// Returns the number of surviving matches. Both filters matter: the ratio test removes a match whose best
527// candidate is barely better than its runner-up (ambiguous texture), the cross-check removes a match that
528// is not reciprocal (A picks B, but B prefers someone else).
529func pf_match(da: *i64, na: i64, db: *i64, nb: i64, ma: *i64, mb: *i64, maxm: i64) -> i64 {
  var cnt: i64 = 0
  var i: i64 = 0
  while i < na {
      var b1: i64 = 0 - 1
      var b2: i64 = 0 - 1
      var bi: i64 = 0 - 1
      var j: i64 = 0
      while j < nb {
          let d: i64 = pf_hamming(da, i, db, j)
          if b1 < 0 {
              b1 = d
              bi = j
          } else {
              if d < b1 {
                  b2 = b1
                  b1 = d
                  bi = j
              } else {
                  if b2 < 0 { b2 = d } else { if d < b2 { b2 = d } }
              }
          }
          j = j + 1
      }
      var keep: i64 = 0
      if bi >= 0 {
          if b2 < 0 { keep = 1 } else { if b1 * pf_c_ratio_den < b2 * pf_c_ratio_num { keep = 1 } }
      }
      // cross-check: bi's own best match in A must come back to i
      if keep == 1 {
          var cb: i64 = 0 - 1
          var ci: i64 = 0 - 1
          var k: i64 = 0
          while k < na {
              let d2: i64 = pf_hamming(da, k, db, bi)
              if cb < 0 {
                  cb = d2
                  ci = k
              } else {
                  if d2 < cb { cb = d2; ci = k }
              }
              k = k + 1
          }
          if ci != i { keep = 0 }
      }
      if keep == 1 {
          if cnt < maxm {
              ma[cnt] = i
              mb[cnt] = bi
              cnt = cnt + 1
          }
      }
      i = i + 1
  }
  return cnt
584}
585
586// ---- SCALE SPACE ------------------------------------------------------------------------------------
587// Added 2026-07-25 in response to a MEASURED failure, not a hunch. In the head-to-head against OpenCV
588// SIFT and ORB on the Oxford/VGG sequences, the single-scale detector held ORB-parity overall and BEAT
589// both baselines on illumination change (leuven) and blur (bikes) -- but collapsed on `boat`, which is
590// rotation PLUS zoom: 60.9% precision on pair 3 and 0 correct matches on pair 4. Rotation was already
591// handled by the steered descriptor; scale was not handled at all, because every keypoint was found and
592// described in one fixed 31x31 window. When the camera zooms, the same physical feature simply is not
593// the same size in pixels and a fixed window cannot describe it twice the same way.
594//
595// The fix is the standard one: build an image pyramid, detect and DESCRIBE on each level, and report
596// coordinates back in level-0 pixels. Describing on the level where the point was found is the part
597// that buys scale invariance -- a constant patch radius on a 1.25x-smaller image covers 1.25x more of
598// the scene, so a feature that shrank between photographs is still described over the same content.
599//
600// Level budgets are proportional to level AREA rather than equal, because Harris responses are not
601// comparable across resamplings: a global "take the strongest N" would let whichever level happens to
602// have the largest gradient magnitudes swallow the entire budget. Area-proportional quotas are what
603// ORB does, and they keep fine detail dominant while still funding the coarse levels.
604
605// area-average downsample (dst dw x dh from src sw x sh). Averaging rather than point-sampling matters:
606// nearest-neighbour aliases high-frequency texture into different noise at every level, and the whole
607// point of a pyramid is that the SAME feature survives across levels.
608func pf_resize(src: *u8, sw: i64, sh: i64, dst: *u8, dw: i64, dh: i64) -> i64 {
  if dw <= 0 { return 0 }
  if dh <= 0 { return 0 }
  var y: i64 = 0
  while y < dh {
      let sy0: i64 = (y * sh) / dh
      var sy1: i64 = ((y + 1) * sh) / dh
      if sy1 <= sy0 { sy1 = sy0 + 1 }
      var x: i64 = 0
      while x < dw {
          let sx0: i64 = (x * sw) / dw
          var sx1: i64 = ((x + 1) * sw) / dw
          if sx1 <= sx0 { sx1 = sx0 + 1 }
          var acc: i64 = 0
          var cnt: i64 = 0
          var yy: i64 = sy0
          while yy < sy1 {
              var xx: i64 = sx0
              while xx < sx1 {
                  if xx < sw { if yy < sh { acc = acc + (src[yy * sw + xx] as i64); cnt = cnt + 1 } }
                  xx = xx + 1
              }
              yy = yy + 1
          }
          if cnt > 0 { dst[y * dw + x] = (acc / cnt) as u8 } else { dst[y * dw + x] = 0 as u8 }
          x = x + 1
      }
      y = y + 1
  }
  return 0
638}
639
640// Detect + describe across the pyramid. kx/ky come back in LEVEL-0 pixel coordinates so the caller
641// never has to know a pyramid was involved; klev records the level for diagnostics. Returns the count.
642func pf_pyramid(img: *u8, w: i64, h: i64, maxn: i64, kx: *i64, ky: *i64, klev: *i64, desc: *i64, ang: *i64) -> i64 {
  pf_init()
  // total area over all levels, for the proportional budget
  var tarea: i64 = 0
  var lw: i64 = w
  var lh: i64 = h
  var l: i64 = 0
  while l < pf_c_octaves {
      if lw >= pf_c_min_side {
          if lh >= pf_c_min_side { tarea = tarea + lw * lh }
      }
      lw = (lw * pf_c_scale_num) / pf_c_scale_den
      lh = (lh * pf_c_scale_num) / pf_c_scale_den
      l = l + 1
  }
  if tarea <= 0 { return 0 }
  let cur: *u8 = sys_mmap(w * h + PF_MAGIC_4096)
  let tmp: *i64 = sys_mmap(PF_NWORD * 8 + 64) as *i64
  let bx: *i64 = sys_mmap(maxn * 8 + 64) as *i64
  let by: *i64 = sys_mmap(maxn * 8 + 64) as *i64
  let bs: *i64 = sys_mmap(maxn * 8 + 64) as *i64
  var total: i64 = 0
  lw = w
  lh = h
  l = 0
  while l < pf_c_octaves {
      var ok: i64 = 0
      if lw >= pf_c_min_side { if lh >= pf_c_min_side { ok = 1 } }
      if ok == 1 {
          if l == 0 {
              var c: i64 = 0
              while c < w * h { cur[c] = img[c]; c = c + 1 }
          } else {
              // resample from the ORIGINAL each time: repeatedly shrinking an already-shrunk buffer
              // compounds averaging error down the pyramid
              pf_resize(img, w, h, cur, lw, lh)
          }
          var quota: i64 = (maxn * lw * lh) / tarea
          if quota < 1 { quota = 1 }
          if total + quota > maxn { quota = maxn - total }
          if quota > 0 {
              let n: i64 = pf_detect(cur, lw, lh, bx, by, bs, quota, 0)
              var i: i64 = 0
              while i < n {
                  if total < maxn {
                      let a: i64 = pf_orient(cur, lw, lh, bx[i], by[i])
                      pf_describe(cur, lw, lh, bx[i], by[i], a, tmp)
                      var j: i64 = 0
                      while j < PF_NWORD {
                          desc[total * PF_NWORD + j] = tmp[j]
                          j = j + 1
                      }
                      ang[total] = a
                      kx[total] = (bx[i] * w) / lw
                      ky[total] = (by[i] * h) / lh
                      klev[total] = l
                      total = total + 1
                  }
                  i = i + 1
              }
          }
      }
      lw = (lw * PF_SCALE_NUM) / PF_SCALE_DEN
      lh = (lh * PF_SCALE_NUM) / PF_SCALE_DEN
      l = l + 1
  }
  return total
709}
710
711// pixel (u,v) -> Q14 NORMALIZED camera coordinates given principal point and focal length in pixels.
712// This is the hand-off contract to the geometry back end, which works entirely in normalized coords.
713func pf_normalize(u: i64, v: i64, cx: i64, cy: i64, f: i64, out2: *i64) -> i64 {
  if f == 0 { out2[0] = 0; out2[1] = 0; return 0 }
  out2[0] = ((u - cx) * PF_Q) / f
  out2[1] = ((v - cy) * PF_Q) / f
  return 0
718}
719
720// ---- 5. VERIFY: RANSAC whose model is the ESSENTIAL MATRIX. c0/c1 are Q14 normalized correspondence
721// arrays (stride 2). inl receives 1/0 per correspondence. outE receives the E re-fit on all inliers.
722// Returns the inlier count. Deterministic: the sample indices come from the fixed-seed LCG, so the same
723// input always yields the same E -- no random seed to explain away a bad run. ----
724func pf_ransac_e(c0: *i64, c1: *i64, n: i64, tol: i64, inl: *i64, outE: *i64) -> i64 {
  pf_init()
  if n < PF_RANSAC_SAMPLE { return 0 }
  let s0: *i64 = sys_mmap(PF_RANSAC_SAMPLE * 2 * 8 + 64) as *i64
  let s1: *i64 = sys_mmap(PF_RANSAC_SAMPLE * 2 * 8 + 64) as *i64
  let idx: *i64 = sys_mmap(PF_RANSAC_SAMPLE * 8 + 64) as *i64
  let etmp: *i64 = sys_mmap(9 * 8 + 64) as *i64
  let ebest: *i64 = sys_mmap(9 * 8 + 64) as *i64
  var bestcnt: i64 = 0
  var st: i64 = PF_LCG_SEED
  var it: i64 = 0
  while it < pf_c_ransac_iters {
      // draw PF_RANSAC_SAMPLE distinct indices
      var picked: i64 = 0
      var guard: i64 = 0
      while picked < PF_RANSAC_SAMPLE {
          if guard > PF_MAGIC_4096 {
              idx[picked] = picked % n
              picked = picked + 1
          } else {
              st = (st * PF_LCG_MUL + PF_LCG_ADD) & PF_LCG_MASK
              let cand: i64 = st % n
              var dup: i64 = 0
              var q: i64 = 0
              while q < picked {
                  if idx[q] == cand { dup = 1 }
                  q = q + 1
              }
              if dup == 0 {
                  idx[picked] = cand
                  picked = picked + 1
              }
              guard = guard + 1
          }
      }
      var k: i64 = 0
      while k < PF_RANSAC_SAMPLE {
          let ix: i64 = idx[k]
          s0[k * 2] = c0[ix * 2]
          s0[k * 2 + 1] = c0[ix * 2 + 1]
          s1[k * 2] = c1[ix * 2]
          s1[k * 2 + 1] = c1[ix * 2 + 1]
          k = k + 1
      }
      pf_estimate_e(s0, s1, PF_RANSAC_SAMPLE, etmp)
      var cnt: i64 = 0
      var i: i64 = 0
      while i < n {
          let r: i64 = e3_residual(etmp, c0[i * 2], c0[i * 2 + 1], c1[i * 2], c1[i * 2 + 1])
          if r <= tol { cnt = cnt + 1 }
          i = i + 1
      }
      if cnt > bestcnt {
          bestcnt = cnt
          var m: i64 = 0
          while m < 9 { ebest[m] = etmp[m]; m = m + 1 }
      }
      it = it + 1
  }
  if bestcnt == 0 { return 0 }
  // mark inliers under the winning model
  var i2: i64 = 0
  while i2 < n {
      let r2: i64 = e3_residual(ebest, c0[i2 * 2], c0[i2 * 2 + 1], c1[i2 * 2], c1[i2 * 2 + 1])
      if r2 <= tol { inl[i2] = 1 } else { inl[i2] = 0 }
      i2 = i2 + 1
  }
  // RE-FIT on every inlier: the minimal-sample E is only a hypothesis; the reported model should use
  // all the evidence that supports it (this is the step that turns RANSAC into an estimator).
  var nin: i64 = 0
  var c: i64 = 0
  while c < n { if inl[c] == 1 { nin = nin + 1 } c = c + 1 }
  if nin >= PF_RANSAC_SAMPLE {
      let f0: *i64 = sys_mmap(nin * 2 * 8 + 64) as *i64
      let f1: *i64 = sys_mmap(nin * 2 * 8 + 64) as *i64
      var w2: i64 = 0
      var i3: i64 = 0
      while i3 < n {
          if inl[i3] == 1 {
              f0[w2 * 2] = c0[i3 * 2]
              f0[w2 * 2 + 1] = c0[i3 * 2 + 1]
              f1[w2 * 2] = c1[i3 * 2]
              f1[w2 * 2 + 1] = c1[i3 * 2 + 1]
              w2 = w2 + 1
          }
          i3 = i3 + 1
      }
      pf_estimate_e(f0, f1, nin, outE)
  } else {
      var m2: i64 = 0
      while m2 < 9 { outE[m2] = ebest[m2]; m2 = m2 + 1 }
  }
  return bestcnt
817}
nx_photogram_front.nx source

`nx_photogram_front.nx` source
