QVisionCore: speed up PSNR/SSIM with SSE2 and integral-image box filter (#10)

Fix qSSIM horizontal window sliding: the update
  win_x += col_x[c + SIZE_N] - col_x[c]
computed the column-sum difference as qu32, which wrapped to a large
positive value when the new column had a smaller sum than the leaving
column.  That wrapped value was then zero-extended into the qu64
accumulator, inflating win_x by ~2^32 and corrupting SSIM (sometimes
producing scores greater than 1).  Fix by evaluating the subtraction
in qu64 order: win_x = win_x + col_x[c + SIZE_N] - col_x[c].

Add a regression test: identical bright-left / dark-right images must
score exactly 1.0 regardless of content, which the old code failed.

Co-Authored-By: Claude Sonnet 4.6 <noreply@anthropic.com>

Author: chammoru

QVisionCore/qimage_metrics.c

@@ -1,25 +1,87 @@
 #include "qimage_metrics.h"
 
+#include <stdlib.h>
+
+#if defined(_M_X64) || defined(_M_IX86) || defined(__x86_64__) || defined(__i386__)
+#include <emmintrin.h>
+#define QM_HAS_SSE2 1
+#else
+#define QM_HAS_SSE2 0
+#endif
+
+#if QM_HAS_SSE2
+static qu64 psnr_row_sse2(qu8 *a, qu8 *b, int n)
+{
+	qu64 sum;
+	qu32 lane[4];
+	__m128i acc = _mm_setzero_si128();
+	__m128i zero = _mm_setzero_si128();
+	int j = 0;
+
+	for (; j + 16 <= n; j += 16) {
+		__m128i va = _mm_loadu_si128((const __m128i *)(a + j));
+		__m128i vb = _mm_loadu_si128((const __m128i *)(b + j));
+		__m128i a_lo = _mm_unpacklo_epi8(va, zero);
+		__m128i a_hi = _mm_unpackhi_epi8(va, zero);
+		__m128i b_lo = _mm_unpacklo_epi8(vb, zero);
+		__m128i b_hi = _mm_unpackhi_epi8(vb, zero);
+		__m128i d_lo = _mm_sub_epi16(a_lo, b_lo);
+		__m128i d_hi = _mm_sub_epi16(a_hi, b_hi);
+		// _mm_madd_epi16(x, x) widens int16 to int32 and pairs them.
+		// Per-row max: (W/16) * 8 lanes * (2 * 255^2) ~ 4.16M * (W/16);
+		// for W <= 8192 this is well under 2^31, so 32-bit accum is safe.
+		acc = _mm_add_epi32(acc, _mm_add_epi32(
+			_mm_madd_epi16(d_lo, d_lo),
+			_mm_madd_epi16(d_hi, d_hi)));
+	}
+
+	_mm_storeu_si128((__m128i *)lane, acc);
+	sum = (qu64)lane[0] + lane[1] + lane[2] + lane[3];
+
+	for (; j < n; j++) {
+		long d = a[j] - b[j];
+		sum += (qu64)(d * d);
+	}
+	return sum;
+}
+#endif
+
 double qPSNR(qu8 *src, qu8 *dst, int w, int h, int stride, int px_w)
 {
 	int i, j;
 	qu64 sum_dxd = 0;
 	double MSE, PMSE;
-	const int MAXi = 255; 
-	int margin = (stride - w) * px_w;
+	const int MAXi = 255;
+
+#if QM_HAS_SSE2
+	if (px_w == 1) {
+		for (i = 0; i < h; i++) {
+			sum_dxd += psnr_row_sse2(src, dst, w);
+			src += stride;
+			dst += stride;
+		}
+		MSE = (double)sum_dxd / ((double)w * h);
+		PMSE = MSE / (MAXi * MAXi);
+		return -10 * log10(PMSE);
+	}
+#endif
 
-	for (i = 0; i < h; i++) {
-		for (j = 0; j < w; j++) {
-			long d = *src - *dst;
+	{
+		int margin = (stride - w) * px_w;
 
-			src += px_w;
-			dst += px_w;
+		for (i = 0; i < h; i++) {
+			for (j = 0; j < w; j++) {
+				long d = *src - *dst;
 
-			sum_dxd += (qu64)(d * d);
-		}
+				src += px_w;
+				dst += px_w;
 
-		src += margin;
-		dst += margin;
+				sum_dxd += (qu64)(d * d);
+			}
+
+			src += margin;
+			dst += margin;
+		}
 	}
 
 	MSE = (double)sum_dxd / ((double)w * h);
@@ -36,109 +98,111 @@ double qPSNR(qu8 *src, qu8 *dst, int w, int h, int stride, int px_w)
 
 static const double window_size_reciprocal = 1 / (double)WINDOW_SIZE;
 
-static double ssim_average(qu8 *src, int real_w, int px_w)
+// Integral-image (box-filter) reformulation of the naive sliding-window SSIM.
+// The previous implementation rescanned each 8x8 window five times, making
+// it O(W*H*64); this version keeps running column sums of s, d, s*s, d*d, s*d
+// over the last SIZE_N rows and slides them horizontally to evaluate each
+// window in O(1), yielding O(W*H) overall. Algebraically the variance and
+// covariance use Sx^2/N - mu^2 instead of Sigma(x - mu)^2 / N: mathematically
+// identical, may differ from the old code in the last few ULPs of double
+// rounding.
+double qSSIM(qu8 *src, qu8 *dst, int w, int h, int stride, int px_w)
 {
-	int total = 0;
-	int i, j;
-	double average;
-
-	for (i = 0; i < SIZE_N; i++) {
-		for (j = 0; j < SIZE_N; j++)
-			total += src[j * px_w];
+	const double c1 = K1 * L * K1 * L;
+	const double c2 = K2 * L * K2 * L;
 
-		src += real_w;
-	}
+	int real_w = stride * px_w;
+	int nW = w - SIZE_N + 1;
+	int nH = h - SIZE_N + 1;
+	int r, c, j;
 
-	average = total * window_size_reciprocal;
+	qu32 *col_s, *col_d, *col_ss, *col_dd, *col_sd;
+	double total_ssim = 0.0;
 
-	return average;
-}
+	if (nW <= 0 || nH <= 0)
+		return 0.0;
 
-static double ssim_variance(qu8 *src, int real_w, double average, int px_w)
-{
-	int i, j;
-	double t, sum_square = 0, variance;
+	col_s  = (qu32 *)malloc(sizeof(qu32) * w);
+	col_d  = (qu32 *)malloc(sizeof(qu32) * w);
+	col_ss = (qu32 *)malloc(sizeof(qu32) * w);
+	col_dd = (qu32 *)malloc(sizeof(qu32) * w);
+	col_sd = (qu32 *)malloc(sizeof(qu32) * w);
+	if (!col_s || !col_d || !col_ss || !col_dd || !col_sd) {
+		free(col_s); free(col_d); free(col_ss); free(col_dd); free(col_sd);
+		return 0.0;
+	}
 
-	for (i = 0; i < SIZE_N; i++) {
-		for (j = 0; j < SIZE_N; j++) {
-			t = src[j * px_w] - average;
-			sum_square += t * t;
+	// Seed column sums over rows [0, SIZE_N).
+	for (j = 0; j < w; j++) {
+		col_s[j] = 0; col_d[j] = 0;
+		col_ss[j] = 0; col_dd[j] = 0;
+		col_sd[j] = 0;
+	}
+	for (r = 0; r < SIZE_N; r++) {
+		qu8 *sp = src + r * real_w;
+		qu8 *dp = dst + r * real_w;
+		for (j = 0; j < w; j++) {
+			qu32 s = sp[j * px_w];
+			qu32 d = dp[j * px_w];
+			col_s[j]  += s;
+			col_d[j]  += d;
+			col_ss[j] += s * s;
+			col_dd[j] += d * d;
+			col_sd[j] += s * d;
 		}
-		src += real_w;
 	}
 
-	variance = sum_square * window_size_reciprocal;
+	for (r = 0; r < nH; r++) {
+		qu64 win_s = 0, win_d = 0, win_ss = 0, win_dd = 0, win_sd = 0;
 
-	return variance;
-}
-
-static double ssim_covariance(qu8 *src, qu8 *dst, int real_w,
-							  double s_average, double d_average, int px_w)
-{
-	int i, j;
-	double total = 0, covariance;
-
-	for (i = 0; i < SIZE_N; i++) {
 		for (j = 0; j < SIZE_N; j++) {
-			int idx = j * px_w;
-
-			total += (src[idx] - s_average) * (dst[idx] - d_average);
+			win_s  += col_s[j];
+			win_d  += col_d[j];
+			win_ss += col_ss[j];
+			win_dd += col_dd[j];
+			win_sd += col_sd[j];
 		}
 
-		src += real_w;
-		dst += real_w;
-	}
-
-	covariance = total * window_size_reciprocal;
-
-	return covariance;
-}
-
-static double ssim_window(qu8 *src, qu8 *dst, int real_w, int px_w)
-{
-	const double c1 = K1 * L * K1 * L;
-	const double c2 = K2 * L * K2 * L;
-
-	double u1 = ssim_average(src, real_w, px_w);
-	double u2 = ssim_average(dst, real_w, px_w);
-
-	double a1 = ssim_variance(src, real_w, u1, px_w);
-	double a2 = ssim_variance(dst, real_w, u2, px_w);
-
-	double covariance = ssim_covariance(src, dst, real_w, u1, u2, px_w);
-
-	double ssim = ((2 * u1 * u2 + c1) * (2 * covariance + c2)) /
-		((u1 * u1 + u2 * u2 + c1) * (a1 + a2 + c2));
-
-	return ssim;
-}
-
-double qSSIM(qu8 *src, qu8 *dst, int w, int h, int stride, int px_w)
-{
-	double mean_ssim, total_ssim = 0;
-	int real_w = stride * px_w;
-	int i, j;
-	int nW = w - SIZE_N + 1;
-	int nH = h - SIZE_N + 1;
-
-	if (nW <= 0 || nH <= 0)
-		return 0.0;
-
-	// SIZE_N x SIZE_N linear sliding window
-	// My implementation of SSIM is certainly inefficient,
-	// but it gives values and is easy to understand.
-	for (i = 0; i < nH; i++) {
-		qu8 *src_t, *dst_t;
-		int h_hop = real_w * i;
-
-		src_t = src + h_hop;
-		dst_t = dst + h_hop;
+		for (c = 0; c < nW; c++) {
+			double u1 = (double)win_s  * window_size_reciprocal;
+			double u2 = (double)win_d  * window_size_reciprocal;
+			double a1 = (double)win_ss * window_size_reciprocal - u1 * u1;
+			double a2 = (double)win_dd * window_size_reciprocal - u2 * u2;
+			double cov = (double)win_sd * window_size_reciprocal - u1 * u2;
+
+			double ssim = ((2 * u1 * u2 + c1) * (2 * cov + c2)) /
+				((u1 * u1 + u2 * u2 + c1) * (a1 + a2 + c2));
+			total_ssim += ssim;
+
+			if (c + SIZE_N < w) {
+				win_s  = win_s  + col_s[c + SIZE_N]  - col_s[c];
+				win_d  = win_d  + col_d[c + SIZE_N]  - col_d[c];
+				win_ss = win_ss + col_ss[c + SIZE_N] - col_ss[c];
+				win_dd = win_dd + col_dd[c + SIZE_N] - col_dd[c];
+				win_sd = win_sd + col_sd[c + SIZE_N] - col_sd[c];
+			}
+		}
 
-		for (j = 0; j < nW; j++, src_t += px_w, dst_t += px_w)
-			total_ssim += ssim_window(src_t, dst_t, real_w, px_w);
+		if (r + 1 < nH) {
+			qu8 *sp_old = src + r * real_w;
+			qu8 *dp_old = dst + r * real_w;
+			qu8 *sp_new = src + (r + SIZE_N) * real_w;
+			qu8 *dp_new = dst + (r + SIZE_N) * real_w;
+			for (j = 0; j < w; j++) {
+				qu32 s_old = sp_old[j * px_w];
+				qu32 d_old = dp_old[j * px_w];
+				qu32 s_new = sp_new[j * px_w];
+				qu32 d_new = dp_new[j * px_w];
+				col_s[j]  += s_new       - s_old;
+				col_d[j]  += d_new       - d_old;
+				col_ss[j] += s_new * s_new - s_old * s_old;
+				col_dd[j] += d_new * d_new - d_old * d_old;
+				col_sd[j] += s_new * d_new - s_old * d_old;
+			}
+		}
 	}
 
-	mean_ssim = total_ssim / (nW * nH);
+	free(col_s); free(col_d); free(col_ss); free(col_dd); free(col_sd);
 
-	return mean_ssim;
+	return total_ssim / (nW * nH);
 }

Tests/CoreRegressionTests.c

@@ -317,6 +317,22 @@ static void test_metrics(void)
 	check_near("SSIM constant reference",
 		qSSIM(zero, ten, 8, 8, 8, 1),
 		6.5025 / 106.5025, 0.000000001);
+
+	{
+		/* When the window slides from a bright region into a dark region the
+		 * column sums decrease.  The sliding update
+		 *   win_x += col_x[c + SIZE_N] - col_x[c]
+		 * computed the difference in 32-bit unsigned, which wrapped to a large
+		 * positive value and inflated win_x, corrupting SSIM (sometimes > 1).
+		 * For identical src/dst the result must be exactly 1.0. */
+		qu8 bright_dark[8 * 16];
+		for (y = 0; y < 8; y++) {
+			for (x = 0; x < 16; x++)
+				bright_dark[y * 16 + x] = (x < 8) ? 255 : 0;
+		}
+		check_near("SSIM bright-to-dark identical images",
+			qSSIM(bright_dark, bright_dark, 16, 8, 16, 1), 1.0, 0.000000001);
+	}
 }
 
 int main(void)