Quality pass: faster DCT, consistent numbers, standard SSIM, accessibility, two inspector features

Makefile

@@ -1,11 +1,11 @@
 # --- CONFIGURATION ---
 WEB_COMPILER = emcc
 # C functions exported to JavaScript (keep in sync with web/src/types/wasm.d.ts)
-WEB_EXPORTS = _init_session _process_image _rd_encode _get_rd_bit_estimate \
+WEB_EXPORTS = _init_session _process_image _preview_image _rd_encode _get_rd_bit_estimate \
               _get_view_ptr _set_view_tint _set_artifact_gain \
               _get_psnr_y _get_psnr_cr _get_psnr_cb \
               _get_ssim_y _get_ssim_cr _get_ssim_cb \
-              _get_last_bit_estimate _inspect_block_data _get_coeff_histogram \
+              _get_last_bit_estimate _inspect_block_data _get_coeff_stats \
               _init_me_session _run_motion_estimation _get_mv_count \
               _get_me_residual_ptr _get_search_steps _get_search_step_count \
               _malloc _free

core/inc/ImageCodec.h

@@ -65,6 +65,22 @@ class ImageCodec {
 
     bool isWavelet() const { return m_transformType != TransformType::DCT; }
 
+    /*
+    * Distribution of a channel's AC / detail coefficients (DC and LL excluded)
+    * under this codec's transform, before quantization. The channel is treated
+    * as luma (DCT uses the luma quantization table).
+    *   histogram[0]   — fraction with |c| < 1
+    *   histogram[b]   — fraction with 2^(b-1) <= |c| < 2^b; the last bin is open-ended
+    *   zeroFraction   — fraction that quantize to zero at this codec's quality
+    */
+    struct CoefficientStats {
+        static constexpr int NUM_BINS = 12;
+        double histogram[NUM_BINS] = {};
+        double zeroFraction = 0.0;
+    };
+
+    CoefficientStats coefficientStats(const Image& channel) const;
+
     /*
     * Downsamples a chroma channel by box-averaging (2×1 for 4:2:2, 2×2 for 4:2:0).
     */

core/src/CodecAnalysis.cpp

@@ -19,6 +19,7 @@
 #include <cmath>
 #include "colorspace.h"
 #include <algorithm>
+#include <vector>
 
 Image CodecAnalysis::computeArtifactMap(
     const Image& original,
@@ -136,111 +137,109 @@ double CodecAnalysis::computePSNR(const Image& I1, const Image& I2) {
     return 10.0 * log10((255.0 * 255.0) / mse);
 }
 
-// Helper for SSIM
-static double getMean(const double* data, int width, int height, int x, int y, int kernelSize) {
+// --- SSIM (Wang, Bovik, Sheikh & Simoncelli, 2004) ---------------------------
+// 11×11 Gaussian window with σ = 1.5, constants K1 = 0.01, K2 = 0.03 for 8-bit
+// data, evaluated at every pixel and averaged over the region where the window
+// fits inside the image (as in the reference ssim_index.m and scikit-image
+// with gaussian_weights=True). Images smaller than the window are handled by
+// clamping at the borders and averaging over every pixel.
+
+static const int SSIM_RADIUS = 5;
+static const double SSIM_SIGMA = 1.5;
+
+static std::vector<double> ssimKernel() {
+    std::vector<double> k(2 * SSIM_RADIUS + 1);
     double sum = 0.0;
-    int count = 0;
-    int half = kernelSize / 2;
-
-    for (int j = -half; j <= half; ++j) {
-        for (int i = -half; i <= half; ++i) {
-            int cx = x + i;
-            int cy = y + j;
-            if (cx >= 0 && cx < width && cy >= 0 && cy < height) {
-                sum += data[cy * width + cx];
-                count++;
-            }
-        }
+    for (int i = -SSIM_RADIUS; i <= SSIM_RADIUS; ++i) {
+        k[i + SSIM_RADIUS] = std::exp(-(i * i) / (2.0 * SSIM_SIGMA * SSIM_SIGMA));
+        sum += k[i + SSIM_RADIUS];
     }
-    return (count > 0) ? (sum / count) : 0.0;
+    for (double& v : k) v /= sum;
+    return k;
 }
 
-static double getVariance(const double* data, int width, int height, int x, int y, int kernelSize, double mean) {
-    double sum = 0.0;
-    int count = 0;
-    int half = kernelSize / 2;
-
-    for (int j = -half; j <= half; ++j) {
-        for (int i = -half; i <= half; ++i) {
-            int cx = x + i;
-            int cy = y + j;
-            if (cx >= 0 && cx < width && cy >= 0 && cy < height) {
-                double diff = data[cy * width + cx] - mean;
-                sum += diff * diff;
-                count++;
-            }
+// One 1D pass of the symmetric (2r+1)-tap Gaussian over n samples spaced
+// `stride` apart, with edge clamping. The kernel is symmetric, so taps are
+// paired (half the multiplies), and clamping is only done near the ends.
+static void gaussian1D(const double* src, double* dst, int n, size_t stride, const std::vector<double>& k) {
+    const int r = SSIM_RADIUS;
+    auto at = [&](int i) { return src[(size_t)std::max(0, std::min(n - 1, i)) * stride]; };
+    for (int x = 0; x < n; ++x) {
+        double acc = k[r] * src[(size_t)x * stride];
+        if (x >= r && x < n - r) {
+            for (int i = 1; i <= r; ++i)
+                acc += k[r + i] * (src[(size_t)(x - i) * stride] + src[(size_t)(x + i) * stride]);
+        } else {
+            for (int i = 1; i <= r; ++i)
+                acc += k[r + i] * (at(x - i) + at(x + i));
         }
+        dst[(size_t)x * stride] = acc;
     }
-    return (count > 0) ? (sum / count) : 0.0;
 }
 
-static double getCovariance(const double* d1, const double* d2, int width, int height, int x, int y, int kernelSize, double mean1, double mean2) {
+// Separable Gaussian blur of a w×h plane with edge clamping.
+static std::vector<double> gaussianBlur(const std::vector<double>& src, int w, int h,
+                                        const std::vector<double>& k) {
+    std::vector<double> tmp(src.size()), out(src.size());
+    for (int y = 0; y < h; ++y)
+        gaussian1D(src.data() + (size_t)y * w, tmp.data() + (size_t)y * w, w, 1, k);
+    for (int x = 0; x < w; ++x)
+        gaussian1D(tmp.data() + x, out.data() + x, h, (size_t)w, k);
+    return out;
+}
+
+static double ssimPlane(const std::vector<double>& a, const std::vector<double>& b, int w, int h) {
+    const double C1 = (0.01 * 255.0) * (0.01 * 255.0);
+    const double C2 = (0.03 * 255.0) * (0.03 * 255.0);
+    const std::vector<double> k = ssimKernel();
+
+    std::vector<double> aa(a.size()), bb(a.size()), ab(a.size());
+    for (size_t i = 0; i < a.size(); ++i) {
+        aa[i] = a[i] * a[i];
+        bb[i] = b[i] * b[i];
+        ab[i] = a[i] * b[i];
+    }
+    const auto muA = gaussianBlur(a, w, h, k);
+    const auto muB = gaussianBlur(b, w, h, k);
+    const auto eAA = gaussianBlur(aa, w, h, k);
+    const auto eBB = gaussianBlur(bb, w, h, k);
+    const auto eAB = gaussianBlur(ab, w, h, k);
+
+    const bool fits = w > 2 * SSIM_RADIUS && h > 2 * SSIM_RADIUS;
+    const int r = fits ? SSIM_RADIUS : 0;
     double sum = 0.0;
-    int count = 0;
-    int half = kernelSize / 2;
-
-    for (int j = -half; j <= half; ++j) {
-        for (int i = -half; i <= half; ++i) {
-            int cx = x + i;
-            int cy = y + j;
-            if (cx >= 0 && cx < width && cy >= 0 && cy < height) {
-                double diff1 = d1[cy * width + cx] - mean1;
-                double diff2 = d2[cy * width + cx] - mean2;
-                sum += diff1 * diff2;
-                count++;
-            }
+    long long count = 0;
+    for (int y = r; y < h - r; ++y)
+        for (int x = r; x < w - r; ++x) {
+            const size_t i = (size_t)y * w + x;
+            const double varA = eAA[i] - muA[i] * muA[i];
+            const double varB = eBB[i] - muB[i] * muB[i];
+            const double cov  = eAB[i] - muA[i] * muB[i];
+            sum += ((2.0 * muA[i] * muB[i] + C1) * (2.0 * cov + C2)) /
+                   ((muA[i] * muA[i] + muB[i] * muB[i] + C1) * (varA + varB + C2));
+            ++count;
         }
-    }
-    return (count > 0) ? (sum / count) : 0.0;
+    return count > 0 ? sum / static_cast<double>(count) : 1.0;
 }
 
 double CodecAnalysis::computeSSIM(const Image& I1, const Image& I2) {
     if (I1.width() != I2.width() || I1.height() != I2.height() || I1.channels() != I2.channels()) {
         return 0.0;
     }
 
-    const double C1 = 6.5025;  // (0.01 * 255)^2
-    const double C2 = 58.5225; // (0.03 * 255)^2
-
-    int width = I1.width();
-    int height = I1.height();
-    const double* d1 = I1.data();
-    const double* d2 = I2.data();
-
-    // Using a simpler block-based approach for performance instead of full sliding window Gaussian
-    // But to be somewhat accurate we'll use a sliding window with stride
-    // For this implementation, let's do a basic sliding window with a small stride for speed
-    // or just checking every Nth pixel.
-
-    // Let's do a proper sliding window but maybe with skip if needed. 
-    // For now, full sliding window.
-
-    // NOTE: This is a simplified SSIM. A full Gaussian weighting is standard but slower.
-    // Using an 8x8 uniform window is a common approximation.
-    int kernelSize = 8;
-    int stride = 4; // Speed up calculation
-
-    double mssim = 0.0;
-    int blocks = 0;
-
-    for (int y = 0; y < height; y += stride) {
-        for (int x = 0; x < width; x += stride) {
-            double ux = getMean(d1, width, height, x, y, kernelSize);
-            double uy = getMean(d2, width, height, x, y, kernelSize);
-
-            double sigx2 = getVariance(d1, width, height, x, y, kernelSize, ux);
-            double sigy2 = getVariance(d2, width, height, x, y, kernelSize, uy);
-            double sigxy = getCovariance(d1, d2, width, height, x, y, kernelSize, ux, uy);
-
-            double num = (2 * ux * uy + C1) * (2 * sigxy + C2);
-            double den = (ux * ux + uy * uy + C1) * (sigx2 + sigy2 + C2);
-
-            mssim += (num / den);
-            blocks++;
+    // Multi-channel input: mean SSIM over channels.
+    const int w = I1.width(), h = I1.height(), ch = I1.channels();
+    const size_t n = static_cast<size_t>(w) * h;
+    double total = 0.0;
+    std::vector<double> a(n), b(n);
+    for (int c = 0; c < ch; ++c) {
+        for (size_t i = 0; i < n; ++i) {
+            a[i] = I1.data()[i * ch + c];
+            b[i] = I2.data()[i * ch + c];
         }
+        total += ssimPlane(a, b, w, h);
     }
-
-    return (blocks > 0) ? (mssim / blocks) : 0.0;
+    return total / ch;
 }
 
 CodecMetrics CodecAnalysis::computeMetrics(const Image& originalBgr, const Image& reconstructedBgr) {

core/src/ImageCodec.cpp

@@ -237,6 +237,62 @@ Image ImageCodec::processChannelDWT(const Image& channel, bool isChroma)
     return result;
 }
 
+// Log2 magnitude bin: 0 for |c| < 1, b for 2^(b-1) <= |c| < 2^b, capped.
+static int magnitudeBin(double c) {
+    const double mag = std::abs(c);
+    if (mag < 1.0) return 0;
+    const int bin = 1 + static_cast<int>(std::floor(std::log2(mag)));
+    return std::min(bin, ImageCodec::CoefficientStats::NUM_BINS - 1);
+}
+
+ImageCodec::CoefficientStats ImageCodec::coefficientStats(const Image& channel) const
+{
+    CoefficientStats stats;
+    long long total = 0;
+    long long zeroed = 0;
+
+    auto record = [&](double coeff, double step) {
+        stats.histogram[magnitudeBin(coeff)] += 1.0;
+        const bool zero = m_enableQuantization ? std::round(coeff / step) == 0.0 : coeff == 0.0;
+        if (zero) ++zeroed;
+        ++total;
+    };
+
+    if (isWavelet()) {
+        std::vector<double> buf;
+        int W = 0, H = 0, levels = 0;
+        transformChannelDWT(channel, buf, W, H, levels);
+        const double baseStep = DWT_BASE_STEP * qualityScale(m_quality);
+        for (int y = 0; y < H; ++y)
+            for (int x = 0; x < W; ++x) {
+                if (dwtSubbandIndex(x, y, W, H, levels) == 0) continue; // LL
+                record(buf[(size_t)y * W + x], dwtQuantStep(x, y, W, H, levels, baseStep));
+            }
+    } else {
+        // Same block extraction (edge replication, level shift) as processChannel.
+        const int width = channel.width();
+        const int height = channel.height();
+        const double* data = channel.data();
+        for (int y = 0; y < height; y += 8)
+            for (int x = 0; x < width; x += 8) {
+                double block[8][8], coeffs[8][8];
+                for (int i = 0; i < 8; ++i)
+                    for (int j = 0; j < 8; ++j)
+                        block[i][j] = data[std::min(y + i, height - 1) * width + std::min(x + j, width - 1)] - 128.0;
+                dct8x8(block, coeffs);
+                for (int i = 0; i < 8; ++i)
+                    for (int j = 0; j < 8; ++j)
+                        if (i != 0 || j != 0) record(coeffs[i][j], m_lumaQuantTable[i][j]);
+            }
+    }
+
+    if (total > 0) {
+        for (double& h : stats.histogram) h /= static_cast<double>(total);
+        stats.zeroFraction = static_cast<double>(zeroed) / static_cast<double>(total);
+    }
+    return stats;
+}
+
 Image ImageCodec::subbandMap(const Image& channel) const
 {
     if (!isWavelet()) return Image();