nx_image.nx source
↩ module page · 428 lines · 15132 B
1// nx_image.nx -- pure-math image primitives.
2//
3// No model files, no learned weights, no VLM. Classical computer
4// vision math: pixel access, grayscale conversion, 3x3 convolution,
5// Sobel gradient, gradient magnitude via integer sqrt, histogram.
6//
7// All i64. Uses nx_i128 muldiv where intermediate products risk
8// overflow (per NUMERIC_TIER_LADDER + S_CLASS_MATH_DOCTRINE).
9//
10// genealogy_id: gonzalez_woods_digital_image_processing + canny_1986
11// + sobel_feldman_1968
12// lineage_id: convolution + gradient + integer_pixel_arithmetic
13// axioms: NX_AX_ALG_DISTRIBUTIVITY, NX_AX_ALG_COMMUTATIVITY
14//
15// nx_safety_envelope:
16// intended_use: "Image container + convolution/gradient/
17// integer-pixel ops. Foundation for all
18// substrate visual processing."
19// sil_target: SIL2 (visual integrity in safety-relevant
20// UI -- e.g. medical imaging,
21// aviation HUDs)
22// asil_target: QM
23// dal_target: DAL C
24// evidence: [Image_canonical_container_type,
25// bounded_pixel_loops,
26// axiom_distributivity_commutativity_tagged]
27// hazard_register: [bug-tape-image-stride-not-validated,
28// bug-tape-OOB-pixel-read-via-bad-coord,
29// bug-tape-channel-count-mismatch-corruption]
30// residual_risk: "Image stride + dimensions are caller-
31// supplied. Substrate validates bounds at
32// function entry; loop bodies trust them.
33// Image format compatibility (RGB vs RGBA
34// vs YUV) is upstream responsibility."
35// verdict: NOT_YET_EVALUATED
36
37import "nx_syscalls.nx"
38import "nx_axioms.nx"
39import "nx_i128.nx"
40
41// ===== image struct ====================================================
42//
43// pixels = byte buffer, row-major. For RGB: 3 bytes per pixel.
44// stride = bytes per row (== width * channels for tightly packed).
45
46struct Image {
47 pixels: *u8,
48 width: i64,
49 height: i64,
50 channels: i64,
51 stride: i64,
52}
53
54const NX_IMG_BYTES: i64 = 40
55
56func nx_image_alloc(width: i64, height: i64, channels: i64) -> *Image {
57 let raw: *u8 = sys_mmap(NX_IMG_BYTES)
58 let img: *Image = raw as *Image
59 img.width = width
60 img.height = height
61 img.channels = channels
62 img.stride = width * channels
63 img.pixels = sys_mmap(width * height * channels + 16)
64 return img
65}
66
67// Allocate a 1-channel signed-grayscale image whose values can be
68// negative (e.g., gradient output). We store as i64 array under
69// the hood (width*height*8 bytes).
70struct ImageS64 {
71 data: *i64,
72 width: i64,
73 height: i64,
74}
75
76const NX_IMGS64_BYTES: i64 = 24
77
78func nx_image_s64_alloc(width: i64, height: i64) -> *ImageS64 {
79 let raw: *u8 = sys_mmap(NX_IMGS64_BYTES)
80 let img: *ImageS64 = raw as *ImageS64
81 img.width = width
82 img.height = height
83 img.data = (sys_mmap(width * height * 8 + 16)) as *i64
84 return img
85}
86
87// ===== pixel access (u8 image) =========================================
88
89func nx_image_get(img: *Image, x: i64, y: i64, c: i64) -> i64 {
90 if x < 0 { return 0 }
91 if y < 0 { return 0 }
92 if x >= img.width { return 0 }
93 if y >= img.height { return 0 }
94 return img.pixels[y * img.stride + x * img.channels + c]
95}
96
97func nx_image_set(img: *Image, x: i64, y: i64, c: i64, v: i64) -> i64 {
98 if x < 0 { return -1 }
99 if y < 0 { return -1 }
100 if x >= img.width { return -1 }
101 if y >= img.height { return -1 }
102 var vc: i64 = v
103 if vc < 0 { vc = 0 }
104 if vc > 255 { vc = 255 }
105 img.pixels[y * img.stride + x * img.channels + c] = vc
106 return 0
107}
108
109// ===== signed-grayscale access ========================================
110
111func nx_image_s64_get(img: *ImageS64, x: i64, y: i64) -> i64 {
112 if x < 0 { return 0 }
113 if y < 0 { return 0 }
114 if x >= img.width { return 0 }
115 if y >= img.height { return 0 }
116 return img.data[y * img.width + x]
117}
118
119func nx_image_s64_set(img: *ImageS64, x: i64, y: i64, v: i64) -> i64 {
120 if x < 0 { return -1 }
121 if y < 0 { return -1 }
122 if x >= img.width { return -1 }
123 if y >= img.height { return -1 }
124 img.data[y * img.width + x] = v
125 return 0
126}
127
128// ===== grayscale conversion (Rec. 709 weights, scaled) =================
129//
130// Y = (218 * R + 732 * G + 74 * B) / 1024
131// (approximates 0.2126 + 0.7152 + 0.0722 with 10-bit precision).
132
133func nx_image_to_grayscale(rgb: *Image) -> *Image {
134 let gray: *Image = nx_image_alloc(rgb.width, rgb.height, 1)
135 var y: i64 = 0
136 while y < rgb.height {
137 var x: i64 = 0
138 while x < rgb.width {
139 let r: i64 = nx_image_get(rgb, x, y, 0)
140 let g: i64 = nx_image_get(rgb, x, y, 1)
141 let b: i64 = nx_image_get(rgb, x, y, 2)
142 let lum: i64 = (218 * r + 732 * g + 74 * b) / 1024
143 nx_image_set(gray, x, y, 0, lum)
144 x = x + 1
145 }
146 y = y + 1
147 }
148 return gray
149}
150
151// ===== 3x3 convolution (kernel is i64[9], output as ImageS64) ==========
152//
153// Output: out[y, x] = sum over (kj, ki) of kernel[(kj+1)*3 + (ki+1)] *
154// input[y+kj, x+ki] for kj, ki in {-1, 0, 1}.
155// Border policy: clamp to edge (uses nx_image_get's bounds check).
156
157func nx_image_convolve_3x3_into(src: *Image, kernel: *i64, channel: i64, out: *ImageS64) -> *ImageS64 {
158 // 'out' is CALLER-OWNED scratch. The allocating wrapper is at the end of this file.
159 var y: i64 = 0
160 while y < src.height {
161 var x: i64 = 0
162 while x < src.width {
163 var sum: i64 = 0
164 var kj: i64 = 0
165 while kj < 3 {
166 var ki: i64 = 0
167 while ki < 3 {
168 let yy: i64 = y + kj - 1
169 let xx: i64 = x + ki - 1
170 let px: i64 = nx_image_get(src, xx, yy, channel)
171 sum = sum + kernel[kj * 3 + ki] * px
172 ki = ki + 1
173 }
174 kj = kj + 1
175 }
176 nx_image_s64_set(out, x, y, sum)
177 x = x + 1
178 }
179 y = y + 1
180 }
181 return out
182}
183
184// ===== Sobel gradient ==================================================
185//
186// Sx = [-1 0 1 ; -2 0 2 ; -1 0 1]
187// Sy = [-1 -2 -1 ; 0 0 0 ; 1 2 1]
188
189// The ALLOCATING form: exactly the contract this file had before the 2026-08-25 refactor, now one
190// line over the scratch variant. Every existing caller in the estate is unchanged BY CONSTRUCTION.
191func nx_image_convolve_3x3(src: *Image, kernel: *i64, channel: i64) -> *ImageS64 {
192 let out: *ImageS64 = nx_image_s64_alloc(src.width, src.height)
193 return nx_image_convolve_3x3_into(src, kernel, channel, out)
194}
195
196// Sobel into caller-owned scratch. The 72-byte kernel is STILL built per call and that is
197// deliberate, not an oversight: at 72 bytes it is four orders of magnitude below the w*h*8 buffer
198// this refactor exists to stop reallocating, so hoisting it would add a parameter to every call
199// site for no measurable gain. Naming the imprecision beats silently implying there is none.
200func nx_image_sobel_x_into(gray: *Image, out: *ImageS64) -> *ImageS64 {
201 let k: *i64 = (sys_mmap(72)) as *i64
202 k[0] = -1; k[1] = 0; k[2] = 1
203 k[3] = -2; k[4] = 0; k[5] = 2
204 k[6] = -1; k[7] = 0; k[8] = 1
205 return nx_image_convolve_3x3_into(gray, k, 0, out)
206}
207
208func nx_image_sobel_y_into(gray: *Image, out: *ImageS64) -> *ImageS64 {
209 let k: *i64 = (sys_mmap(72)) as *i64
210 k[0] = -1; k[1] = -2; k[2] = -1
211 k[3] = 0; k[4] = 0; k[5] = 0
212 k[6] = 1; k[7] = 2; k[8] = 1
213 return nx_image_convolve_3x3_into(gray, k, 0, out)
214}
215
216func nx_image_sobel_x(gray: *Image) -> *ImageS64 {
217 let k: *i64 = (sys_mmap(72)) as *i64
218 k[0] = -1; k[1] = 0; k[2] = 1
219 k[3] = -2; k[4] = 0; k[5] = 2
220 k[6] = -1; k[7] = 0; k[8] = 1
221 return nx_image_convolve_3x3(gray, k, 0)
222}
223
224func nx_image_sobel_y(gray: *Image) -> *ImageS64 {
225 let k: *i64 = (sys_mmap(72)) as *i64
226 k[0] = -1; k[1] = -2; k[2] = -1
227 k[3] = 0; k[4] = 0; k[5] = 0
228 k[6] = 1; k[7] = 2; k[8] = 1
229 return nx_image_convolve_3x3(gray, k, 0)
230}
231
232// Gradient magnitude: sqrt(gx^2 + gy^2) via nx_th_isqrt-style.
233// Uses muldiv to avoid overflow in gx^2 + gy^2.
234func nx_image_gradient_magnitude(gx: *ImageS64, gy: *ImageS64) -> *ImageS64 {
235 let out: *ImageS64 = nx_image_s64_alloc(gx.width, gx.height)
236 var y: i64 = 0
237 while y < gx.height {
238 var x: i64 = 0
239 while x < gx.width {
240 let gxv: i64 = nx_image_s64_get(gx, x, y)
241 let gyv: i64 = nx_image_s64_get(gy, x, y)
242 let s: i64 = nx_muldiv_i64(gxv, gxv, 1) + nx_muldiv_i64(gyv, gyv, 1)
243 // integer sqrt via Newton's method
244 var mag: i64 = 0
245 if s > 0 {
246 var z: i64 = s
247 var w: i64 = (z + 1) / 2
248 while w < z {
249 z = w
250 w = (z + s / z) / 2
251 }
252 mag = z
253 }
254 nx_image_s64_set(out, x, y, mag)
255 x = x + 1
256 }
257 y = y + 1
258 }
259 return out
260}
261
262// ===== histogram =======================================================
263//
264// 256-bin intensity histogram for a single-channel image. Output is
265// caller-allocated i64[256] filled with counts.
266
267func nx_image_histogram_256(img: *Image, channel: i64, out_bins: *i64) -> i64 {
268 var i: i64 = 0
269 while i < 256 { out_bins[i] = 0; i = i + 1 }
270 var y: i64 = 0
271 while y < img.height {
272 var x: i64 = 0
273 while x < img.width {
274 let v: i64 = nx_image_get(img, x, y, channel)
275 if v >= 0 {
276 if v < 256 {
277 out_bins[v] = out_bins[v] + 1
278 }
279 }
280 x = x + 1
281 }
282 y = y + 1
283 }
284 return 0
285}
286
287// Otsu's method: find the threshold that maximizes between-class
288// variance. Returns threshold in [0, 255]. Pure classical statistics.
289//
290// genealogy_id: otsu_1979
291// lineage_id: histogram + variance + maximization
292// axioms: NX_AX_PROB_NONNEGATIVITY
293
294func nx_image_otsu_threshold(img: *Image, channel: i64) -> i64 {
295 let bins: *i64 = (sys_mmap(256 * 8)) as *i64
296 nx_image_histogram_256(img, channel, bins)
297 var total: i64 = 0
298 var sum_all: i64 = 0
299 var i: i64 = 0
300 while i < 256 {
301 total = total + bins[i]
302 sum_all = sum_all + i * bins[i]
303 i = i + 1
304 }
305 if total == 0 { return 0 }
306 var sum_back: i64 = 0
307 var w_back: i64 = 0
308 var max_var: i64 = -1
309 var best_t: i64 = 0
310 var t: i64 = 0
311 while t < 256 {
312 w_back = w_back + bins[t]
313 if w_back > 0 {
314 let w_fore: i64 = total - w_back
315 if w_fore > 0 {
316 sum_back = sum_back + t * bins[t]
317 let mean_back: i64 = sum_back / w_back
318 let mean_fore: i64 = (sum_all - sum_back) / w_fore
319 let diff: i64 = mean_back - mean_fore
320 let var_between: i64 = nx_muldiv_i64(
321 nx_muldiv_i64(w_back, w_fore, 1),
322 nx_muldiv_i64(diff, diff, 1), 1)
323 if var_between > max_var {
324 max_var = var_between
325 best_t = t
326 }
327 }
328 }
329 t = t + 1
330 }
331 return best_t
332}
333
334// ===== threshold to binary mask =======================================
335
336func nx_image_threshold(img: *Image, channel: i64, thresh: i64) -> *Image {
337 let mask: *Image = nx_image_alloc(img.width, img.height, 1)
338 var y: i64 = 0
339 while y < img.height {
340 var x: i64 = 0
341 while x < img.width {
342 let v: i64 = nx_image_get(img, x, y, channel)
343 var out: i64 = 0
344 if v >= thresh { out = 255 }
345 nx_image_set(mask, x, y, 0, out)
346 x = x + 1
347 }
348 y = y + 1
349 }
350 return mask
351}
352
353// ===== ImageS64 -> normalized u8 image ================================
354//
355// Maps the dynamic range of an ImageS64 to [0, 255] for display /
356// further processing. Saturating clamp.
357
358// ===== PYRAMID PRIMITIVES (2026-08-25) =================================
359//
360// Added for coarse-to-fine optical flow. Lucas-Kanade's brightness-constancy derivation assumes
361// motion SMALL relative to its integration window; beyond about a pixel per frame the linearisation
362// stops holding and the estimate collapses toward zero. Halving the image halves the motion, so a
363// pyramid turns a large displacement into a small one at some level. These two primitives are what
364// a pyramid needs and nothing more.
365
366func nx_image_downsample2_into(src: *Image, dst: *Image) -> *Image {
367 // 2x2 BOX average. Box, not nearest: nearest aliases a moving edge in and out between levels
368 // and manufactures motion that is sampling rather than surface.
369 var y: i64 = 0
370 while y < dst.height {
371 var x: i64 = 0
372 while x < dst.width {
373 let a: i64 = nx_image_get(src, x * 2, y * 2, 0)
374 let b: i64 = nx_image_get(src, x * 2 + 1, y * 2, 0)
375 let c: i64 = nx_image_get(src, x * 2, y * 2 + 1, 0)
376 let d: i64 = nx_image_get(src, x * 2 + 1, y * 2 + 1, 0)
377 nx_image_set(dst, x, y, 0, (a + b + c + d) / 4)
378 x = x + 1
379 }
380 y = y + 1
381 }
382 return dst
383}
384
385func nx_image_shift_into(src: *Image, dx: i64, dy: i64, dst: *Image) -> *Image {
386 // Integer translate. Out-of-range reads return 0 via nx_image_get's own bounds guard, which is
387 // the honest behaviour here: a warped-in border has no data and must not borrow a neighbour's.
388 var y: i64 = 0
389 while y < dst.height {
390 var x: i64 = 0
391 while x < dst.width {
392 nx_image_set(dst, x, y, 0, nx_image_get(src, x + dx, y + dy, 0))
393 x = x + 1
394 }
395 y = y + 1
396 }
397 return dst
398}
399
400func nx_image_s64_to_u8(src: *ImageS64) -> *Image {
401 let out: *Image = nx_image_alloc(src.width, src.height, 1)
402 // Find min + max.
403 var min_v: i64 = src.data[0]
404 var max_v: i64 = src.data[0]
405 var i: i64 = 0
406 while i < src.width * src.height {
407 let v: i64 = src.data[i]
408 if v < min_v { min_v = v }
409 if v > max_v { max_v = v }
410 i = i + 1
411 }
412 let range: i64 = max_v - min_v
413 if range <= 0 {
414 return out
415 }
416 var y: i64 = 0
417 while y < src.height {
418 var x: i64 = 0
419 while x < src.width {
420 let v: i64 = nx_image_s64_get(src, x, y)
421 let scaled: i64 = nx_muldiv_i64(v - min_v, 255, range)
422 nx_image_set(out, x, y, 0, scaled)
423 x = x + 1
424 }
425 y = y + 1
426 }
427 return out
428}