nx_self_similarity.nx source
↩ module page · 253 lines · 8762 B
1// nx_self_similarity.nx -- box-counting fractal dimension primitive.
2//
3// Cross-modal aesthetic kernel. Box-counting estimates the fractal
4// dimension D of an edge-image; for natural images D is empirically
5// in [1.0, 2.0] where:
6// D = 1.0 pure line (low fractal complexity)
7// D ~ 1.4 Spehar-2003 peak aesthetic preference
8// D = 2.0 filled plane (maximally complex)
9//
10// Mandelbrot 1967, Falconer "Fractal Geometry" 1990, Spehar et al.
11// 2003 "Universal aesthetic of fractals." Generalizes from images
12// to any 2D signal: spectrograms, attention heatmaps, terrain.
13//
14// Output: Q10 of (D - 1), so 0 = line, 1024 = filled plane. Plus an
15// aesthetic_q10 derived from Spehar-fit Gaussian-shaped peak at
16// D~1.4 (Q10 dim ~410).
17//
18// Why Q10 of (D - 1) instead of Q10 of D directly: 1.0 < D < 2.0 in
19// every realistic case, so the [1, 2] range maps cleanly to [0, Q10]
20// without wasting integer bits encoding the always-present 1.
21//
22// genealogy_id: mandelbrot_1967_fractal + falconer_1990_geometry +
23// spehar_2003_universal_aesthetic_fractals
24// lineage_id: box_counting_fractal_dim_q10
25
26// nx_safety_envelope:
27// intended_use: AUTO_APPLIED -- primitive-specific tuning queued
28// sil_target: SIL1
29// evidence: [bulk_applied_2026-05-16, see-file-comment-for-detail]
30// verdict: NOT_YET_EVALUATED
31
32import "nx_syscalls.nx"
33import "nx_tier.nx"
34import "nx_image.nx"
35
36struct SelfSimilarityReport {
37 fractal_dim_q10: nx_int, // Q10 of (D - 1)
38 aesthetic_q10: nx_int, // Spehar-fit aesthetic preference
39 edge_density_q10: nx_int, // diagnostic; below floor = low fidelity
40 fidelity_q10: nx_int, // [0, 1024]
41}
42
43// ===== Helpers =========================================================
44
45// Integer floor + linear-interp fractional log2 in Q10.
46// log2_q10(2) = 1024
47// log2_q10(4) = 2048
48// log2_q10(8) = 3072
49// log2_q10(16) = 4096
50// log2_q10(3) ~ 1536 (linear-interp; true log2(3) ~ 1623)
51// 4-15% error in the fractional part; sufficient for slope estimation
52// across decades of box counts (where the integer part dominates).
53func _log2_q10(n: nx_int) -> nx_int {
54 if n <= 0 { return 0 }
55 var k: nx_int = 0
56 var x: nx_int = n
57 while x > 1 {
58 x = x / 2
59 k = k + 1
60 }
61 // k = floor(log2(n))
62 var p: nx_int = 1
63 var i: nx_int = 0
64 while i < k {
65 p = p * 2
66 i = i + 1
67 }
68 // n_frac = n - p in [0, p), fraction = n_frac / p in [0, 1)
69 let frac_q10: nx_int = ((n - p) * 1024) / p
70 return k * 1024 + frac_q10
71}
72
73// Sobel-magnitude edge mask: returns 1D buffer of 0/1 per pixel
74// (row-major, w*h entries) above a threshold. Threshold is
75// dynamic-range-based so the function is scale-agnostic.
76func _edge_mask_alloc(gray: *Image) -> *i64 {
77 let w: nx_int = gray.width
78 let h: nx_int = gray.height
79 let n: nx_int = w * h
80 let mask: *i64 = (sys_mmap(n * NX_SIZEOF_NX_INT)) as *i64
81
82 let gx: *ImageS64 = nx_image_sobel_x(gray)
83 let gy: *ImageS64 = nx_image_sobel_y(gray)
84 let mag: *ImageS64 = nx_image_gradient_magnitude(gx, gy)
85
86 // Find dynamic range to set threshold.
87 var lo: nx_int = 0
88 var hi: nx_int = 0
89 var first: nx_int = 1
90 var yy: nx_int = 0
91 while yy < h {
92 var xx: nx_int = 0
93 while xx < w {
94 let v: nx_int = nx_image_s64_get(mag, xx, yy)
95 if first == 1 {
96 lo = v
97 hi = v
98 first = 0
99 } else {
100 if v < lo { lo = v }
101 if v > hi { hi = v }
102 }
103 xx = xx + 1
104 }
105 yy = yy + 1
106 }
107 let range: nx_int = hi - lo
108 let thr: nx_int = lo + (range * 4) / 10 // 40% of range
109
110 var y: nx_int = 0
111 while y < h {
112 var x: nx_int = 0
113 while x < w {
114 let v: nx_int = nx_image_s64_get(mag, x, y)
115 var bit: nx_int = 0
116 if v > thr { bit = 1 }
117 mask[y * w + x] = bit
118 x = x + 1
119 }
120 y = y + 1
121 }
122 return mask
123}
124
125// Count boxes of size s that contain at least one edge pixel.
126func _box_count(mask: *i64, w: nx_int, h: nx_int, s: nx_int) -> nx_int {
127 let nx: nx_int = w / s
128 let ny: nx_int = h / s
129 var occupied: nx_int = 0
130 var by: nx_int = 0
131 while by < ny {
132 var bx: nx_int = 0
133 while bx < nx {
134 // Scan the sxs box at (bx*s, by*s); count if ANY edge pixel.
135 var found: nx_int = 0
136 var dy: nx_int = 0
137 while dy < s {
138 var dx: nx_int = 0
139 while dx < s {
140 let px: nx_int = bx * s + dx
141 let py: nx_int = by * s + dy
142 if mask[py * w + px] == 1 { found = 1 }
143 dx = dx + 1
144 }
145 dy = dy + 1
146 }
147 if found == 1 { occupied = occupied + 1 }
148 bx = bx + 1
149 }
150 by = by + 1
151 }
152 return occupied
153}
154
155// Edge density: fraction of mask cells that are 1. Q10.
156func _edge_density_q10(mask: *i64, w: nx_int, h: nx_int) -> nx_int {
157 let n: nx_int = w * h
158 if n == 0 { return 0 }
159 var sum: nx_int = 0
160 var i: nx_int = 0
161 while i < n {
162 sum = sum + mask[i]
163 i = i + 1
164 }
165 return (sum * 1024) / n
166}
167
168// ===== Public composite ================================================
169//
170// Box-count at four scales (s = 2, 4, 8, 16), regress log(N) on log(s),
171// negative slope = fractal dimension. Q10 throughout.
172//
173// Requires image >= 16x16 (so largest box size fits at least once).
174
175func nx_self_similarity_compute(gray: *Image, report: *SelfSimilarityReport) -> nx_int {
176 let w: nx_int = gray.width
177 let h: nx_int = gray.height
178
179 // Fidelity gate: tiny images can't host a box-counting regression.
180 if w < 16 { return 1 }
181 if h < 16 { return 1 }
182
183 let mask: *i64 = _edge_mask_alloc(gray)
184 let ed_q10: nx_int = _edge_density_q10(mask, w, h)
185 report.edge_density_q10 = ed_q10
186
187 // Below 1% edge density (q10=10) the box counts collapse and the
188 // slope is dominated by quantization noise. Refuse to claim.
189 if ed_q10 < 10 {
190 report.fractal_dim_q10 = 0
191 report.aesthetic_q10 = 0
192 report.fidelity_q10 = 0
193 return 2
194 }
195
196 let n2: nx_int = _box_count(mask, w, h, 2)
197 let n4: nx_int = _box_count(mask, w, h, 4)
198 let n8: nx_int = _box_count(mask, w, h, 8)
199 let n16: nx_int = _box_count(mask, w, h, 16)
200
201 // If any scale collapses to zero, the slope is undefined.
202 if n2 <= 0 { return 3 }
203 if n4 <= 0 { return 3 }
204 if n8 <= 0 { return 3 }
205 if n16 <= 0 { return 3 }
206
207 // log2 in Q10. log2(s) for s=2,4,8,16 -> 1024, 2048, 3072, 4096.
208 let l_n2: nx_int = _log2_q10(n2)
209 let l_n16: nx_int = _log2_q10(n16)
210 // log2(2/16) -> -3 (in log2 units) -> -3*Q10 = -3072 in Q10 units.
211 // D = -slope = (log N(s_min) - log N(s_max)) / (log s_max - log s_min)
212 // = (l_n2 - l_n16) / (4 - 1) in log2 units
213 // = (l_n2_q10 - l_n16_q10) / 3 in Q10 of log2
214 // That's already Q10 of D.
215 let D_q10: nx_int = (l_n2 - l_n16) / 3
216
217 // Clamp D into [1024, 2048] (i.e., D in [1, 2]) -- everything
218 // outside this is either degenerate (D<1) or noise (D>2).
219 var D_clamped: nx_int = D_q10
220 if D_clamped < 1024 { D_clamped = 1024 }
221 if D_clamped > 2048 { D_clamped = 2048 }
222 // Report Q10 of (D - 1) in [0, 1024]
223 let dim_minus_one_q10: nx_int = D_clamped - 1024
224 report.fractal_dim_q10 = dim_minus_one_q10
225
226 // Spehar 2003: aesthetic preference peaks at D ~ 1.4
227 // Q10(D-1) ~ 410 is peak. Triangular falloff: aesthetic =
228 // 1024 - 4 * |dim - 410|, clamped to [0, 1024].
229 var dev: nx_int = dim_minus_one_q10 - 410
230 if dev < 0 { dev = -dev }
231 var aesthetic: nx_int = 1024 - dev * 4
232 if aesthetic < 0 { aesthetic = 0 }
233 if aesthetic > 1024 { aesthetic = 1024 }
234 report.aesthetic_q10 = aesthetic
235
236 // Fidelity: high-edge-density + non-degenerate box counts -> high
237 // confidence. Scale fidelity by edge density (capped at 0.5).
238 var fid: nx_int = ed_q10 * 2
239 if fid > 1024 { fid = 1024 }
240 report.fidelity_q10 = fid
241
242 return 0
243}
244
245// Convenience cross-modal predicate: did this signal hit the
246// Spehar-aesthetic peak band [1.3, 1.5]? Q10 of (D-1) in [307, 512].
247// True (return 1) when fractal_dim_q10 in that band; false (0) otherwise.
248func nx_self_similarity_in_peak_band(report: *SelfSimilarityReport) -> nx_int {
249 let d: nx_int = report.fractal_dim_q10
250 if d < 307 { return 0 }
251 if d > 512 { return 0 }
252 return 1
253}