nx_imgkp.nx source
↩ module page · 408 lines · 18832 B
1// nx_imgkp.nx -- LOCAL FEATURES: FAST corners + BRIEF binary descriptors, integer, deterministic.
2//
3// The global descriptors (dHash, edge-orientation layout) hit a hard ceiling the ruler MEASURED:
4// crop-30pct = 562, letterbox = 656 permille. That ceiling is structural -- a global fingerprint
5// summarises the WHOLE frame, so removing a third of the pixels or padding a border changes it no
6// matter how the summary is computed. TinEye clears exactly this by matching LOCAL regions and then
7// checking they agree GEOMETRICALLY: even if half the picture is gone, the surviving corners still
8// line up under one consistent shift.
9//
10// Every piece here is integer + zero-float + deterministic, which is the whole reason it belongs in
11// the sovereign engine rather than an OpenCV shell-out:
12// * FAST-9 corner test (Rosten & Drummond 2006) -- pure intensity comparisons on a Bresenham circle.
13// * BRIEF-256 descriptor (Calonder 2010) -- 256 intensity-pair comparisons -> a 256-bit string whose
14// distance is popcount Hamming, i.e. the nx_simhash kernel the engine already rides.
15// * A FIXED sampling pattern regenerated from one seed every call, so index-time and query-time
16// descriptors are bit-identical without any stored state.
17//
18// HONEST SCOPE, stated so it is not mistaken for more: plain BRIEF is NOT rotation-invariant. That is
19// deliberate and correct here -- the orientation axis (mirror / quarter-turn) is already owned by the
20// dihedral tier at the global level, and the two classes this rung targets (crop, letterbox) are pure
21// translations. Oriented-BRIEF (ORB) to combine rotation WITH occlusion is a later refinement, filed.
22//
23// PACKED LAYOUT (one image's feature set as a flat i64 block, so it fits the tier vtable's fixed
24// descriptor slot): out[0] = keypoint count; then KP_MAX records of 6 i64 each = [x, y, d0,d1,d2,d3].
25// genealogy_id: rosten_2006_fast + calonder_2010_brief. license_tier: ORIGINAL
26import "syscalls.nx"
27import "nx_simhash.nx"
28
29const KP_MAX: i64 = 24 // keypoints kept per image (spatially spread by NMS)
30const KP_WORDS: i64 = 4 // 256-bit BRIEF descriptor = 4 x i64
31const KP_REC: i64 = 6 // per keypoint: x, y, 4 descriptor words
32const KP_PACK: i64 = 145 // 1 + KP_MAX*KP_REC -- the tier idx_dim (UNCHANGED by the pyramid)
33const KP_FAST_N: i64 = 9 // contiguous arc length for a FAST-9 corner
34const KP_FAST_T: i64 = 16 // intensity threshold (of 0..255)
35const KP_PATCH: i64 = 6 // BRIEF samples drawn within +/- this of the keypoint
36const KP_NMS_R: i64 = 3 // non-max suppression radius (spread keypoints out)
37const KP_CAND_PEROCT: i64 = 8192 // per-octave candidate cap
38const KP_CAND_MAX: i64 = 8192 // total candidate cap = KP_CAND_PEROCT * KP_OCTAVES (bounds per-call memory)
39const KP_BRIEF_BITS: i64 = 256
40// SCALE INVARIANCE -- SOLVED ON THE QUERY SIDE, NOT THE INDEX SIDE (2026-08-05).
41//
42// HISTORY, MEASURED 2026-07-24: a NAIVE 3-octave INDEX pyramid REGRESSED the working same-scale
43// classes (crop-30 937->593, crop-20 1000->906) and still did not recover crop+rescale (500->437).
44// TWO mechanisms, both structural to putting the pyramid in the index:
45// 1. BUDGET COMPETITION -- the fixed KP_MAX records are shared, so every octave-1/2 keypoint
46// admitted EVICTS an octave-0 keypoint that the same-scale classes depend on.
47// 2. CROSS-OCTAVE SUPPRESSION -- selection NMS runs in BASE coordinates, so a coarse-octave
48// keypoint silently suppresses fine-octave neighbours within KP_NMS_R. Real ORB suppresses
49// WITHIN a level, never across them.
50// Both are competitions for one budget, so no amount of tuning removes them: the index descriptor
51// is a fixed-size summary and scale is one more thing demanding room in it.
52//
53// THE FIX IS TO MOVE THE SCALE WORK TO THE SIDE THAT HAS NO BUDGET PRESSURE. The tier vtable already
54// separates idx_dim from qry_dim, and orient-d4 is the precedent: it stores ONE hash and computes the
55// query's EIGHT-element D4 orbit, taking the best. Scale is the same shape of problem. So the INDEX
56// stays single-scale and bit-identical (no reindex, no storage growth over 10^8 images), and the
57// QUERY is described at KP_QSCALES scales, each getting its OWN FULL KP_MAX budget -- no competition,
58// no cross-octave NMS, because each scale is an independent detector run on its own image.
59//
60// WHY THIS CANNOT REGRESS, structurally rather than hopefully: scale 1/1 is IN the sweep and is
61// byte-identical to today's descriptor, and the tier takes the BEST (minimum) distance over the
62// sweep. A minimum that includes today's value can never exceed today's value.
63//
64// The ladder is half-octave: 1/2, 7/10, 1/1, 10/7, 2/1 (2^-1 .. 2^+1). It is two-sided because the
65// query may be either larger or smaller than the corpus image. 7/10 is exact for the ruler's
66// crop+rescale class (centre-70% crop resampled back = content at 1/0.7), and 10/7 is its inverse.
67const KP_OCTAVES: i64 = 1 // INDEX stays single-scale, deliberately -- see above
68const KP_QSCALES: i64 = 5 // query-side scale sweep width
69const KP_QPIXCAP: i64 = 40000000 // skip a sweep scale whose resampled image exceeds this
70
71// sweep scale i as a rational num/den. Explicit if-chains, NOT a const array: const-array indexing
72// has a known miscompile trap in this toolchain (same reason kp_circle is written out).
73func kp_qscale_num(i: i64) -> i64 {
74 if i == 0 { return 1 }
75 if i == 1 { return 7 }
76 if i == 2 { return 10 }
77 if i == 3 { return 1 }
78 return 2
79}
80func kp_qscale_den(i: i64) -> i64 {
81 if i == 0 { return 1 }
82 if i == 1 { return 10 }
83 if i == 2 { return 7 }
84 if i == 3 { return 2 }
85 return 1
86}
87
88func kp_pack_dim() -> i64 { return KP_PACK }
89
90// the 16 Bresenham-circle offsets (radius 3), written into caller buffers. Kept as explicit stores
91// rather than a const array -- const-array indexing has a known miscompile trap in this toolchain.
92func kp_circle(cx: *i64, cy: *i64) -> i64 {
93 cx[0]=0; cy[0]=0-3
94 cx[1]=1; cy[1]=0-3
95 cx[2]=2; cy[2]=0-2
96 cx[3]=3; cy[3]=0-1
97 cx[4]=3; cy[4]=0
98 cx[5]=3; cy[5]=1
99 cx[6]=2; cy[6]=2
100 cx[7]=1; cy[7]=3
101 cx[8]=0; cy[8]=3
102 cx[9]=0-1; cy[9]=3
103 cx[10]=0-2; cy[10]=2
104 cx[11]=0-3; cy[11]=1
105 cx[12]=0-3; cy[12]=0
106 cx[13]=0-3; cy[13]=0-1
107 cx[14]=0-2; cy[14]=0-2
108 cx[15]=0-1; cy[15]=0-3
109 return 16
110}
111
112func kp_rng(s: *i64) -> i64 { let x: i64 = s[0] * 0x5851F42D4C957F2D + 0x14057B7EF767814F; s[0] = x; return x }
113
114// fill the BRIEF sampling pattern: 256 pairs of (ax,ay,bx,by), each coordinate in [-KP_PATCH,KP_PATCH].
115// FIXED seed -> identical pattern on every call -> descriptors are comparable across images with no
116// stored state. pat holds 256*4 ints.
117func kp_brief_pattern(pat: *i64) -> i64 {
118 let s: *i64 = sys_mmap(8) as *i64
119 s[0] = 0x1D872B41A6F3C5E9
120 let span: i64 = 2 * KP_PATCH + 1
121 var i: i64 = 0
122 while i < KP_BRIEF_BITS {
123 var r: i64 = kp_rng(s) >> 11
124 if r < 0 { r = 0 - r }
125 pat[i*4] = (r % span) - KP_PATCH
126 r = kp_rng(s) >> 11; if r < 0 { r = 0 - r }
127 pat[i*4+1] = (r % span) - KP_PATCH
128 r = kp_rng(s) >> 11; if r < 0 { r = 0 - r }
129 pat[i*4+2] = (r % span) - KP_PATCH
130 r = kp_rng(s) >> 11; if r < 0 { r = 0 - r }
131 pat[i*4+3] = (r % span) - KP_PATCH
132 i = i + 1
133 }
134 return 0
135}
136
137func kp_abs(v: i64) -> i64 { if v < 0 { return 0 - v } return v }
138
139// FAST-9 corner test at (x,y): is there an arc of >=9 contiguous circle pixels all brighter than
140// Ip+t or all darker than Ip-t? Writes the corner SCORE (sum |Ii-Ip| over the circle) to score[0].
141// bright/dark are CALLER-OWNED 16-slot scratch buffers -- hoisted out so this per-pixel hot path makes
142// ZERO syscalls (an earlier version sys_mmap'd them every call: millions of syscalls, the real cost).
143func kp_is_corner(gray: *u8, w: i64, h: i64, x: i64, y: i64, cx: *i64, cy: *i64, bright: *i64, dark: *i64, score: *i64) -> i64 {
144 let ip: i64 = gray[y*w+x] as i64
145 var sc: i64 = 0
146 var i: i64 = 0
147 while i < 16 {
148 let ii: i64 = gray[(y+cy[i])*w + (x+cx[i])] as i64
149 if ii > ip + KP_FAST_T { bright[i] = 1 } else { bright[i] = 0 }
150 if ii < ip - KP_FAST_T { dark[i] = 1 } else { dark[i] = 0 }
151 sc = sc + kp_abs(ii - ip)
152 i = i + 1
153 }
154 score[0] = sc
155 // contiguous run of >=9 over the circular 16 (scan 16+8 with wraparound)
156 var runb: i64 = 0
157 var rund: i64 = 0
158 var corner: i64 = 0
159 i = 0
160 while i < 16 + KP_FAST_N - 1 {
161 let idx: i64 = i % 16
162 if bright[idx] == 1 { runb = runb + 1 } else { runb = 0 }
163 if dark[idx] == 1 { rund = rund + 1 } else { rund = 0 }
164 if runb >= KP_FAST_N { corner = 1 }
165 if rund >= KP_FAST_N { corner = 1 }
166 i = i + 1
167 }
168 return corner
169}
170
171// compute the 256-bit BRIEF descriptor at (x,y) into desc[0..4). bit set if I(p+a) < I(p+b).
172func kp_brief_at(gray: *u8, w: i64, h: i64, x: i64, y: i64, pat: *i64, desc: *i64) -> i64 {
173 desc[0] = 0; desc[1] = 0; desc[2] = 0; desc[3] = 0
174 var b: i64 = 0
175 while b < KP_BRIEF_BITS {
176 let ia: i64 = gray[(y+pat[b*4+1])*w + (x+pat[b*4])] as i64
177 let ib: i64 = gray[(y+pat[b*4+3])*w + (x+pat[b*4+2])] as i64
178 if ia < ib {
179 let word: i64 = b >> 6
180 let bit: i64 = b & 63
181 desc[word] = desc[word] | (1 << bit)
182 }
183 b = b + 1
184 }
185 return 0
186}
187
188// nearest-neighbour downscale of a gray image into a caller buffer (the pyramid builder).
189func kp_downscale(src: *u8, sw: i64, sh: i64, dst: *u8, dw: i64, dh: i64) -> i64 {
190 var oy: i64 = 0
191 while oy < dh {
192 let sy: i64 = (oy * sh) / dh
193 var ox: i64 = 0
194 while ox < dw {
195 let sx: i64 = (ox * sw) / dw
196 dst[oy * dw + ox] = src[sy * sw + sx]
197 ox = ox + 1
198 }
199 oy = oy + 1
200 }
201 return 0
202}
203
204// target width/height of octave o for a base w x h: o=0 full, o=1 /1.4 (5/7), o=2 /2. Geometric.
205func kp_oct_dim(o: i64, base: i64) -> i64 {
206 if o == 0 { return base }
207 if o == 1 { return base * 5 / 7 }
208 return base / 2
209}
210
211// DETECT + DESCRIBE across a SCALE PYRAMID: fill out[0..KP_PACK) = [count, (x,y,d0..d3) x KP_MAX].
212// FAST corners are found at KP_OCTAVES scales; a keypoint's coordinates are stored in BASE resolution
213// but its BRIEF descriptor is computed on the octave image where it was found, so a zoomed copy of the
214// image still yields matching descriptors at the octave whose patch covers the same physical content.
215// Top-KP_MAX by corner score after base-coordinate NMS. Zero-fills unused records. Returns the count.
216func nx_imgkp_describe(gray: *u8, w: i64, h: i64, out: *i64) -> i64 {
217 var i: i64 = 0
218 while i < KP_PACK { out[i] = 0; i = i + 1 }
219
220 let margin: i64 = 3 + KP_PATCH // room for the circle AND the BRIEF patch (per octave)
221
222 let cx: *i64 = sys_mmap(8*16) as *i64
223 let cy: *i64 = sys_mmap(8*16) as *i64
224 kp_circle(cx, cy)
225 let bright: *i64 = sys_mmap(8*16) as *i64 // hoisted scratch for kp_is_corner (allocated ONCE)
226 let dark: *i64 = sys_mmap(8*16) as *i64
227
228 // candidate accumulators: base coords + score + octave + octave-local coords (to describe later)
229 let candbx: *i64 = sys_mmap(8*KP_CAND_MAX) as *i64
230 let candby: *i64 = sys_mmap(8*KP_CAND_MAX) as *i64
231 let cands: *i64 = sys_mmap(8*KP_CAND_MAX) as *i64
232 let cand_oct: *i64 = sys_mmap(8*KP_CAND_MAX) as *i64
233 let cand_ox: *i64 = sys_mmap(8*KP_CAND_MAX) as *i64
234 let cand_oy: *i64 = sys_mmap(8*KP_CAND_MAX) as *i64
235 let scorebox: *i64 = sys_mmap(8) as *i64
236 var ncand: i64 = 0
237
238 // per-octave image pointers + dims (octave 0 = the input itself, no copy)
239 let oimg: *i64 = sys_mmap(8*KP_OCTAVES) as *i64
240 let otw: *i64 = sys_mmap(8*KP_OCTAVES) as *i64
241 let oth: *i64 = sys_mmap(8*KP_OCTAVES) as *i64
242
243 var o: i64 = 0
244 while o < KP_OCTAVES {
245 var tw: i64 = kp_oct_dim(o, w)
246 var th: i64 = kp_oct_dim(o, h)
247 var og: *u8 = gray
248 if o != 0 {
249 if tw < 1 { tw = 1 }
250 if th < 1 { th = 1 }
251 og = sys_mmap(tw * th)
252 kp_downscale(gray, w, h, og, tw, th)
253 }
254 oimg[o] = og as i64; otw[o] = tw; oth[o] = th
255 var octn: i64 = 0 // per-octave count so a dense octave-0 can't starve octaves 1..2
256 if tw >= 2*margin + 1 { if th >= 2*margin + 1 {
257 var y: i64 = margin
258 while y < th - margin {
259 var x: i64 = margin
260 while x < tw - margin {
261 if ncand < KP_CAND_MAX { if octn < KP_CAND_PEROCT {
262 if kp_is_corner(og, tw, th, x, y, cx, cy, bright, dark, scorebox) == 1 {
263 candbx[ncand] = x * w / tw // map octave coord -> base resolution
264 candby[ncand] = y * h / th
265 cands[ncand] = scorebox[0]
266 cand_oct[ncand] = o; cand_ox[ncand] = x; cand_oy[ncand] = y
267 ncand = ncand + 1; octn = octn + 1
268 }
269 } }
270 x = x + 1
271 }
272 y = y + 1
273 }
274 } }
275 o = o + 1
276 }
277
278 let pat: *i64 = sys_mmap(8 * KP_BRIEF_BITS * 4) as *i64
279 kp_brief_pattern(pat)
280 let taken: *u8 = sys_mmap(KP_CAND_MAX)
281 i = 0
282 while i < ncand { taken[i] = 0 as u8; i = i + 1 }
283
284 // greedy top-KP_MAX by score with base-coordinate NMS (SEPARATE loop flag -- never overload nkp as
285 // the break signal; that once made a flat image report 24 all-zero keypoints, caught by the gate).
286 var nkp: i64 = 0
287 var go: i64 = 1
288 while go == 1 {
289 if nkp >= KP_MAX { go = 0 } else {
290 var best: i64 = 0 - 1
291 var bestsc: i64 = 0 - 1
292 var c: i64 = 0
293 while c < ncand {
294 if taken[c] == (0 as u8) {
295 if cands[c] > bestsc { bestsc = cands[c]; best = c }
296 }
297 c = c + 1
298 }
299 if best < 0 { go = 0 } else {
300 let bx: i64 = candbx[best]
301 let by: i64 = candby[best]
302 let rec: i64 = 1 + nkp * KP_REC
303 out[rec] = bx
304 out[rec+1] = by
305 let desc: *i64 = ((out as i64) + (rec+2)*8) as *i64
306 // describe on the OCTAVE image where the corner was found (scale-appropriate patch)
307 let oo: i64 = cand_oct[best]
308 kp_brief_at(oimg[oo] as *u8, otw[oo], oth[oo], cand_ox[best], cand_oy[best], pat, desc)
309 nkp = nkp + 1
310 c = 0
311 while c < ncand {
312 if taken[c] == (0 as u8) {
313 if kp_abs(candbx[c]-bx) <= KP_NMS_R { if kp_abs(candby[c]-by) <= KP_NMS_R { taken[c] = 1 as u8 } }
314 }
315 c = c + 1
316 }
317 }
318 }
319 }
320 out[0] = nkp
321 return nkp
322}
323
324func nx_imgkp_count(pk: *i64) -> i64 { return pk[0] }
325
326// Hamming distance between two packed keypoints' 256-bit descriptors (sum of per-word popcounts).
327func kp_desc_hamming(pk_a: *i64, ia: i64, pk_b: *i64, ib: i64) -> i64 {
328 let ba: i64 = 1 + ia*KP_REC + 2
329 let bb: i64 = 1 + ib*KP_REC + 2
330 var d: i64 = 0
331 var w: i64 = 0
332 while w < KP_WORDS { d = d + nx_simhash_hamming(pk_a[ba+w], pk_b[bb+w]); w = w + 1 }
333 return d
334}
335
336// MATCH the query keypoints to the corpus keypoints: for each query kp find its nearest corpus kp by
337// descriptor Hamming, keep it only if it passes Lowe's ratio test (best clearly better than second).
338// Fills the corresponding coordinate lists; returns the match count. This is what feeds RANSAC.
339const KP_MATCH_MAXHAM: i64 = 96 // absolute Hamming ceiling for a plausible match (of 256 bits)
340func nx_imgkp_match(q: *i64, d: *i64, out_qx: *i64, out_qy: *i64, out_dx: *i64, out_dy: *i64, cap: i64) -> i64 {
341 let nq: i64 = q[0]
342 let nd: i64 = d[0]
343 var m: i64 = 0
344 var i: i64 = 0
345 while i < nq {
346 var best: i64 = 999
347 var second: i64 = 999
348 var bestj: i64 = 0 - 1
349 var j: i64 = 0
350 while j < nd {
351 let hh: i64 = kp_desc_hamming(q, i, d, j)
352 if hh < best { second = best; best = hh; bestj = j } else { if hh < second { second = hh } }
353 j = j + 1
354 }
355 // ratio test: best/second < 0.8, i.e. best*5 < second*4 -- plus an absolute ceiling
356 if bestj >= 0 { if best <= KP_MATCH_MAXHAM { if best * 5 < second * 4 {
357 if m < cap {
358 out_qx[m] = q[1 + i*KP_REC]
359 out_qy[m] = q[1 + i*KP_REC + 1]
360 out_dx[m] = d[1 + bestj*KP_REC]
361 out_dy[m] = d[1 + bestj*KP_REC + 1]
362 m = m + 1
363 }
364 } } }
365 i = i + 1
366 }
367 return m
368}
369
370// ===== QUERY-SIDE SCALE SWEEP =====================================
371//
372// Fills KP_QSCALES consecutive packs: out[i*KP_PACK .. (i+1)*KP_PACK). Pack 0 is the identity scale
373// and is BYTE-IDENTICAL to nx_imgkp_describe on the same image -- that identity is what makes the
374// tier's minimum-over-sweep provably non-regressive.
375//
376// Keypoint coordinates stay in EACH SCALE'S OWN frame, deliberately: at the scale that undoes the
377// query's rescale, the query and the corpus image are at the same effective resolution, so a pure
378// TRANSLATION relates them and the cheap translation-RANSAC confirms the match. Mapping coordinates
379// back to base resolution would instead force every scale through the O(n^3) similarity search.
380func nx_imgkp_query_dim() -> i64 { return KP_PACK * KP_QSCALES }
381
382func nx_imgkp_describe_query(gray: *u8, w: i64, h: i64, out: *i64) -> i64 {
383 let minside: i64 = 2*(3 + KP_PATCH) + 1
384 var i: i64 = 0
385 while i < KP_QSCALES {
386 let dest: *i64 = ((out as i64) + i*KP_PACK*8) as *i64
387 var k: i64 = 0
388 while k < KP_PACK { dest[k] = 0; k = k + 1 }
389 let num: i64 = kp_qscale_num(i)
390 let den: i64 = kp_qscale_den(i)
391 if num == den {
392 nx_imgkp_describe(gray, w, h, dest)
393 } else {
394 let tw: i64 = w * num / den
395 let th: i64 = h * num / den
396 if tw >= minside { if th >= minside { if tw * th <= KP_QPIXCAP {
397 let sc: *u8 = sys_mmap(tw * th + 64)
398 kp_downscale(gray, w, h, sc, tw, th)
399 nx_imgkp_describe(sc, tw, th, dest)
400 sys_munmap(sc, tw * th + 64)
401 } } }
402 // a scale that is too small, too large, or unbuildable leaves count=0, which the tier
403 // skips -- an absent scale must never be a match, and never a crash
404 }
405 i = i + 1
406 }
407 return KP_QSCALES
408}