nx_quality_metric.nx source
↩ module page · 104 lines · 4632 B
1// nx_quality_metric.nx -- Computes video quality metrics like MSE, RMSE, PSNR, and SSIM for rate-distortion analysis in codecs.
2const QM_MAGIC_65536: i64 = 65536
3const QM_MAGIC_131072: i64 = 131072
4const QM_MAGIC_32768: i64 = 32768
5const QM_MAGIC_9900: i64 = 9900
6const QM_MAGIC_65025: i64 = 65025
7const QM_MAGIC_217707: i64 = 217707
8// nx_quality_metric.nx -- sovereign objective video-quality metrics (PSNR / MSE / RMSE) so every codec -- ours
9// AND the competitors' (x264/VP9/AV1 reconstructions) -- lands on ONE measured rate-distortion axis. This is the
10// core of the live benchmark harness: distortion(MSE) vs rate(bytes) gives the RD curve that proves exceed/parity.
11// Integer-only (no FPU): MSE/RMSE are exact; PSNR uses a fixed-point log2 (good to ~0.1 dB). license_tier: ORIGINAL
12
13// mean squared error between two 8bpp planes of n pixels (exact)
14func qm_mse(a: *u8, b: *u8, n: i64) -> i64 {
15 if n <= 0 { return 0 }
16 var s: i64 = 0
17 var i: i64 = 0
18 while i < n { let d: i64 = (a[i] as i64) - (b[i] as i64); s = s + d*d; i = i + 1 }
19 return s / n
20}
21// integer sqrt (Newton)
22func qm_isqrt(v: i64) -> i64 {
23 if v <= 0 { return 0 }
24 var x: i64 = v
25 var y: i64 = (x + 1) / 2
26 var go: i64 = 1
27 while go == 1 { if y < x { x = y; y = (x + v/x) / 2 } else { go = 0 } }
28 return x
29}
30// RMSE = sqrt(MSE): average pixel-error magnitude (human-readable distortion, ranks like PSNR)
31func qm_rmse(a: *u8, b: *u8, n: i64) -> i64 { return qm_isqrt(qm_mse(a, b, n)) }
32
33// log2(x) in Q16 fixed point (x >= 1): integer part by MSB, fractional part by iterative squaring of the mantissa.
34func qm_log2_q16(x: i64) -> i64 {
35 if x <= 1 { return 0 }
36 var n: i64 = 0
37 var t: i64 = x
38 while t > 1 { t = t / 2; n = n + 1 } // n = floor(log2(x))
39 var m: i64 = x << 16
40 var i: i64 = 0
41 while i < n { m = m / 2; i = i + 1 } // mantissa in Q16, range [QM_MAGIC_65536, QM_MAGIC_131072) = [1,2)
42 var frac: i64 = 0
43 var bit: i64 = QM_MAGIC_32768 // 0.5 in Q16
44 var k: i64 = 0
45 while k < 16 {
46 m = (m * m) >> 16
47 if m >= QM_MAGIC_131072 { m = m / 2; frac = frac + bit }
48 bit = bit / 2
49 k = k + 1
50 }
51 return (n << 16) + frac
52}
53// PSNR in centi-dB (x100). PSNR = 10*log10(MAX^2/MSE), MAX=255 -> MAX^2=65025. Identical -> capped 99.00 dB.
54func qm_psnr_cdb(a: *u8, b: *u8, n: i64) -> i64 {
55 let mse: i64 = qm_mse(a, b, n)
56 if mse <= 0 { return QM_MAGIC_9900 }
57 let diff: i64 = qm_log2_q16(QM_MAGIC_65025) - qm_log2_q16(mse) // log2(QM_MAGIC_65025/MSE) in Q16
58 if diff <= 0 { return 0 }
59 return diff * 1000 / QM_MAGIC_217707 // /(QM_MAGIC_65536*log2(10)) -> 1000*log10(ratio) = PSNR*100
60}
61
62// SSIM x1000 (1000 = identical), averaged over 8x8 windows. Unlike PSNR/MSE, SSIM measures STRUCTURE (luminance,
63// contrast, correlation), so a uniform brightness shift barely dents it while real structural damage tanks it --
64// it tracks human perception far better, and is what modern codec RD comparisons report. C1/C2 are the standard
65// (0.01*255)^2 / (0.03*255)^2 constants. Integer fixed-point via the windowed-sums identity.
66const QM_SSIM_C1: i64 = 6
67const QM_SSIM_C2: i64 = 59
68func qm_ssim_milli(a: *u8, b: *u8, W: i64, H: i64) -> i64 {
69 let WIN: i64 = 8
70 let N: i64 = WIN * WIN
71 var total: i64 = 0
72 var nwin: i64 = 0
73 var wy: i64 = 0
74 while wy + WIN <= H {
75 var wx: i64 = 0
76 while wx + WIN <= W {
77 var Sx: i64 = 0; var Sy: i64 = 0; var Sxx: i64 = 0; var Syy: i64 = 0; var Sxy: i64 = 0
78 var yy: i64 = 0
79 while yy < WIN {
80 var xx: i64 = 0
81 while xx < WIN {
82 let px: i64 = a[(wy+yy)*W + (wx+xx)] as i64
83 let py: i64 = b[(wy+yy)*W + (wx+xx)] as i64
84 Sx = Sx + px; Sy = Sy + py
85 Sxx = Sxx + px*px; Syy = Syy + py*py; Sxy = Sxy + px*py
86 xx = xx + 1
87 }
88 yy = yy + 1
89 }
90 let A1: i64 = 2*Sx*Sy + QM_SSIM_C1*N*N
91 let B1: i64 = Sx*Sx + Sy*Sy + QM_SSIM_C1*N*N
92 let A2: i64 = 2*(N*Sxy - Sx*Sy) + QM_SSIM_C2*N*N
93 let B2: i64 = (N*Sxx - Sx*Sx) + (N*Syy - Sy*Sy) + QM_SSIM_C2*N*N
94 var s: i64 = 0
95 if B1 != 0 { if B2 != 0 { s = A1*1000/B1 * A2 / B2 } } // (A1/B1)*(A2/B2)*1000, overflow-safe order
96 total = total + s
97 nwin = nwin + 1
98 wx = wx + WIN
99 }
100 wy = wy + WIN
101 }
102 if nwin == 0 { return 0 }
103 return total / nwin
104}