diff --git a/Makefile b/Makefile index 978e993..e074abb 100644 --- a/Makefile +++ b/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 diff --git a/core/inc/ImageCodec.h b/core/inc/ImageCodec.h index 0e93f87..947e050 100644 --- a/core/inc/ImageCodec.h +++ b/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). */ diff --git a/core/src/CodecAnalysis.cpp b/core/src/CodecAnalysis.cpp index e403893..e3f023c 100644 --- a/core/src/CodecAnalysis.cpp +++ b/core/src/CodecAnalysis.cpp @@ -19,6 +19,7 @@ #include #include "colorspace.h" #include +#include Image CodecAnalysis::computeArtifactMap( const Image& original, @@ -136,62 +137,89 @@ 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 ssimKernel() { + std::vector 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& 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 gaussianBlur(const std::vector& src, int w, int h, + const std::vector& k) { + std::vector 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& a, const std::vector& 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 k = ssimKernel(); + + std::vector 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(count) : 1.0; } double CodecAnalysis::computeSSIM(const Image& I1, const Image& I2) { @@ -199,48 +227,19 @@ double CodecAnalysis::computeSSIM(const Image& I1, const Image& I2) { 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(w) * h; + double total = 0.0; + std::vector 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) { diff --git a/core/src/ImageCodec.cpp b/core/src/ImageCodec.cpp index ea0f33b..de9da69 100644 --- a/core/src/ImageCodec.cpp +++ b/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(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 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(total); + stats.zeroFraction = static_cast(zeroed) / static_cast(total); + } + return stats; +} + Image ImageCodec::subbandMap(const Image& channel) const { if (!isWavelet()) return Image(); diff --git a/core/src/transform.cpp b/core/src/transform.cpp index da79d48..1f0b645 100644 --- a/core/src/transform.cpp +++ b/core/src/transform.cpp @@ -20,52 +20,73 @@ const double PI = 3.14159265358979323846; - /* -* Helper function to compute the normalization factor C(u) used in DCT and IDCT. +* Orthonormal 8-point DCT-II basis: B[u][x] = 1/2 · C(u) · cos((2x+1)uπ/16), +* with C(0) = 1/√2 and C(u) = 1 otherwise. Computed once on first use. +* +* The 2D DCT is separable, F = B · f · Bᵀ, and the IDCT is f = Bᵀ · F · B. +* Applying it as two 1D passes costs 2·8³ = 1024 multiply-adds per block, +* instead of 8⁴ = 4096 terms with two cos() calls each for the direct formula +* F(u,v) = 1/4 C(u) C(v) Σx Σy f(x,y) cos((2x+1)uπ/16) cos((2y+1)vπ/16). */ -static double C(int u) { - return (u == 0) ? (1.0 / std::sqrt(2.0)) : 1.0; +struct DctBasis { + double m[8][8]; + DctBasis() { + for (int u = 0; u < 8; ++u) { + const double c = (u == 0) ? (1.0 / std::sqrt(2.0)) : 1.0; + for (int x = 0; x < 8; ++x) + m[u][x] = 0.5 * c * std::cos(((2 * x + 1) * u * PI) / 16.0); + } + } +}; + +static const double (&dctBasis())[8][8] { + static const DctBasis basis; // thread-safe one-time initialisation + return basis.m; } /* -* Performs the Discrete Cosine Transform (DCT) on an 8x8 block. -* The input block should be of type CV_64F and the output will also be CV_64F. -* 2D DCT is computed using the formula: -* F(u, v) = 1/4 * C(u) * C(v) * sum_{x=0}^{7} sum_{y=0}^{7} f(x, y) * cos((2x+1)uπ/16) * cos((2y+1)vπ/16) +* Performs the 2D Discrete Cosine Transform (DCT-II) on an 8x8 block. +* src is indexed [x][y] (row, column); dst is indexed [u][v]. */ void dct8x8(const double src[8][8], double dst[8][8]) { - for (int u = 0; u < 8; ++u) { + const double (&B)[8][8] = dctBasis(); + double tmp[8][8]; // tmp = B · src (transform along the first index) + + for (int u = 0; u < 8; ++u) + for (int y = 0; y < 8; ++y) { + double sum = 0.0; + for (int x = 0; x < 8; ++x) sum += B[u][x] * src[x][y]; + tmp[u][y] = sum; + } + + for (int u = 0; u < 8; ++u) for (int v = 0; v < 8; ++v) { double sum = 0.0; - for (int x = 0; x < 8; ++x) { - for (int y = 0; y < 8; ++y) { - sum += src[x][y] * - std::cos(((2 * x + 1) * u * PI) / 16.0) * - std::cos(((2 * y + 1) * v * PI) / 16.0); - } - } - dst[u][v] = 0.25 * C(u) * C(v) * sum; + for (int y = 0; y < 8; ++y) sum += tmp[u][y] * B[v][y]; + dst[u][v] = sum; } - } } /* * Performs the Inverse Discrete Cosine Transform (IDCT) on an 8x8 block. -* The input block should be of type CV_64F and the output will also be CV_64F. +* src is indexed [u][v]; dst is indexed [x][y]. */ void idct8x8(const double src[8][8], double dst[8][8]) { - for (int x = 0; x < 8; ++x) { + const double (&B)[8][8] = dctBasis(); + double tmp[8][8]; // tmp = Bᵀ · src + + for (int x = 0; x < 8; ++x) + for (int v = 0; v < 8; ++v) { + double sum = 0.0; + for (int u = 0; u < 8; ++u) sum += B[u][x] * src[u][v]; + tmp[x][v] = sum; + } + + for (int x = 0; x < 8; ++x) for (int y = 0; y < 8; ++y) { double sum = 0.0; - for (int u = 0; u < 8; ++u) { - for (int v = 0; v < 8; ++v) { - sum += C(u) * C(v) * src[u][v] * - std::cos(((2 * x + 1) * u * PI) / 16.0) * - std::cos(((2 * y + 1) * v * PI) / 16.0); - } - } - dst[x][y] = 0.25 * sum; + for (int v = 0; v < 8; ++v) sum += tmp[x][v] * B[v][y]; + dst[x][y] = sum; } - } } diff --git a/core/tests/test_codecanalysis.cpp b/core/tests/test_codecanalysis.cpp index 387bdea..5842759 100644 --- a/core/tests/test_codecanalysis.cpp +++ b/core/tests/test_codecanalysis.cpp @@ -108,3 +108,53 @@ TEST(CodecAnalysisTest, StructuralMapsAcceptSingleChannel) { EXPECT_GT(CodecAnalysis::computeEdgeDistortionMap(a, b).at(8, 8, 0), 0.0); EXPECT_GT(CodecAnalysis::computeBlockingMap(b).at(8, 8, 0), 0.0); } + +// --- SSIM against an independent reference ------------------------------- +// Expected values come from a NumPy implementation of Wang et al. (2004): +// 11×11 Gaussian window, σ = 1.5, K1 = 0.01, K2 = 0.03, L = 255, mean over +// the 'valid' region — the formulation of ssim_index.m and scikit-image. + +static Image lcgImage(int w, int h, unsigned seed) { + Image img(w, h, 1); + unsigned s = seed; + for (int y = 0; y < h; ++y) + for (int x = 0; x < w; ++x) { + s = s * 1103515245u + 12345u; + img.at(x, y, 0) = static_cast((s >> 16) % 256); + } + return img; +} + +static Image smoothImage(int w, int h) { + Image img(w, h, 1); + for (int y = 0; y < h; ++y) + for (int x = 0; x < w; ++x) + img.at(x, y, 0) = 128.0 + 90.0 * std::sin(x * 0.21) * std::cos(y * 0.17); + return img; +} + +TEST(CodecAnalysisTest, SSIM_MatchesReferenceImplementation) { + const int W = 40, H = 32; + const Image noise = lcgImage(W, H, 12345u); + const Image smooth = smoothImage(W, H); + + Image noisy = smooth, shifted = smooth; + for (int y = 0; y < H; ++y) + for (int x = 0; x < W; ++x) { + noisy.at(x, y, 0) += (noise.at(x, y, 0) - 128.0) / 8.0; + shifted.at(x, y, 0) += 20.0; + } + + EXPECT_NEAR(CodecAnalysis::computeSSIM(smooth, noisy), 0.868245474080, 1e-9); + EXPECT_NEAR(CodecAnalysis::computeSSIM(noise, smooth), 0.007581148316, 1e-9); + EXPECT_NEAR(CodecAnalysis::computeSSIM(smooth, shifted), 0.985743727949, 1e-9); +} + +TEST(CodecAnalysisTest, SSIM_SmallImagesFallBackGracefully) { + // Smaller than the 11×11 window: clamped borders, still in range. + Image a = smoothImage(6, 5), b = lcgImage(6, 5, 7u); + const double ssim = CodecAnalysis::computeSSIM(a, b); + EXPECT_GE(ssim, -1.0); + EXPECT_LE(ssim, 1.0); + EXPECT_NEAR(CodecAnalysis::computeSSIM(a, a), 1.0, 1e-12); +} diff --git a/core/tests/test_imagecodec.cpp b/core/tests/test_imagecodec.cpp index 23996c8..fc91c22 100644 --- a/core/tests/test_imagecodec.cpp +++ b/core/tests/test_imagecodec.cpp @@ -440,3 +440,40 @@ TEST(ImageCodecTest, InspectBlockPadsPartialEdgeBlocksByReplication) { EXPECT_DOUBLE_EQ(debug.original[i][j], channel.at(x, y, 0)) << i << "," << j; } } + +TEST(ImageCodecTest, CoefficientStatsAreNormalisedAndTrackQuality) { + Image noise = createNoiseImage(64, 48); + Image luma(64, 48, 1); + for (int y = 0; y < 48; ++y) + for (int x = 0; x < 64; ++x) + luma.at(x, y, 0) = noise.at(x, y, 0); + + for (auto t : ALL_TRANSFORMS) { + ImageCodec low(10.0, true, ImageCodec::ChromaSubsampling::CS_444, t); + ImageCodec high(90.0, true, ImageCodec::ChromaSubsampling::CS_444, t); + const auto lo = low.coefficientStats(luma); + const auto hi = high.coefficientStats(luma); + + double sum = 0.0; + for (double h : lo.histogram) sum += h; + EXPECT_NEAR(sum, 1.0, 1e-9) << "transform " << static_cast(t); + + // Same transform → same histogram; coarser quantization zeroes more. + for (int b = 0; b < ImageCodec::CoefficientStats::NUM_BINS; ++b) + EXPECT_DOUBLE_EQ(lo.histogram[b], hi.histogram[b]); + EXPECT_GT(lo.zeroFraction, hi.zeroFraction); + } +} + +TEST(ImageCodecTest, CoefficientStatsWithoutQuantizationCountOnlyExactZeros) { + Image flat(16, 16, 1); + for (size_t i = 0; i < flat.size(); ++i) flat.data()[i] = 77.0; + ImageCodec codec(50.0, false); + // A flat block has all-zero AC coefficients (up to rounding). + EXPECT_GT(codec.coefficientStats(flat).histogram[0], 0.99); + + Image noise(16, 16, 1); + unsigned st = 3u; + for (size_t i = 0; i < noise.size(); ++i) { st = st * 1103515245u + 12345u; noise.data()[i] = (st >> 16) % 256; } + EXPECT_LT(codec.coefficientStats(noise).zeroFraction, 0.05); +} diff --git a/core/tests/transform_test.cpp b/core/tests/transform_test.cpp index 54f06b3..efc65c9 100644 --- a/core/tests/transform_test.cpp +++ b/core/tests/transform_test.cpp @@ -56,3 +56,47 @@ TEST(TransformTest, RoundTrip) { } } } + +// Reference: the direct (non-separable) DCT-II definition the fast version +// must match exactly (up to floating-point rounding). +static void referenceDct(const double src[8][8], double dst[8][8]) { + const double pi = 3.14159265358979323846; + for (int u = 0; u < 8; ++u) + for (int v = 0; v < 8; ++v) { + double sum = 0.0; + for (int x = 0; x < 8; ++x) + for (int y = 0; y < 8; ++y) + sum += src[x][y] * std::cos((2 * x + 1) * u * pi / 16.0) + * std::cos((2 * y + 1) * v * pi / 16.0); + const double cu = u == 0 ? 1.0 / std::sqrt(2.0) : 1.0; + const double cv = v == 0 ? 1.0 / std::sqrt(2.0) : 1.0; + dst[u][v] = 0.25 * cu * cv * sum; + } +} + +TEST(TransformTest, MatchesDirectDefinition) { + double src[8][8], fast[8][8], ref[8][8]; + unsigned state = 7u; + for (int i = 0; i < 8; ++i) + for (int j = 0; j < 8; ++j) { + state = state * 1103515245u + 12345u; + src[i][j] = static_cast((state >> 16) % 256) - 128.0; + } + dct8x8(src, fast); + referenceDct(src, ref); + for (int u = 0; u < 8; ++u) + for (int v = 0; v < 8; ++v) + ASSERT_NEAR(fast[u][v], ref[u][v], 1e-9) << u << "," << v; +} + +TEST(TransformTest, AsymmetricBlockKeepsOrientation) { + // A horizontal ramp (varies along the column index y) must put its + // energy in row u = 0, i.e. dst[0][v], not in column dst[u][0]. + double src[8][8], dst[8][8]; + for (int x = 0; x < 8; ++x) + for (int y = 0; y < 8; ++y) + src[x][y] = y * 10.0; + dct8x8(src, dst); + EXPECT_GT(std::abs(dst[0][1]), 1.0); + EXPECT_NEAR(dst[1][0], 0.0, 1e-9); +} diff --git a/readme.md b/readme.md index 2d8d215..afb6506 100644 --- a/readme.md +++ b/readme.md @@ -11,10 +11,12 @@ An interactive codec laboratory to visualize how transform-based image compressi - DCT vs. wavelet (Haar and CDF 9/7, the JPEG 2000 wavelet) comparison - Wavelet subband pyramid view - Rate–distortion curves using a per-subband entropy estimate of bits -- Interactive block inspector +- Interactive block inspector, with "what if?" coefficient dropping (keep the first N in zig-zag order, or toggle individual coefficients) - Zig-zag and zero-run visualization -- Artifact analysis maps -- PSNR-based quality metrics +- Artifact analysis maps — click a hotspot to open that block in the inspector +- PSNR and SSIM (standard 11×11 Gaussian window) quality metrics +- Coefficient sparsity comparison across DCT, Haar and CDF 9/7 +- Live preview while dragging the quality slider; keyboard- and screen-reader-accessible inspectors; phone-friendly layout - Shareable settings via the URL (`#q=30&t=cdf97&cs=420&img=3`), PNG export, and keyboard shortcuts (`[` / `]` quality, `T` transform) ## 🛠 Technologies diff --git a/web/cpp/codec_web.cpp b/web/cpp/codec_web.cpp index 8673083..514a9d8 100644 --- a/web/cpp/codec_web.cpp +++ b/web/cpp/codec_web.cpp @@ -95,22 +95,34 @@ void init_session(uint8_t* rgba_input, int width, int height) { g_session.initialized = true; } -EMSCRIPTEN_KEEPALIVE -void process_image(int quality, int cs_mode, int transform_mode) { - if (!g_session.initialized) return; - - auto cs = map_cs_mode(cs_mode); - auto transform = map_transform_mode(transform_mode); - ImageCodec codec(quality, true, cs, transform); +// Encodes the session image and stores the reconstruction (no metrics). +static void encodeIntoSession(int quality, int cs_mode, int transform_mode) { + ImageCodec codec(quality, true, map_cs_mode(cs_mode), map_transform_mode(transform_mode)); g_session.processedBgr = codec.process(g_session.originalImage); g_session.lastBitEstimate = codec.getLastBitEstimate(); - g_session.metrics = CodecAnalysis::computeMetrics(g_session.originalImage, g_session.processedBgr); g_session.processedYCrCb = bgrToYCrCb(g_session.processedBgr); g_session.quality = quality; g_session.csMode = cs_mode; g_session.transformMode = transform_mode; } +EMSCRIPTEN_KEEPALIVE +void process_image(int quality, int cs_mode, int transform_mode) { + if (!g_session.initialized) return; + + encodeIntoSession(quality, cs_mode, transform_mode); + g_session.metrics = CodecAnalysis::computeMetrics(g_session.originalImage, g_session.processedBgr); +} + +// Fast path for live preview while a slider is dragged: updates the displayed +// reconstruction and bit estimate but skips PSNR/SSIM (left from the last +// process_image call). Follow up with process_image once the value settles. +EMSCRIPTEN_KEEPALIVE +void preview_image(int quality, int cs_mode, int transform_mode) { + if (!g_session.initialized) return; + encodeIntoSession(quality, cs_mode, transform_mode); +} + // Lightweight encode for rate–distortion sweeps: returns luma PSNR and records // the bit estimate (read via get_rd_bit_estimate) without touching the // session's processed image or computing SSIM. @@ -271,65 +283,28 @@ double get_last_bit_estimate() { return g_session.initialized ? g_session.lastBitEstimate : 0.0; } -// Returns a malloc'd array of 2*num_bins doubles: [dct_bins | dwt_bins]. -// Each series is normalized by the total AC coefficient count. -// Scans all 8×8 blocks of the Y channel, applies both DCT and DWT per block, -// and bins the absolute magnitudes of AC coefficients (skipping DC [0][0]). -// Caller must free the returned pointer. +// Coefficient statistics of the original luma for every transform at the +// given quality, computed with the codec's own transform and quantizer. +// Returns a malloc'd array of 3 × (NUM_BINS + 1) doubles, one record per +// transform mode (0 = DCT, 1 = Haar, 2 = CDF 9/7): +// [histogram[0..NUM_BINS-1], zeroFraction] +// See ImageCodec::CoefficientStats for the bin layout. Caller must free. EMSCRIPTEN_KEEPALIVE -double* get_coeff_histogram(int num_bins, double max_val) { - if (!g_session.initialized || num_bins <= 0 || max_val <= 0.0) return nullptr; - - const Image& ycrcb = g_session.originalYCrCb; - const int w = ycrcb.width(); - const int h = ycrcb.height(); - const int blocksX = w / 8; - const int blocksY = h / 8; - if (blocksX == 0 || blocksY == 0) return nullptr; +double* get_coeff_stats(int quality) { + if (!g_session.initialized) return nullptr; - double* out = (double*)calloc(2 * num_bins, sizeof(double)); + constexpr int BINS = ImageCodec::CoefficientStats::NUM_BINS; + constexpr int RECORD = BINS + 1; + double* out = (double*)malloc(3 * RECORD * sizeof(double)); if (!out) return nullptr; - double* dctBins = out; - double* dwtBins = out + num_bins; - long long totalCount = 0; - - double src[8][8], dctCoeffs[8][8], dwtCoeffs[8][8]; - const double* ycrcbData = ycrcb.data(); - - for (int by = 0; by < blocksY; by++) { - for (int bx = 0; bx < blocksX; bx++) { - // Extract 8×8 Y block with level shift matching processChannel - for (int row = 0; row < 8; row++) { - for (int col = 0; col < 8; col++) { - src[row][col] = ycrcbData[((by * 8 + row) * w + (bx * 8 + col)) * 3] - 128.0; - } - } - - dct8x8(src, dctCoeffs); - dwt8x8(src, dwtCoeffs); - - // Bin AC coefficients — skip [0][0] (DC / LL) - for (int i = 0; i < 8; i++) { - for (int j = 0; j < 8; j++) { - if (i == 0 && j == 0) continue; - - int dctBin = std::min(num_bins - 1, (int)(std::abs(dctCoeffs[i][j]) / max_val * num_bins)); - dctBins[dctBin]++; - - int dwtBin = std::min(num_bins - 1, (int)(std::abs(dwtCoeffs[i][j]) / max_val * num_bins)); - dwtBins[dwtBin]++; - - totalCount++; - } - } - } - } - - if (totalCount > 0) { - for (int i = 0; i < 2 * num_bins; i++) out[i] /= totalCount; + const Image luma = extractChannel(g_session.originalYCrCb, 0); + for (int mode = 0; mode < 3; ++mode) { + ImageCodec codec(quality, true, ImageCodec::ChromaSubsampling::CS_444, map_transform_mode(mode)); + const auto stats = codec.coefficientStats(luma); + std::copy(stats.histogram, stats.histogram + BINS, out + mode * RECORD); + out[mode * RECORD + BINS] = stats.zeroFraction; } - return out; } diff --git a/web/public/codec.js b/web/public/codec.js index 85a8c67..2668789 100644 --- a/web/public/codec.js +++ b/web/public/codec.js @@ -1 +1 @@ -var Module=typeof Module!="undefined"?Module:{};var ENVIRONMENT_IS_WEB=!!globalThis.window;var ENVIRONMENT_IS_WORKER=!!globalThis.WorkerGlobalScope;var ENVIRONMENT_IS_NODE=globalThis.process?.versions?.node&&globalThis.process?.type!="renderer";var arguments_=[];var thisProgram="./this.program";var quit_=(status,toThrow)=>{throw toThrow};var _scriptName=globalThis.document?.currentScript?.src;if(typeof __filename!="undefined"){_scriptName=__filename}else if(ENVIRONMENT_IS_WORKER){_scriptName=self.location.href}var scriptDirectory="";function locateFile(path){if(Module["locateFile"]){return Module["locateFile"](path,scriptDirectory)}return scriptDirectory+path}var readAsync,readBinary;if(ENVIRONMENT_IS_NODE){var fs=require("node:fs");scriptDirectory=__dirname+"/";readBinary=filename=>{filename=isFileURI(filename)?new URL(filename):filename;var ret=fs.readFileSync(filename);return ret};readAsync=async(filename,binary=true)=>{filename=isFileURI(filename)?new URL(filename):filename;var ret=fs.readFileSync(filename,binary?undefined:"utf8");return ret};if(process.argv.length>1){thisProgram=process.argv[1].replace(/\\/g,"/")}arguments_=process.argv.slice(2);if(typeof module!="undefined"){module["exports"]=Module}quit_=(status,toThrow)=>{process.exitCode=status;throw toThrow}}else if(ENVIRONMENT_IS_WEB||ENVIRONMENT_IS_WORKER){try{scriptDirectory=new URL(".",_scriptName).href}catch{}{if(ENVIRONMENT_IS_WORKER){readBinary=url=>{var xhr=new XMLHttpRequest;xhr.open("GET",url,false);xhr.responseType="arraybuffer";xhr.send(null);return new Uint8Array(xhr.response)}}readAsync=async url=>{if(isFileURI(url)){return new Promise((resolve,reject)=>{var xhr=new XMLHttpRequest;xhr.open("GET",url,true);xhr.responseType="arraybuffer";xhr.onload=()=>{if(xhr.status==200||xhr.status==0&&xhr.response){resolve(xhr.response);return}reject(xhr.status)};xhr.onerror=reject;xhr.send(null)})}var response=await fetch(url,{credentials:"same-origin"});if(response.ok){return response.arrayBuffer()}throw new Error(response.status+" : "+response.url)}}}else{}var out=console.log.bind(console);var err=console.error.bind(console);var wasmBinary;var ABORT=false;var isFileURI=filename=>filename.startsWith("file://");var HEAP8,HEAPU8,HEAP16,HEAPU16,HEAP32,HEAPU32,HEAPF32,HEAPF64;var HEAP64,HEAPU64;var runtimeInitialized=false;function updateMemoryViews(){var b=wasmMemory.buffer;HEAP8=new Int8Array(b);HEAP16=new Int16Array(b);Module["HEAPU8"]=HEAPU8=new Uint8Array(b);HEAPU16=new Uint16Array(b);HEAP32=new Int32Array(b);HEAPU32=new Uint32Array(b);HEAPF32=new Float32Array(b);HEAPF64=new Float64Array(b);HEAP64=new BigInt64Array(b);HEAPU64=new BigUint64Array(b)}function preRun(){if(Module["preRun"]){if(typeof Module["preRun"]=="function")Module["preRun"]=[Module["preRun"]];while(Module["preRun"].length){addOnPreRun(Module["preRun"].shift())}}callRuntimeCallbacks(onPreRuns)}function initRuntime(){runtimeInitialized=true;wasmExports["e"]()}function postRun(){if(Module["postRun"]){if(typeof Module["postRun"]=="function")Module["postRun"]=[Module["postRun"]];while(Module["postRun"].length){addOnPostRun(Module["postRun"].shift())}}callRuntimeCallbacks(onPostRuns)}function abort(what){Module["onAbort"]?.(what);what="Aborted("+what+")";err(what);ABORT=true;what+=". Build with -sASSERTIONS for more info.";var e=new WebAssembly.RuntimeError(what);throw e}var wasmBinaryFile;function findWasmBinary(){return locateFile("codec.wasm")}function getBinarySync(file){if(file==wasmBinaryFile&&wasmBinary){return new Uint8Array(wasmBinary)}if(readBinary){return readBinary(file)}throw"both async and sync fetching of the wasm failed"}async function getWasmBinary(binaryFile){if(!wasmBinary){try{var response=await readAsync(binaryFile);return new Uint8Array(response)}catch{}}return getBinarySync(binaryFile)}async function instantiateArrayBuffer(binaryFile,imports){try{var binary=await getWasmBinary(binaryFile);var instance=await WebAssembly.instantiate(binary,imports);return instance}catch(reason){err(`failed to asynchronously prepare wasm: ${reason}`);abort(reason)}}async function instantiateAsync(binary,binaryFile,imports){if(!binary&&!isFileURI(binaryFile)&&!ENVIRONMENT_IS_NODE){try{var response=fetch(binaryFile,{credentials:"same-origin"});var instantiationResult=await WebAssembly.instantiateStreaming(response,imports);return instantiationResult}catch(reason){err(`wasm streaming compile failed: ${reason}`);err("falling back to ArrayBuffer instantiation")}}return instantiateArrayBuffer(binaryFile,imports)}function getWasmImports(){var imports={a:wasmImports};return imports}async function createWasm(){function receiveInstance(instance,module){wasmExports=instance.exports;assignWasmExports(wasmExports);updateMemoryViews();removeRunDependency("wasm-instantiate");return wasmExports}addRunDependency("wasm-instantiate");function receiveInstantiationResult(result){return receiveInstance(result["instance"])}var info=getWasmImports();if(Module["instantiateWasm"]){return new Promise((resolve,reject)=>{Module["instantiateWasm"](info,(inst,mod)=>{resolve(receiveInstance(inst,mod))})})}wasmBinaryFile??=findWasmBinary();var result=await instantiateAsync(wasmBinary,wasmBinaryFile,info);var exports=receiveInstantiationResult(result);return exports}class ExitStatus{name="ExitStatus";constructor(status){this.message=`Program terminated with exit(${status})`;this.status=status}}var callRuntimeCallbacks=callbacks=>{while(callbacks.length>0){callbacks.shift()(Module)}};var onPostRuns=[];var addOnPostRun=cb=>onPostRuns.push(cb);var onPreRuns=[];var addOnPreRun=cb=>onPreRuns.push(cb);var runDependencies=0;var dependenciesFulfilled=null;var removeRunDependency=id=>{runDependencies--;Module["monitorRunDependencies"]?.(runDependencies);if(runDependencies==0){if(dependenciesFulfilled){var callback=dependenciesFulfilled;dependenciesFulfilled=null;callback()}}};var addRunDependency=id=>{runDependencies++;Module["monitorRunDependencies"]?.(runDependencies)};var noExitRuntime=true;var stackRestore=val=>__emscripten_stack_restore(val);var stackSave=()=>_emscripten_stack_get_current();class ExceptionInfo{constructor(excPtr){this.excPtr=excPtr;this.ptr=excPtr-24}set_type(type){HEAPU32[this.ptr+4>>2]=type}get_type(){return HEAPU32[this.ptr+4>>2]}set_destructor(destructor){HEAPU32[this.ptr+8>>2]=destructor}get_destructor(){return HEAPU32[this.ptr+8>>2]}set_caught(caught){caught=caught?1:0;HEAP8[this.ptr+12]=caught}get_caught(){return HEAP8[this.ptr+12]!=0}set_rethrown(rethrown){rethrown=rethrown?1:0;HEAP8[this.ptr+13]=rethrown}get_rethrown(){return HEAP8[this.ptr+13]!=0}init(type,destructor){this.set_adjusted_ptr(0);this.set_type(type);this.set_destructor(destructor)}set_adjusted_ptr(adjustedPtr){HEAPU32[this.ptr+16>>2]=adjustedPtr}get_adjusted_ptr(){return HEAPU32[this.ptr+16>>2]}}var exceptionLast=0;var uncaughtExceptionCount=0;var ___cxa_throw=(ptr,type,destructor)=>{var info=new ExceptionInfo(ptr);info.init(type,destructor);exceptionLast=ptr;uncaughtExceptionCount++;throw exceptionLast};var __abort_js=()=>abort("");var getHeapMax=()=>2147483648;var alignMemory=(size,alignment)=>Math.ceil(size/alignment)*alignment;var growMemory=size=>{var oldHeapSize=wasmMemory.buffer.byteLength;var pages=(size-oldHeapSize+65535)/65536|0;try{wasmMemory.grow(pages);updateMemoryViews();return 1}catch(e){}};var _emscripten_resize_heap=requestedSize=>{var oldSize=HEAPU8.length;requestedSize>>>=0;var maxHeapSize=getHeapMax();if(requestedSize>maxHeapSize){return false}for(var cutDown=1;cutDown<=4;cutDown*=2){var overGrownHeapSize=oldSize*(1+.2/cutDown);overGrownHeapSize=Math.min(overGrownHeapSize,requestedSize+100663296);var newSize=Math.min(maxHeapSize,alignMemory(Math.max(requestedSize,overGrownHeapSize),65536));var replacement=growMemory(newSize);if(replacement){return true}}return false};var getCFunc=ident=>{var func=Module["_"+ident];return func};var writeArrayToMemory=(array,buffer)=>{HEAP8.set(array,buffer)};var lengthBytesUTF8=str=>{var len=0;for(var i=0;i=55296&&c<=57343){len+=4;++i}else{len+=3}}return len};var stringToUTF8Array=(str,heap,outIdx,maxBytesToWrite)=>{if(!(maxBytesToWrite>0))return 0;var startIdx=outIdx;var endIdx=outIdx+maxBytesToWrite-1;for(var i=0;i=endIdx)break;heap[outIdx++]=u}else if(u<=2047){if(outIdx+1>=endIdx)break;heap[outIdx++]=192|u>>6;heap[outIdx++]=128|u&63}else if(u<=65535){if(outIdx+2>=endIdx)break;heap[outIdx++]=224|u>>12;heap[outIdx++]=128|u>>6&63;heap[outIdx++]=128|u&63}else{if(outIdx+3>=endIdx)break;heap[outIdx++]=240|u>>18;heap[outIdx++]=128|u>>12&63;heap[outIdx++]=128|u>>6&63;heap[outIdx++]=128|u&63;i++}}heap[outIdx]=0;return outIdx-startIdx};var stringToUTF8=(str,outPtr,maxBytesToWrite)=>stringToUTF8Array(str,HEAPU8,outPtr,maxBytesToWrite);var stackAlloc=sz=>__emscripten_stack_alloc(sz);var stringToUTF8OnStack=str=>{var size=lengthBytesUTF8(str)+1;var ret=stackAlloc(size);stringToUTF8(str,ret,size);return ret};var UTF8Decoder=globalThis.TextDecoder&&new TextDecoder;var findStringEnd=(heapOrArray,idx,maxBytesToRead,ignoreNul)=>{var maxIdx=idx+maxBytesToRead;if(ignoreNul)return maxIdx;while(heapOrArray[idx]&&!(idx>=maxIdx))++idx;return idx};var UTF8ArrayToString=(heapOrArray,idx=0,maxBytesToRead,ignoreNul)=>{var endPtr=findStringEnd(heapOrArray,idx,maxBytesToRead,ignoreNul);if(endPtr-idx>16&&heapOrArray.buffer&&UTF8Decoder){return UTF8Decoder.decode(heapOrArray.subarray(idx,endPtr))}var str="";while(idx>10,56320|ch&1023)}}return str};var UTF8ToString=(ptr,maxBytesToRead,ignoreNul)=>ptr?UTF8ArrayToString(HEAPU8,ptr,maxBytesToRead,ignoreNul):"";var ccall=(ident,returnType,argTypes,args,opts)=>{var toC={string:str=>{var ret=0;if(str!==null&&str!==undefined&&str!==0){ret=stringToUTF8OnStack(str)}return ret},array:arr=>{var ret=stackAlloc(arr.length);writeArrayToMemory(arr,ret);return ret}};function convertReturnValue(ret){if(returnType==="string"){return UTF8ToString(ret)}if(returnType==="boolean")return Boolean(ret);return ret}var func=getCFunc(ident);var cArgs=[];var stack=0;if(args){for(var i=0;i{var numericArgs=!argTypes||argTypes.every(type=>type==="number"||type==="boolean");var numericRet=returnType!=="string";if(numericRet&&numericArgs&&!opts){return getCFunc(ident)}return(...args)=>ccall(ident,returnType,argTypes,args,opts)};{if(Module["noExitRuntime"])noExitRuntime=Module["noExitRuntime"];if(Module["print"])out=Module["print"];if(Module["printErr"])err=Module["printErr"];if(Module["wasmBinary"])wasmBinary=Module["wasmBinary"];if(Module["arguments"])arguments_=Module["arguments"];if(Module["thisProgram"])thisProgram=Module["thisProgram"];if(Module["preInit"]){if(typeof Module["preInit"]=="function")Module["preInit"]=[Module["preInit"]];while(Module["preInit"].length>0){Module["preInit"].shift()()}}}Module["ccall"]=ccall;Module["cwrap"]=cwrap;var _init_session,_process_image,_rd_encode,_get_rd_bit_estimate,_get_view_ptr,_malloc,_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,_get_coeff_histogram,_inspect_block_data,_init_me_session,_run_motion_estimation,_get_mv_count,_get_me_residual_ptr,_get_search_steps,_get_search_step_count,_free,__emscripten_stack_restore,__emscripten_stack_alloc,_emscripten_stack_get_current,memory,__indirect_function_table,wasmMemory;function assignWasmExports(wasmExports){_init_session=Module["_init_session"]=wasmExports["f"];_process_image=Module["_process_image"]=wasmExports["g"];_rd_encode=Module["_rd_encode"]=wasmExports["h"];_get_rd_bit_estimate=Module["_get_rd_bit_estimate"]=wasmExports["i"];_get_view_ptr=Module["_get_view_ptr"]=wasmExports["j"];_malloc=Module["_malloc"]=wasmExports["k"];_set_view_tint=Module["_set_view_tint"]=wasmExports["l"];_set_artifact_gain=Module["_set_artifact_gain"]=wasmExports["m"];_get_psnr_y=Module["_get_psnr_y"]=wasmExports["n"];_get_psnr_cr=Module["_get_psnr_cr"]=wasmExports["o"];_get_psnr_cb=Module["_get_psnr_cb"]=wasmExports["p"];_get_ssim_y=Module["_get_ssim_y"]=wasmExports["q"];_get_ssim_cr=Module["_get_ssim_cr"]=wasmExports["r"];_get_ssim_cb=Module["_get_ssim_cb"]=wasmExports["s"];_get_last_bit_estimate=Module["_get_last_bit_estimate"]=wasmExports["t"];_get_coeff_histogram=Module["_get_coeff_histogram"]=wasmExports["u"];_inspect_block_data=Module["_inspect_block_data"]=wasmExports["v"];_init_me_session=Module["_init_me_session"]=wasmExports["w"];_run_motion_estimation=Module["_run_motion_estimation"]=wasmExports["x"];_get_mv_count=Module["_get_mv_count"]=wasmExports["y"];_get_me_residual_ptr=Module["_get_me_residual_ptr"]=wasmExports["z"];_get_search_steps=Module["_get_search_steps"]=wasmExports["A"];_get_search_step_count=Module["_get_search_step_count"]=wasmExports["B"];_free=Module["_free"]=wasmExports["C"];__emscripten_stack_restore=wasmExports["D"];__emscripten_stack_alloc=wasmExports["E"];_emscripten_stack_get_current=wasmExports["F"];memory=wasmMemory=wasmExports["d"];__indirect_function_table=wasmExports["__indirect_function_table"]}var wasmImports={a:___cxa_throw,b:__abort_js,c:_emscripten_resize_heap};function run(){if(runDependencies>0){dependenciesFulfilled=run;return}preRun();if(runDependencies>0){dependenciesFulfilled=run;return}function doRun(){Module["calledRun"]=true;if(ABORT)return;initRuntime();Module["onRuntimeInitialized"]?.();postRun()}if(Module["setStatus"]){Module["setStatus"]("Running...");setTimeout(()=>{setTimeout(()=>Module["setStatus"](""),1);doRun()},1)}else{doRun()}}var wasmExports;createWasm();run(); +var Module=typeof Module!="undefined"?Module:{};var ENVIRONMENT_IS_WEB=!!globalThis.window;var ENVIRONMENT_IS_WORKER=!!globalThis.WorkerGlobalScope;var ENVIRONMENT_IS_NODE=globalThis.process?.versions?.node&&globalThis.process?.type!="renderer";var arguments_=[];var thisProgram="./this.program";var quit_=(status,toThrow)=>{throw toThrow};var _scriptName=globalThis.document?.currentScript?.src;if(typeof __filename!="undefined"){_scriptName=__filename}else if(ENVIRONMENT_IS_WORKER){_scriptName=self.location.href}var scriptDirectory="";function locateFile(path){if(Module["locateFile"]){return Module["locateFile"](path,scriptDirectory)}return scriptDirectory+path}var readAsync,readBinary;if(ENVIRONMENT_IS_NODE){var fs=require("node:fs");scriptDirectory=__dirname+"/";readBinary=filename=>{filename=isFileURI(filename)?new URL(filename):filename;var ret=fs.readFileSync(filename);return ret};readAsync=async(filename,binary=true)=>{filename=isFileURI(filename)?new URL(filename):filename;var ret=fs.readFileSync(filename,binary?undefined:"utf8");return ret};if(process.argv.length>1){thisProgram=process.argv[1].replace(/\\/g,"/")}arguments_=process.argv.slice(2);if(typeof module!="undefined"){module["exports"]=Module}quit_=(status,toThrow)=>{process.exitCode=status;throw toThrow}}else if(ENVIRONMENT_IS_WEB||ENVIRONMENT_IS_WORKER){try{scriptDirectory=new URL(".",_scriptName).href}catch{}{if(ENVIRONMENT_IS_WORKER){readBinary=url=>{var xhr=new XMLHttpRequest;xhr.open("GET",url,false);xhr.responseType="arraybuffer";xhr.send(null);return new Uint8Array(xhr.response)}}readAsync=async url=>{if(isFileURI(url)){return new Promise((resolve,reject)=>{var xhr=new XMLHttpRequest;xhr.open("GET",url,true);xhr.responseType="arraybuffer";xhr.onload=()=>{if(xhr.status==200||xhr.status==0&&xhr.response){resolve(xhr.response);return}reject(xhr.status)};xhr.onerror=reject;xhr.send(null)})}var response=await fetch(url,{credentials:"same-origin"});if(response.ok){return response.arrayBuffer()}throw new Error(response.status+" : "+response.url)}}}else{}var out=console.log.bind(console);var err=console.error.bind(console);var wasmBinary;var ABORT=false;var isFileURI=filename=>filename.startsWith("file://");var HEAP8,HEAPU8,HEAP16,HEAPU16,HEAP32,HEAPU32,HEAPF32,HEAPF64;var HEAP64,HEAPU64;var runtimeInitialized=false;function updateMemoryViews(){var b=wasmMemory.buffer;HEAP8=new Int8Array(b);HEAP16=new Int16Array(b);Module["HEAPU8"]=HEAPU8=new Uint8Array(b);HEAPU16=new Uint16Array(b);HEAP32=new Int32Array(b);HEAPU32=new Uint32Array(b);HEAPF32=new Float32Array(b);HEAPF64=new Float64Array(b);HEAP64=new BigInt64Array(b);HEAPU64=new BigUint64Array(b)}function preRun(){if(Module["preRun"]){if(typeof Module["preRun"]=="function")Module["preRun"]=[Module["preRun"]];while(Module["preRun"].length){addOnPreRun(Module["preRun"].shift())}}callRuntimeCallbacks(onPreRuns)}function initRuntime(){runtimeInitialized=true;wasmExports["e"]()}function postRun(){if(Module["postRun"]){if(typeof Module["postRun"]=="function")Module["postRun"]=[Module["postRun"]];while(Module["postRun"].length){addOnPostRun(Module["postRun"].shift())}}callRuntimeCallbacks(onPostRuns)}function abort(what){Module["onAbort"]?.(what);what="Aborted("+what+")";err(what);ABORT=true;what+=". Build with -sASSERTIONS for more info.";var e=new WebAssembly.RuntimeError(what);throw e}var wasmBinaryFile;function findWasmBinary(){return locateFile("codec.wasm")}function getBinarySync(file){if(file==wasmBinaryFile&&wasmBinary){return new Uint8Array(wasmBinary)}if(readBinary){return readBinary(file)}throw"both async and sync fetching of the wasm failed"}async function getWasmBinary(binaryFile){if(!wasmBinary){try{var response=await readAsync(binaryFile);return new Uint8Array(response)}catch{}}return getBinarySync(binaryFile)}async function instantiateArrayBuffer(binaryFile,imports){try{var binary=await getWasmBinary(binaryFile);var instance=await WebAssembly.instantiate(binary,imports);return instance}catch(reason){err(`failed to asynchronously prepare wasm: ${reason}`);abort(reason)}}async function instantiateAsync(binary,binaryFile,imports){if(!binary&&!isFileURI(binaryFile)&&!ENVIRONMENT_IS_NODE){try{var response=fetch(binaryFile,{credentials:"same-origin"});var instantiationResult=await WebAssembly.instantiateStreaming(response,imports);return instantiationResult}catch(reason){err(`wasm streaming compile failed: ${reason}`);err("falling back to ArrayBuffer instantiation")}}return instantiateArrayBuffer(binaryFile,imports)}function getWasmImports(){var imports={a:wasmImports};return imports}async function createWasm(){function receiveInstance(instance,module){wasmExports=instance.exports;assignWasmExports(wasmExports);updateMemoryViews();removeRunDependency("wasm-instantiate");return wasmExports}addRunDependency("wasm-instantiate");function receiveInstantiationResult(result){return receiveInstance(result["instance"])}var info=getWasmImports();if(Module["instantiateWasm"]){return new Promise((resolve,reject)=>{Module["instantiateWasm"](info,(inst,mod)=>{resolve(receiveInstance(inst,mod))})})}wasmBinaryFile??=findWasmBinary();var result=await instantiateAsync(wasmBinary,wasmBinaryFile,info);var exports=receiveInstantiationResult(result);return exports}class ExitStatus{name="ExitStatus";constructor(status){this.message=`Program terminated with exit(${status})`;this.status=status}}var callRuntimeCallbacks=callbacks=>{while(callbacks.length>0){callbacks.shift()(Module)}};var onPostRuns=[];var addOnPostRun=cb=>onPostRuns.push(cb);var onPreRuns=[];var addOnPreRun=cb=>onPreRuns.push(cb);var runDependencies=0;var dependenciesFulfilled=null;var removeRunDependency=id=>{runDependencies--;Module["monitorRunDependencies"]?.(runDependencies);if(runDependencies==0){if(dependenciesFulfilled){var callback=dependenciesFulfilled;dependenciesFulfilled=null;callback()}}};var addRunDependency=id=>{runDependencies++;Module["monitorRunDependencies"]?.(runDependencies)};var noExitRuntime=true;var stackRestore=val=>__emscripten_stack_restore(val);var stackSave=()=>_emscripten_stack_get_current();class ExceptionInfo{constructor(excPtr){this.excPtr=excPtr;this.ptr=excPtr-24}set_type(type){HEAPU32[this.ptr+4>>2]=type}get_type(){return HEAPU32[this.ptr+4>>2]}set_destructor(destructor){HEAPU32[this.ptr+8>>2]=destructor}get_destructor(){return HEAPU32[this.ptr+8>>2]}set_caught(caught){caught=caught?1:0;HEAP8[this.ptr+12]=caught}get_caught(){return HEAP8[this.ptr+12]!=0}set_rethrown(rethrown){rethrown=rethrown?1:0;HEAP8[this.ptr+13]=rethrown}get_rethrown(){return HEAP8[this.ptr+13]!=0}init(type,destructor){this.set_adjusted_ptr(0);this.set_type(type);this.set_destructor(destructor)}set_adjusted_ptr(adjustedPtr){HEAPU32[this.ptr+16>>2]=adjustedPtr}get_adjusted_ptr(){return HEAPU32[this.ptr+16>>2]}}var exceptionLast=0;var uncaughtExceptionCount=0;var ___cxa_throw=(ptr,type,destructor)=>{var info=new ExceptionInfo(ptr);info.init(type,destructor);exceptionLast=ptr;uncaughtExceptionCount++;throw exceptionLast};var __abort_js=()=>abort("");var getHeapMax=()=>2147483648;var alignMemory=(size,alignment)=>Math.ceil(size/alignment)*alignment;var growMemory=size=>{var oldHeapSize=wasmMemory.buffer.byteLength;var pages=(size-oldHeapSize+65535)/65536|0;try{wasmMemory.grow(pages);updateMemoryViews();return 1}catch(e){}};var _emscripten_resize_heap=requestedSize=>{var oldSize=HEAPU8.length;requestedSize>>>=0;var maxHeapSize=getHeapMax();if(requestedSize>maxHeapSize){return false}for(var cutDown=1;cutDown<=4;cutDown*=2){var overGrownHeapSize=oldSize*(1+.2/cutDown);overGrownHeapSize=Math.min(overGrownHeapSize,requestedSize+100663296);var newSize=Math.min(maxHeapSize,alignMemory(Math.max(requestedSize,overGrownHeapSize),65536));var replacement=growMemory(newSize);if(replacement){return true}}return false};var getCFunc=ident=>{var func=Module["_"+ident];return func};var writeArrayToMemory=(array,buffer)=>{HEAP8.set(array,buffer)};var lengthBytesUTF8=str=>{var len=0;for(var i=0;i=55296&&c<=57343){len+=4;++i}else{len+=3}}return len};var stringToUTF8Array=(str,heap,outIdx,maxBytesToWrite)=>{if(!(maxBytesToWrite>0))return 0;var startIdx=outIdx;var endIdx=outIdx+maxBytesToWrite-1;for(var i=0;i=endIdx)break;heap[outIdx++]=u}else if(u<=2047){if(outIdx+1>=endIdx)break;heap[outIdx++]=192|u>>6;heap[outIdx++]=128|u&63}else if(u<=65535){if(outIdx+2>=endIdx)break;heap[outIdx++]=224|u>>12;heap[outIdx++]=128|u>>6&63;heap[outIdx++]=128|u&63}else{if(outIdx+3>=endIdx)break;heap[outIdx++]=240|u>>18;heap[outIdx++]=128|u>>12&63;heap[outIdx++]=128|u>>6&63;heap[outIdx++]=128|u&63;i++}}heap[outIdx]=0;return outIdx-startIdx};var stringToUTF8=(str,outPtr,maxBytesToWrite)=>stringToUTF8Array(str,HEAPU8,outPtr,maxBytesToWrite);var stackAlloc=sz=>__emscripten_stack_alloc(sz);var stringToUTF8OnStack=str=>{var size=lengthBytesUTF8(str)+1;var ret=stackAlloc(size);stringToUTF8(str,ret,size);return ret};var UTF8Decoder=globalThis.TextDecoder&&new TextDecoder;var findStringEnd=(heapOrArray,idx,maxBytesToRead,ignoreNul)=>{var maxIdx=idx+maxBytesToRead;if(ignoreNul)return maxIdx;while(heapOrArray[idx]&&!(idx>=maxIdx))++idx;return idx};var UTF8ArrayToString=(heapOrArray,idx=0,maxBytesToRead,ignoreNul)=>{var endPtr=findStringEnd(heapOrArray,idx,maxBytesToRead,ignoreNul);if(endPtr-idx>16&&heapOrArray.buffer&&UTF8Decoder){return UTF8Decoder.decode(heapOrArray.subarray(idx,endPtr))}var str="";while(idx>10,56320|ch&1023)}}return str};var UTF8ToString=(ptr,maxBytesToRead,ignoreNul)=>ptr?UTF8ArrayToString(HEAPU8,ptr,maxBytesToRead,ignoreNul):"";var ccall=(ident,returnType,argTypes,args,opts)=>{var toC={string:str=>{var ret=0;if(str!==null&&str!==undefined&&str!==0){ret=stringToUTF8OnStack(str)}return ret},array:arr=>{var ret=stackAlloc(arr.length);writeArrayToMemory(arr,ret);return ret}};function convertReturnValue(ret){if(returnType==="string"){return UTF8ToString(ret)}if(returnType==="boolean")return Boolean(ret);return ret}var func=getCFunc(ident);var cArgs=[];var stack=0;if(args){for(var i=0;i{var numericArgs=!argTypes||argTypes.every(type=>type==="number"||type==="boolean");var numericRet=returnType!=="string";if(numericRet&&numericArgs&&!opts){return getCFunc(ident)}return(...args)=>ccall(ident,returnType,argTypes,args,opts)};{if(Module["noExitRuntime"])noExitRuntime=Module["noExitRuntime"];if(Module["print"])out=Module["print"];if(Module["printErr"])err=Module["printErr"];if(Module["wasmBinary"])wasmBinary=Module["wasmBinary"];if(Module["arguments"])arguments_=Module["arguments"];if(Module["thisProgram"])thisProgram=Module["thisProgram"];if(Module["preInit"]){if(typeof Module["preInit"]=="function")Module["preInit"]=[Module["preInit"]];while(Module["preInit"].length>0){Module["preInit"].shift()()}}}Module["ccall"]=ccall;Module["cwrap"]=cwrap;var _init_session,_process_image,_preview_image,_rd_encode,_get_rd_bit_estimate,_get_view_ptr,_malloc,_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,_get_coeff_stats,_inspect_block_data,_init_me_session,_run_motion_estimation,_get_mv_count,_get_me_residual_ptr,_get_search_steps,_get_search_step_count,_free,__emscripten_stack_restore,__emscripten_stack_alloc,_emscripten_stack_get_current,memory,__indirect_function_table,wasmMemory;function assignWasmExports(wasmExports){_init_session=Module["_init_session"]=wasmExports["f"];_process_image=Module["_process_image"]=wasmExports["g"];_preview_image=Module["_preview_image"]=wasmExports["h"];_rd_encode=Module["_rd_encode"]=wasmExports["i"];_get_rd_bit_estimate=Module["_get_rd_bit_estimate"]=wasmExports["j"];_get_view_ptr=Module["_get_view_ptr"]=wasmExports["k"];_malloc=Module["_malloc"]=wasmExports["l"];_set_view_tint=Module["_set_view_tint"]=wasmExports["m"];_set_artifact_gain=Module["_set_artifact_gain"]=wasmExports["n"];_get_psnr_y=Module["_get_psnr_y"]=wasmExports["o"];_get_psnr_cr=Module["_get_psnr_cr"]=wasmExports["p"];_get_psnr_cb=Module["_get_psnr_cb"]=wasmExports["q"];_get_ssim_y=Module["_get_ssim_y"]=wasmExports["r"];_get_ssim_cr=Module["_get_ssim_cr"]=wasmExports["s"];_get_ssim_cb=Module["_get_ssim_cb"]=wasmExports["t"];_get_last_bit_estimate=Module["_get_last_bit_estimate"]=wasmExports["u"];_get_coeff_stats=Module["_get_coeff_stats"]=wasmExports["v"];_inspect_block_data=Module["_inspect_block_data"]=wasmExports["w"];_init_me_session=Module["_init_me_session"]=wasmExports["x"];_run_motion_estimation=Module["_run_motion_estimation"]=wasmExports["y"];_get_mv_count=Module["_get_mv_count"]=wasmExports["z"];_get_me_residual_ptr=Module["_get_me_residual_ptr"]=wasmExports["A"];_get_search_steps=Module["_get_search_steps"]=wasmExports["B"];_get_search_step_count=Module["_get_search_step_count"]=wasmExports["C"];_free=Module["_free"]=wasmExports["D"];__emscripten_stack_restore=wasmExports["E"];__emscripten_stack_alloc=wasmExports["F"];_emscripten_stack_get_current=wasmExports["G"];memory=wasmMemory=wasmExports["d"];__indirect_function_table=wasmExports["__indirect_function_table"]}var wasmImports={a:___cxa_throw,b:__abort_js,c:_emscripten_resize_heap};function run(){if(runDependencies>0){dependenciesFulfilled=run;return}preRun();if(runDependencies>0){dependenciesFulfilled=run;return}function doRun(){Module["calledRun"]=true;if(ABORT)return;initRuntime();Module["onRuntimeInitialized"]?.();postRun()}if(Module["setStatus"]){Module["setStatus"]("Running...");setTimeout(()=>{setTimeout(()=>Module["setStatus"](""),1);doRun()},1)}else{doRun()}}var wasmExports;createWasm();run(); diff --git a/web/public/codec.wasm b/web/public/codec.wasm index c967231..6e72878 100755 Binary files a/web/public/codec.wasm and b/web/public/codec.wasm differ diff --git a/web/src/App.svelte b/web/src/App.svelte index dc970c7..4a0d9e4 100644 --- a/web/src/App.svelte +++ b/web/src/App.svelte @@ -1,5 +1,5 @@ + +
+ +
+ {#each PRESETS as preset} + + {/each} +
+ +
diff --git a/web/src/lib/components/viewer/CodecControls.svelte b/web/src/lib/components/viewer/CodecControls.svelte index 8b2c72d..65850b3 100644 --- a/web/src/lib/components/viewer/CodecControls.svelte +++ b/web/src/lib/components/viewer/CodecControls.svelte @@ -1,6 +1,7 @@ -
- -
- {#each [{ q: 10, label: 'Low', desc: 'Thumbnails' }, { q: 50, label: 'Medium', desc: 'Web / Social' }, { q: 90, label: 'High', desc: 'Print / Archive' }] as preset} - - {/each} -
- -
+
View Mode
diff --git a/web/src/lib/components/viewer/CoeffHistogram.svelte b/web/src/lib/components/viewer/CoeffHistogram.svelte index 48d144e..15549f0 100644 --- a/web/src/lib/components/viewer/CoeffHistogram.svelte +++ b/web/src/lib/components/viewer/CoeffHistogram.svelte @@ -1,37 +1,60 @@ @@ -39,116 +62,97 @@
-

Coefficient Histogram Beta

-

AC coefficient magnitude distribution — DWT concentrates more energy near zero, explaining its advantage at low bitrates.

+

Coefficient Sparsity

+

Where each transform puts the luma channel's detail energy. Coefficients that quantize to zero are nearly free to code.

- {#if histLoading} + {#if loading} Computing… - {:else if histData} - {appState.imgWidth / 8 | 0}×{appState.imgHeight / 8 | 0} blocks {/if} - {#if appState.originalImageData && !histLoading} - + {#if appState.originalImageData && !loading && !stats} + {/if}
{#if !appState.originalImageData} -
Load an image to generate a coefficient histogram.
- {:else if histLoading} -
Computing transform coefficients for all blocks…
- {:else if !histData} -
Click Generate to compare DCT vs DWT coefficient sparsity.
+
Load an image to compare coefficient sparsity.
+ {:else if loading} +
Transforming the image with DCT, Haar and CDF 9/7…
+ {:else if !stats} +
Click Generate to compare how sparse each transform makes the image.
{:else} + {@const series = TRANSFORMS.map(t => ({ info: t, data: stats![t.value] }))} {@const svgW = 700} {@const svgH = 240} {@const padL = 52} - {@const padR = 24} - {@const padT = 18} + {@const padR = 16} + {@const padT = 30} {@const padB = 48} {@const chartW = svgW - padL - padR} {@const chartH = svgH - padT - padB} - {@const maxY = Math.max(...histData.dct, ...histData.dwt) * 1.15} - {@const groupW = chartW / HIST_BINS} - {@const barW = (groupW - 2) / 2} - {@const binWidth = HIST_MAX / HIST_BINS} - {@const mapBarY = (v: number) => padT + chartH - (v / maxY) * chartH} + {@const maxY = Math.max(...series.flatMap(s => s.data.histogram)) * 1.1 || 1} + {@const groupW = chartW / COEFF_BINS} + {@const barW = (groupW - 6) / series.length}
- - - - - - - + {#each [0.25, 0.5, 0.75, 1.0] as frac} {@const gy = padT + chartH * (1 - frac)} - {(maxY * frac * 100).toFixed(1)}% + {(maxY * frac * 100).toFixed(0)}% {/each} - - - {#each histData.dct as val, i} - {@const x = padL + i * groupW + 1} + {#each series as s, k} + {#each s.data.histogram as val, b} {@const barH = Math.max(0, (val / maxY) * chartH)} - + + {s.info.label}: {(val * 100).toFixed(1)}% of coefficients with |c| in bin {binLabel(b)} + {/each} - - - {#each histData.dwt as val, i} - {@const x = padL + i * groupW + 1 + barW} - {@const barH = Math.max(0, (val / maxY) * chartH)} - {/each} - - {#each {length: Math.floor(HIST_BINS / 8) + 1} as _, k} - {@const binIdx = k * 8} - {#if binIdx <= HIST_BINS} - - {binIdx === HIST_BINS ? `${Math.round(HIST_MAX)}+` : Math.round(binIdx * binWidth)} - - {/if} + {#each { length: COEFF_BINS } as _, b} + {binLabel(b)} {/each} - - AC Coefficient Magnitude - Fraction + |AC coefficient| (doubling-width bins) + Share - - - DCT - - DWT (Haar) + {#each series as s, k} + + {s.info.label} + {/each}
- - {@const dctNZ = (histData.dct[0] * 100)} - {@const dwtNZ = (histData.dwt[0] * 100)} - {@const gain = dwtNZ - dctNZ}
-
- DCT near-zero (<{Math.round(binWidth)}) - {dctNZ.toFixed(1)}% -
-
-
- DWT near-zero (<{Math.round(binWidth)}) - {dwtNZ.toFixed(1)}% -
-
-
- DWT sparsity gain - 0} class:stat-poor={gain < 0}> - {gain >= 0 ? '+' : ''}{gain.toFixed(1)} pp - + {#each series as s, k} + {#if k > 0}
{/if} +
+ {s.info.label} zeroed at Q{statsQuality} + {(s.data.zeroFraction * 100).toFixed(1)}%
+ {/each}
+

+ Current transform: {transformInfo(appState.transformType).label}. + DC / LL coefficients are excluded. The zeroed share follows the quality slider. +

{/if}
+ + diff --git a/web/src/lib/components/viewer/ImageSourcePanel.svelte b/web/src/lib/components/viewer/ImageSourcePanel.svelte index 46b7d83..bbfc5b9 100644 --- a/web/src/lib/components/viewer/ImageSourcePanel.svelte +++ b/web/src/lib/components/viewer/ImageSourcePanel.svelte @@ -111,7 +111,7 @@ {#each {length: TEST_IMAGE_COUNT} as _, i} {:else}