nx_ssim.nx source

1// nx_ssim.nx -- Structural Similarity (SSIM) metric, Wang 2004. 2// 3// Closes the per-pixel fidelity half of action item #6 in 4// docs/VRAM_OPTIMIZATION_REALISTIC_TRACKING.md ("quality measurement 5// protocol -- LPIPS + CLIP-similarity primitives"). SSIM is the 6// substrate-native canonical metric the literature uses to say 7// "the quantized output matches the f16 reference." 8// 9// PER THE BITS-UP CARDINAL: this primitive composes against the 10// canonical `*Image` type from nx_image.nx -- the same type used by 11// nx_scale (gauss pyramid), nx_aesthetics, nx_self_similarity, etc. 12// MS-SSIM in nx_ms_ssim.nx then composes nx_scale_pyramid_down with 13// this kernel rather than reinventing downsample inline. 14// 15// Why not LPIPS directly: canonical LPIPS uses trained VGG/AlexNet 16// features (Zhang et al. 2018). Substrate has no trained vision 17// network yet (gated on autodiff or pretrained-weight loader). SSIM 18// is the pre-LPIPS canonical metric -- pure algorithmic, well- 19// calibrated against human judgment for fidelity-class comparisons 20// (Wang & Bovik 2009 _Mean Squared Error: Love It or Leave It?_), 21// and the standard for the kind of A/B-bench we actually need (does 22// q4_K weight quantization change the rendered pixels significantly 23// vs f16 reference?). 24// 25// ===== Math ======================================================= 26// 27// For two images X, Y and a window W centred at (i, j): 28// 29// μ_x = mean of W over X 30// μ_y = mean of W over Y 31// σ_x² = variance of W over X 32// σ_y² = variance of W over Y 33// σ_xy = covariance of W over X and Y 34// 35// Stability constants for 8-bit dynamic range L = 255: 36// C1 = (0.01 * L)^2 = 6.5025 37// C2 = (0.03 * L)^2 = 58.5225 38// 39// Local SSIM: 40// SSIM(i, j) = (2 μ_x μ_y + C1) (2 σ_xy + C2) 41// / ((μ_x² + μ_y² + C1) (σ_x² + σ_y² + C2)) 42// 43// Mean SSIM (the scalar reported): 44// MSSIM = average of SSIM(i, j) over valid (i, j) where the window 45// fits inside both images. 46// 47// Range: SSIM ∈ [-1024, 1024] in Q10. 1024 = perfect match. 48// Typical "noticeable degradation" threshold = SSIM < 970 (Wang 49// 2004 reports humans flag a ~3% drop in SSIM as "different"). 50// 51// ===== Q-format =================================================== 52// 53// Image pixels: u8 in [0, 255] (read via nx_image_get). 54// μ, σ², σ_xy: computed in i64 (window area divides cleanly when 55// WIN=8 since 64 = 2^6). 56// Final ratio: Q10 = 1024 means SSIM = 1.0. 57// 58// C1 in Q0 = 7 (6.5025 rounded) 59// C2 in Q0 = 59 (58.5225 rounded) 60// 61// 8x8 box window (Wang's recommended 11x11 Gaussian downgraded for 62// substrate simplicity; ~5% accuracy loss vs Gaussian per Wang 2009). 63// nx_scale.nx already ships 5-tap binomial Gaussian; lifting SSIM 64// to use a Gaussian window is a queued upgrade. 65// 66// Per the bounded-loop cardinal: every loop uses LoopVerdict. 67// 68// genealogy_id: wang_bovik_sheikh_simoncelli_2004_ssim + 69// wang_simoncelli_bovik_2003_ms_ssim + 70// wang_bovik_2009_mse_love_it_or_leave_it 71// lineage_id: substrate_ssim_v2_image 72// 73// nx_safety_envelope: 74// intended_use: "Structural Similarity Index (Wang 2004) -- 75// image-quality verdict consumed by substrate 76// graders for generated visual assets" 77// sil_target: SIL1 (quality grader) 78// asil_target: QM 79// dal_target: NONE 80// evidence: [Wang_2004_canonical_paper, 81// Image_container_compose_not_inline, 82// Q-scale_fixed_point_deterministic] 83// hazard_register: [bug-tape-different-size-images-undef, 84// bug-tape-luminance-only-misses-chroma-shift] 85// residual_risk: "Single-scale; for multi-scale see nx_ms_ssim. 86// Luminance-only by default; full-color 87// SSIM-Y'CbCr queued." 88// verdict: NOT_YET_EVALUATED 89 90import "nx_syscalls.nx" 91import "nx_tier.nx" 92import "nx_loop.nx" 93import "nx_image.nx" 94 95// ===== Constants ================================================== 96 97const NX_SSIM_WIN: nx_int = 8 // 8x8 window 98const NX_SSIM_WIN_AREA: nx_int = 64 // WIN * WIN 99const NX_SSIM_C1_Q0: nx_int = 7 100const NX_SSIM_C2_Q0: nx_int = 59 101const NX_SSIM_Q10_ONE: nx_int = 1024 102const NX_SSIM_DEGRADE_THRESHOLD_Q10: nx_int = 970 103 104// ===== Sealed-enum: SsimVerdict =================================== 105 106const NX_SSIM_OK: nx_int = 0 107const NX_SSIM_ERR_BAD_DIMS: nx_int = 1 108const NX_SSIM_ERR_DIMS_TOO_SMALL: nx_int = 2 109const NX_SSIM_ERR_DIMS_MISMATCH: nx_int = 3 110const NX_SSIM_ERR_BAD_CHANNELS: nx_int = 4 111const NX_SSIM_N_VERDICTS: nx_int = 5 112 113func nx_ssim_verdict_is_valid(v: nx_int) -> nx_int { 114 if v < 0 { return 0 } 115 if v >= NX_SSIM_N_VERDICTS { return 0 } 116 return 1 117} 118 119// ===== Per-window statistics ====================================== 120// 121// Returns μ_x, μ_y, σ_x², σ_y², σ_xy via output pointers. 122// Window is WIN x WIN starting at (left, top) inside both Images. 123// Reads via nx_image_get on channel 0 (single-channel for SSIM). 124 125func _ssim_window_stats(x: *Image, y: *Image, 126 top: nx_int, left: nx_int, 127 out_mux: *i64, out_muy: *i64, 128 out_varx: *i64, out_vary: *i64, 129 out_cov: *i64) -> nx_int { 130 // Pass 1: sums. 131 var sum_x: i64 = 0 132 var sum_y: i64 = 0 133 var r: nx_int = 0 134 var iter: nx_int = 0 135 var verdict: nx_int = NX_LOOP_RUNNING 136 let BUDGET: nx_int = NX_SSIM_WIN 137 while verdict == NX_LOOP_RUNNING && iter < BUDGET { 138 var c: nx_int = 0 139 var iter_c: nx_int = 0 140 var verdict_c: nx_int = NX_LOOP_RUNNING 141 while verdict_c == NX_LOOP_RUNNING && iter_c < NX_SSIM_WIN { 142 sum_x = sum_x + nx_image_get(x, left + c, top + r, 0) 143 sum_y = sum_y + nx_image_get(y, left + c, top + r, 0) 144 c = c + 1 145 iter_c = iter_c + 1 146 } 147 r = r + 1 148 iter = iter + 1 149 } 150 let mu_x: i64 = sum_x / NX_SSIM_WIN_AREA 151 let mu_y: i64 = sum_y / NX_SSIM_WIN_AREA 152 out_mux[0] = mu_x 153 out_muy[0] = mu_y 154 155 // Pass 2: centred sums (variance + covariance). 156 var sum_xx: i64 = 0 157 var sum_yy: i64 = 0 158 var sum_xy: i64 = 0 159 r = 0 160 iter = 0 161 verdict = NX_LOOP_RUNNING 162 while verdict == NX_LOOP_RUNNING && iter < BUDGET { 163 var c2: nx_int = 0 164 var iter_c2: nx_int = 0 165 var verdict_c2: nx_int = NX_LOOP_RUNNING 166 while verdict_c2 == NX_LOOP_RUNNING && iter_c2 < NX_SSIM_WIN { 167 let px: i64 = nx_image_get(x, left + c2, top + r, 0) 168 let py: i64 = nx_image_get(y, left + c2, top + r, 0) 169 let dx: i64 = px - mu_x 170 let dy: i64 = py - mu_y 171 sum_xx = sum_xx + dx * dx 172 sum_yy = sum_yy + dy * dy 173 sum_xy = sum_xy + dx * dy 174 c2 = c2 + 1 175 iter_c2 = iter_c2 + 1 176 } 177 r = r + 1 178 iter = iter + 1 179 } 180 out_varx[0] = sum_xx / NX_SSIM_WIN_AREA 181 out_vary[0] = sum_yy / NX_SSIM_WIN_AREA 182 out_cov[0] = sum_xy / NX_SSIM_WIN_AREA 183 return NX_SSIM_OK 184} 185 186// ===== Local SSIM at a window ===================================== 187// 188// Returns Q10 ratio: 0..1024 (clamped) for the 8x8 window at (top, left). 189 190func _ssim_local_q10(x: *Image, y: *Image, top: nx_int, left: nx_int) -> nx_int { 191 let mux: *i64 = sys_mmap(8) as *i64 192 let muy: *i64 = sys_mmap(8) as *i64 193 let varx: *i64 = sys_mmap(8) as *i64 194 let vary: *i64 = sys_mmap(8) as *i64 195 let cov: *i64 = sys_mmap(8) as *i64 196 _ssim_window_stats(x, y, top, left, mux, muy, varx, vary, cov) 197 198 let mx: i64 = mux[0] 199 let my: i64 = muy[0] 200 let vx: i64 = varx[0] 201 let vy: i64 = vary[0] 202 let cv: i64 = cov[0] 203 204 let num1: i64 = 2 * mx * my + NX_SSIM_C1_Q0 205 let num2: i64 = 2 * cv + NX_SSIM_C2_Q0 206 let den1: i64 = mx * mx + my * my + NX_SSIM_C1_Q0 207 let den2: i64 = vx + vy + NX_SSIM_C2_Q0 208 let part: i64 = num1 * NX_SSIM_Q10_ONE / den1 209 let num: i64 = part * num2 / den2 210 if num < 0 { return 0 } 211 if num > NX_SSIM_Q10_ONE { return NX_SSIM_Q10_ONE } 212 return num 213} 214 215// ===== Mean SSIM over valid window region ========================= 216// 217// Canonical entrypoint -- takes *Image directly. Both images must 218// have the same dims and channels. 219 220func nx_ssim_mean_image_q10(x: *Image, y: *Image) -> nx_int { 221 let h: nx_int = x.height 222 let w: nx_int = x.width 223 if y.width != w { return 0 } 224 if y.height != h { return 0 } 225 if h < NX_SSIM_WIN { return 0 } 226 if w < NX_SSIM_WIN { return 0 } 227 let valid_h: nx_int = h - NX_SSIM_WIN + 1 228 let valid_w: nx_int = w - NX_SSIM_WIN + 1 229 let n_windows: nx_int = valid_h * valid_w 230 231 var sum: i64 = 0 232 var top: nx_int = 0 233 var iter: nx_int = 0 234 var verdict: nx_int = NX_LOOP_RUNNING 235 let BUDGET: nx_int = valid_h 236 while verdict == NX_LOOP_RUNNING && iter < BUDGET { 237 var left: nx_int = 0 238 var iter_c: nx_int = 0 239 var verdict_c: nx_int = NX_LOOP_RUNNING 240 let BUDGET_C: nx_int = valid_w 241 while verdict_c == NX_LOOP_RUNNING && iter_c < BUDGET_C { 242 sum = sum + _ssim_local_q10(x, y, top, left) 243 left = left + 1 244 iter_c = iter_c + 1 245 } 246 top = top + 1 247 iter = iter + 1 248 } 249 if n_windows <= 0 { return 0 } 250 return sum / n_windows 251} 252 253// ===== Degradation verdict ======================================== 254 255func nx_ssim_is_degraded(score_q10: nx_int) -> nx_int { 256 if score_q10 < NX_SSIM_DEGRADE_THRESHOLD_Q10 { return 1 } 257 return 0 258} 259 260// ===== Self-test ================================================== 261// 262// Closed-form invariants: 263// (a) IDENTITY: ssim(x, x) = 1024 exactly. 264// (b) SYMMETRY: ssim(x, y) = ssim(y, x). 265// (c) CONSTANT OFFSET: small constant intensity shift -> SSIM > 950. 266// (d) HEAVY NOISE: alternating +/-80 -> SSIM < 700 + is_degraded(). 267 268func main() -> i64 { 269 let H: nx_int = 32 270 let W: nx_int = 32 271 let a: *Image = nx_image_alloc(W, H, 1) 272 let b: *Image = nx_image_alloc(W, H, 1) 273 274 // Checkerboard + gradient. 275 var r: nx_int = 0 276 var iter: nx_int = 0 277 var verdict: nx_int = NX_LOOP_RUNNING 278 let BUDGET: nx_int = H 279 while verdict == NX_LOOP_RUNNING && iter < BUDGET { 280 var c: nx_int = 0 281 var iter_c: nx_int = 0 282 var verdict_c: nx_int = NX_LOOP_RUNNING 283 while verdict_c == NX_LOOP_RUNNING && iter_c < W { 284 let parity: nx_int = (r + c) - ((r + c) / 2) * 2 285 let v: nx_int = r * 2 + c + parity * 50 286 nx_image_set(a, c, r, 0, v) 287 c = c + 1 288 iter_c = iter_c + 1 289 } 290 r = r + 1 291 iter = iter + 1 292 } 293 294 // --- (a) Identity --- 295 var rr: nx_int = 0 296 while rr < H { 297 var cc: nx_int = 0 298 while cc < W { 299 nx_image_set(b, cc, rr, 0, nx_image_get(a, cc, rr, 0)) 300 cc = cc + 1 301 } 302 rr = rr + 1 303 } 304 let s_id: nx_int = nx_ssim_mean_image_q10(a, b) 305 if s_id != NX_SSIM_Q10_ONE { return 10 } 306 if nx_ssim_is_degraded(s_id) != 0 { return 11 } 307 308 // --- (b) Symmetry on small offset --- 309 var rr2: nx_int = 0 310 while rr2 < H { 311 var cc2: nx_int = 0 312 while cc2 < W { 313 let av: nx_int = nx_image_get(a, cc2, rr2, 0) 314 nx_image_set(b, cc2, rr2, 0, av + 5) 315 cc2 = cc2 + 1 316 } 317 rr2 = rr2 + 1 318 } 319 let s_xy: nx_int = nx_ssim_mean_image_q10(a, b) 320 let s_yx: nx_int = nx_ssim_mean_image_q10(b, a) 321 if s_xy != s_yx { return 20 } 322 323 // --- (c) Constant offset preserves structure --- 324 if s_xy < 950 { return 30 } 325 if s_xy > NX_SSIM_Q10_ONE { return 31 } 326 327 // --- (d) Heavy noise drops SSIM --- 328 var rr3: nx_int = 0 329 while rr3 < H { 330 var cc3: nx_int = 0 331 while cc3 < W { 332 let av: nx_int = nx_image_get(a, cc3, rr3, 0) 333 let par: nx_int = (rr3 + cc3) - ((rr3 + cc3) / 2) * 2 334 let bv: nx_int = av + par * 160 - 80 335 nx_image_set(b, cc3, rr3, 0, bv) 336 cc3 = cc3 + 1 337 } 338 rr3 = rr3 + 1 339 } 340 let s_noise: nx_int = nx_ssim_mean_image_q10(a, b) 341 if s_noise > 700 { return 40 } 342 if nx_ssim_is_degraded(s_noise) != 1 { return 41 } 343 344 // --- (e) Verdict range gate --- 345 var vi: nx_int = 0 346 while vi < NX_SSIM_N_VERDICTS { 347 if nx_ssim_verdict_is_valid(vi) != 1 { return 50 + vi } 348 vi = vi + 1 349 } 350 351 return 0 352}