nx_geom.nx source

1// nx_geom.nx -- geometric primitives (V4 of vision rollout). 2// 3// Hough transform line detection (Hough 1962, Duda+Hart 1972) and 4// circle detection (Ballard 1981 gradient method). Pure i64. 5// 6// What this unlocks: 7// + room geometry (walls, doors, window frames as line peaks) 8// + architectural analysis (vanishing points via line intersection) 9// + vehicle/animal anatomy (wheels, eyes as circles) 10// + table-top objects (plates, cups, bottles as circles) 11// + rule-of-thirds + leading-line aesthetics (line-angle histogram) 12// 13// All trig is Q10 fixed-point. cos/sin tables hardcoded for 36 angles 14// (5 degree resolution) -- sufficient for major-line detection; 15// finer-grained vision would interpolate or use larger tables. 16// 17// genealogy_id: hough_1962_patent + duda_hart_1972_generalized 18// + ballard_1981_generalized_circle 19// lineage_id: voting_accumulator + gradient_direction_voting 20// axioms: NX_AX_GEOM_LINE_HOMOGENEOUS (rho = x cos t + y sin t) 21 22// nx_safety_envelope: 23// intended_use: AUTO_APPLIED -- primitive-specific tuning queued 24// sil_target: SIL1 25// evidence: [bulk_applied_2026-05-16, see-file-comment-for-detail] 26// verdict: NOT_YET_EVALUATED 27 28import "syscalls.nx" 29import "nx_axioms.nx" 30import "nx_image.nx" 31import "nx_math.nx" 32 33// 36 angle bins at 5 degree resolution. 34const NX_GEOM_NTHETAS: i64 = 36 35const NX_GEOM_Q: i64 = 1024 36 37// Integer sqrt is canonical in nx_math.nx (Newton's method). 38 39// ===== Trig tables (Q10) ================================================ 40// 41// cos_t[i] = round(1024 * cos(i * 5 degrees)) for i in 0..35 42// sin_t[i] = round(1024 * sin(i * 5 degrees)) for i in 0..35 43// 44// These are recomputed per call (cheap, 36 stores) to keep the 45// primitive stateless. 46 47func nx_geom_fill_trig(cos_t: *i64, sin_t: *i64) -> i64 { 48 cos_t[0] = 1024 49 cos_t[1] = 1020 50 cos_t[2] = 1008 51 cos_t[3] = 989 52 cos_t[4] = 962 53 cos_t[5] = 928 54 cos_t[6] = 887 55 cos_t[7] = 839 56 cos_t[8] = 784 57 cos_t[9] = 724 58 cos_t[10] = 658 59 cos_t[11] = 587 60 cos_t[12] = 512 61 cos_t[13] = 433 62 cos_t[14] = 350 63 cos_t[15] = 265 64 cos_t[16] = 178 65 cos_t[17] = 89 66 cos_t[18] = 0 67 cos_t[19] = -89 68 cos_t[20] = -178 69 cos_t[21] = -265 70 cos_t[22] = -350 71 cos_t[23] = -433 72 cos_t[24] = -512 73 cos_t[25] = -587 74 cos_t[26] = -658 75 cos_t[27] = -724 76 cos_t[28] = -784 77 cos_t[29] = -839 78 cos_t[30] = -887 79 cos_t[31] = -928 80 cos_t[32] = -962 81 cos_t[33] = -989 82 cos_t[34] = -1008 83 cos_t[35] = -1020 84 85 sin_t[0] = 0 86 sin_t[1] = 89 87 sin_t[2] = 178 88 sin_t[3] = 265 89 sin_t[4] = 350 90 sin_t[5] = 433 91 sin_t[6] = 512 92 sin_t[7] = 587 93 sin_t[8] = 658 94 sin_t[9] = 724 95 sin_t[10] = 784 96 sin_t[11] = 839 97 sin_t[12] = 887 98 sin_t[13] = 928 99 sin_t[14] = 962 100 sin_t[15] = 989 101 sin_t[16] = 1008 102 sin_t[17] = 1020 103 sin_t[18] = 1024 104 sin_t[19] = 1020 105 sin_t[20] = 1008 106 sin_t[21] = 989 107 sin_t[22] = 962 108 sin_t[23] = 928 109 sin_t[24] = 887 110 sin_t[25] = 839 111 sin_t[26] = 784 112 sin_t[27] = 724 113 sin_t[28] = 658 114 sin_t[29] = 587 115 sin_t[30] = 512 116 sin_t[31] = 433 117 sin_t[32] = 350 118 sin_t[33] = 265 119 sin_t[34] = 178 120 sin_t[35] = 89 121 return 0 122} 123 124// ===== Hough line transform ============================================= 125// 126// Polar parametrization: rho = x*cos(theta) + y*sin(theta). 127// For theta in [0, 180) and rho in [-diag, +diag], each line maps to 128// a unique (rho, theta) point. Each edge pixel votes for all (theta, 129// rho) pairs consistent with the polar line equation. 130// 131// Returns up to max_lines lines satisfying votes >= min_votes and 132// being a 3x3 local maximum in the accumulator. 133 134struct HoughLines { 135 n_lines: i64, 136 rhos: *i64, // length n_lines 137 thetas: *i64, // degrees (0..175 step 5) 138 votes: *i64, 139} 140 141func nx_geom_hough_line(edge_mask: *Image, min_votes: i64, 142 max_lines: i64) -> *HoughLines { 143 let cos_t: *i64 = (sys_mmap(36 * 8)) as *i64 144 let sin_t: *i64 = (sys_mmap(36 * 8)) as *i64 145 nx_geom_fill_trig(cos_t, sin_t) 146 147 let w: i64 = edge_mask.width 148 let h: i64 = edge_mask.height 149 let d2: i64 = w * w + h * h 150 var diag: i64 = nx_math_isqrt(d2) + 1 151 let n_rhos: i64 = 2 * diag + 1 152 let acc_size: i64 = 36 * n_rhos 153 let acc: *i64 = (sys_mmap(acc_size * 8)) as *i64 154 var zi: i64 = 0 155 while zi < acc_size { 156 acc[zi] = 0 157 zi = zi + 1 158 } 159 160 // Accumulate votes. 161 var y: i64 = 0 162 while y < h { 163 var x: i64 = 0 164 while x < w { 165 if nx_image_get(edge_mask, x, y, 0) > 0 { 166 var t: i64 = 0 167 while t < 36 { 168 let rho_q: i64 = x * cos_t[t] + y * sin_t[t] 169 var rho: i64 = rho_q / 1024 170 let idx: i64 = t * n_rhos + (rho + diag) 171 if idx >= 0 { 172 if idx < acc_size { 173 acc[idx] = acc[idx] + 1 174 } 175 } 176 t = t + 1 177 } 178 } 179 x = x + 1 180 } 181 y = y + 1 182 } 183 184 let lines: *HoughLines = (sys_mmap(32)) as *HoughLines 185 lines.rhos = (sys_mmap(max_lines * 8)) as *i64 186 lines.thetas = (sys_mmap(max_lines * 8)) as *i64 187 lines.votes = (sys_mmap(max_lines * 8)) as *i64 188 lines.n_lines = 0 189 190 // Find peaks (3x3 local maximum, >= min_votes). 191 var t2: i64 = 0 192 while t2 < 36 { 193 var ri: i64 = 0 194 while ri < n_rhos { 195 let v: i64 = acc[t2 * n_rhos + ri] 196 if v >= min_votes { 197 var is_peak: i64 = 1 198 var dt: i64 = -1 199 while dt <= 1 { 200 var dr: i64 = -1 201 while dr <= 1 { 202 var is_center: i64 = 0 203 if dt == 0 { if dr == 0 { is_center = 1 } } 204 if is_center == 0 { 205 let nt: i64 = t2 + dt 206 let nr: i64 = ri + dr 207 if nt >= 0 { 208 if nt < 36 { 209 if nr >= 0 { 210 if nr < n_rhos { 211 if acc[nt * n_rhos + nr] > v { 212 is_peak = 0 213 } 214 } 215 } 216 } 217 } 218 } 219 dr = dr + 1 220 } 221 dt = dt + 1 222 } 223 if is_peak == 1 { 224 if lines.n_lines < max_lines { 225 lines.thetas[lines.n_lines] = t2 * 5 226 lines.rhos[lines.n_lines] = ri - diag 227 lines.votes[lines.n_lines] = v 228 lines.n_lines = lines.n_lines + 1 229 } 230 } 231 } 232 ri = ri + 1 233 } 234 t2 = t2 + 1 235 } 236 return lines 237} 238 239// ===== Hough circle (fixed radius, gradient-direction voting) =========== 240// 241// Ballard 1981 trick: an edge pixel at (x, y) with gradient (gx, gy) 242// lies on a circle of radius r whose center is at distance r along 243// the gradient direction (or opposite). This reduces the 3D Hough 244// accumulator to 2D, giving huge speedup. 245// 246// Inputs: 247// edge: binary mask 248// gx, gy: gradient images (from nx_image_sobel_x/y) 249// r: known radius to search for 250// min_votes: peak threshold 251// max_circles: cap on returned circles 252// 253// Returns up to max_circles centers (cx, cy) at radius r. 254 255struct HoughCircles { 256 n_circles: i64, 257 cxs: *i64, 258 cys: *i64, 259 rs: *i64, 260 votes: *i64, 261} 262 263func nx_geom_hough_circle_r(edge: *Image, gx: *ImageS64, gy: *ImageS64, 264 r: i64, min_votes: i64, 265 max_circles: i64) -> *HoughCircles { 266 let w: i64 = edge.width 267 let h: i64 = edge.height 268 let acc_n: i64 = w * h 269 let acc: *i64 = (sys_mmap(acc_n * 8)) as *i64 270 var zi: i64 = 0 271 while zi < acc_n { 272 acc[zi] = 0 273 zi = zi + 1 274 } 275 276 // Vote for centers along the gradient direction. 277 var y: i64 = 0 278 while y < h { 279 var x: i64 = 0 280 while x < w { 281 if nx_image_get(edge, x, y, 0) > 0 { 282 let gxv: i64 = nx_image_s64_get(gx, x, y) 283 let gyv: i64 = nx_image_s64_get(gy, x, y) 284 let mag2: i64 = gxv * gxv + gyv * gyv 285 if mag2 > 0 { 286 let mag: i64 = nx_math_isqrt(mag2) 287 if mag > 0 { 288 let dx: i64 = (r * gxv) / mag 289 let dy: i64 = (r * gyv) / mag 290 let cx1: i64 = x + dx 291 let cy1: i64 = y + dy 292 let cx2: i64 = x - dx 293 let cy2: i64 = y - dy 294 if cx1 >= 0 { 295 if cx1 < w { 296 if cy1 >= 0 { 297 if cy1 < h { 298 let i1: i64 = cy1 * w + cx1 299 acc[i1] = acc[i1] + 1 300 } 301 } 302 } 303 } 304 if cx2 >= 0 { 305 if cx2 < w { 306 if cy2 >= 0 { 307 if cy2 < h { 308 let i2: i64 = cy2 * w + cx2 309 acc[i2] = acc[i2] + 1 310 } 311 } 312 } 313 } 314 } 315 } 316 } 317 x = x + 1 318 } 319 y = y + 1 320 } 321 322 let circles: *HoughCircles = (sys_mmap(40)) as *HoughCircles 323 circles.cxs = (sys_mmap(max_circles * 8)) as *i64 324 circles.cys = (sys_mmap(max_circles * 8)) as *i64 325 circles.rs = (sys_mmap(max_circles * 8)) as *i64 326 circles.votes = (sys_mmap(max_circles * 8)) as *i64 327 circles.n_circles = 0 328 329 var py: i64 = 0 330 while py < h { 331 var px: i64 = 0 332 while px < w { 333 let v: i64 = acc[py * w + px] 334 if v >= min_votes { 335 var is_peak: i64 = 1 336 var dy2: i64 = -1 337 while dy2 <= 1 { 338 var dx2: i64 = -1 339 while dx2 <= 1 { 340 var is_center: i64 = 0 341 if dx2 == 0 { if dy2 == 0 { is_center = 1 } } 342 if is_center == 0 { 343 let nx2: i64 = px + dx2 344 let ny2: i64 = py + dy2 345 if nx2 >= 0 { 346 if nx2 < w { 347 if ny2 >= 0 { 348 if ny2 < h { 349 if acc[ny2 * w + nx2] > v { 350 is_peak = 0 351 } 352 } 353 } 354 } 355 } 356 } 357 dx2 = dx2 + 1 358 } 359 dy2 = dy2 + 1 360 } 361 if is_peak == 1 { 362 if circles.n_circles < max_circles { 363 circles.cxs[circles.n_circles] = px 364 circles.cys[circles.n_circles] = py 365 circles.rs[circles.n_circles] = r 366 circles.votes[circles.n_circles] = v 367 circles.n_circles = circles.n_circles + 1 368 } 369 } 370 } 371 px = px + 1 372 } 373 py = py + 1 374 } 375 return circles 376} 377 378// ===== Line-orientation histogram ======================================= 379// 380// Aggregate Hough line votes by orientation bin. Useful for: 381// - room detection (peaks at 0 and 90 degrees -> rectilinear room) 382// - diagonal composition (peak at 45 -> dynamic composition) 383// - radial composition (broad distribution -> radial layout) 384// 385// Returns histogram[36] of total votes per angle bin. 386 387func nx_geom_line_orientation_hist(lines: *HoughLines, hist: *i64) -> i64 { 388 var i: i64 = 0 389 while i < 36 { 390 hist[i] = 0 391 i = i + 1 392 } 393 var li: i64 = 0 394 while li < lines.n_lines { 395 let theta: i64 = lines.thetas[li] 396 let bin: i64 = theta / 5 397 if bin >= 0 { 398 if bin < 36 { 399 hist[bin] = hist[bin] + lines.votes[li] 400 } 401 } 402 li = li + 1 403 } 404 return 0 405}