CLI: PSNR, NRMSE, and opt-in windowed SSIM

The round-trip report now always includes PSNR and NRMSE (double
accumulation alongside the existing error scan), and --ssim adds the
windowed structural similarity: 7-point windows, stride 2 per axis,
window-local luminance range, averaged over windows, for float32 and
float64 and any recorded shape.

Verified against independent implementations on Miranda at rel 1e-3:
numpy reproduces PSNR 71.39 dB and NRMSE 2.693698e-04 to every digit,
and the SSIM of 0.900335 is identical in all six decimals to the value
computed by cuSZp's own computeSSIM on the same data, so quality numbers
are directly comparable across tools.

Author: JiajunHuang1999

README.md

@@ -122,6 +122,7 @@ fsz -t f64 -i data.d64 -eb rel 1e-3                   any mode on float64 data
 | `--bare` | with `-z`: write the raw bitstream without the 56-byte header. Decompressing a bare stream needs `-n` (or `-d`) and `-eb abs`. |
 | `-n N` | element count, for decompressing bare streams. |
 | `--csv` | additionally print a one-line machine-readable record. |
+| `--ssim` | additionally compute the windowed structural similarity of the reconstruction (slower; uses the `-d` shape). PSNR and NRMSE are always reported. |
 
 Exit codes: `0` success (error check passed), `1` bound violated, `2` usage or
 file errors.

tools/fsz.cu

@@ -18,6 +18,7 @@
 
 #include <cerrno>
 #include <cmath>
+#include <limits>
 #include <cstdarg>
 #include <cstdint>
 #include <cstdio>
@@ -122,6 +123,8 @@ void usage(std::FILE* out) {
         "  -n N                element count, for decompressing a bare bitstream\n"
         "  --bare              with -z, write the bitstream with no file header\n"
         "  --csv               also print one machine-readable csv line\n"
+        "  --ssim              also compute the windowed structural similarity\n"
+        "                      of the reconstruction (slower; uses -d shape)\n"
         "  -h, --help          this message\n"
         "\n"
         "Input is a raw float32 or float64 array of any dimensionality; up to\n"
@@ -299,6 +302,7 @@ struct Options {
     bool                  force_x    = false;
     bool                  bare       = false;
     bool                  csv        = false;
+    bool                  ssim       = false;
 };
 
 void parse_dim_list(int argc, char** argv, int& i, std::vector<uint64_t>& dims) {
@@ -358,6 +362,8 @@ Options parse_args(int argc, char** argv) {
             o.bare = true;
         } else if (std::strcmp(a, "--csv") == 0) {
             o.csv = true;
+        } else if (std::strcmp(a, "--ssim") == 0) {
+            o.ssim = true;
         } else if (std::strcmp(a, "-i") == 0) {
             if (i + 1 >= argc) die_usage("-i needs a path");
             o.in_path = argv[++i];
@@ -438,6 +444,66 @@ void resolve_shape(const Options& o, uint64_t n_file, const std::string& what,
     n = prod;
 }
 
+// Windowed structural similarity over the recorded shape: 7-point windows
+// with stride 2 per axis, the window-local range of the original data as the
+// luminance scale, and the mean over all windows. This matches the windowed
+// convention used by the common compression quality tools, so values are
+// directly comparable.
+template <class T>
+double windowed_ssim(const std::vector<T>& ori, const std::vector<T>& dec,
+                     const uint64_t dims[3], uint32_t ndims) {
+    size_t s2 = 1, s1 = 1, s0 = 1;
+    if (ndims >= 3)      { s2 = dims[0]; s1 = dims[1]; s0 = dims[2]; }
+    else if (ndims == 2) { s1 = dims[0]; s0 = dims[1]; }
+    else                 { s0 = dims[0]; }
+    const size_t w2 = (s2 > 1) ? 7 : 1;
+    const size_t w1 = (s1 > 1) ? 7 : 1;
+    const size_t w0 = 7;
+    if (w0 > s0 || w1 > s1 || w2 > s2)
+        die("--ssim needs every recorded dimension to hold at least 7 values");
+    double sum = 0.0;
+    size_t nw = 0;
+    const double np = (double)(w0 * w1 * w2);
+    for (size_t o2 = 0; o2 + w2 <= s2; o2 += 2) {
+        for (size_t o1 = 0; o1 + w1 <= s1; o1 += 2) {
+            for (size_t o0 = 0; o0 + w0 <= s0; o0 += 2) {
+                double xmin = (double)ori[o0 + s0 * (o1 + s1 * o2)];
+                double xmax = xmin, xs = 0.0, ys = 0.0;
+                for (size_t i2 = 0; i2 < w2; ++i2)
+                    for (size_t i1 = 0; i1 < w1; ++i1)
+                        for (size_t i0 = 0; i0 < w0; ++i0) {
+                            const size_t idx =
+                                (o0 + i0) + s0 * ((o1 + i1) + s1 * (o2 + i2));
+                            const double x = (double)ori[idx];
+                            const double y = (double)dec[idx];
+                            if (x < xmin) xmin = x;
+                            if (x > xmax) xmax = x;
+                            xs += x; ys += y;
+                        }
+                const double xm = xs / np, ym = ys / np;
+                double vx = 0.0, vy = 0.0, vxy = 0.0;
+                for (size_t i2 = 0; i2 < w2; ++i2)
+                    for (size_t i1 = 0; i1 < w1; ++i1)
+                        for (size_t i0 = 0; i0 < w0; ++i0) {
+                            const size_t idx =
+                                (o0 + i0) + s0 * ((o1 + i1) + s1 * (o2 + i2));
+                            const double dx = (double)ori[idx] - xm;
+                            const double dy = (double)dec[idx] - ym;
+                            vx += dx * dx; vy += dy * dy; vxy += dx * dy;
+                        }
+                vx /= np; vy /= np; vxy /= np;
+                const double L  = xmax - xmin;
+                const double c1 = (L == 0.0) ? 0.01 * 0.01 : 0.01 * 0.01 * L * L;
+                const double c2 = (L == 0.0) ? 0.03 * 0.03 : 0.03 * 0.03 * L * L;
+                sum += ((2.0 * xm * ym + c1) * (2.0 * vxy + c2))
+                     / ((xm * xm + ym * ym + c1) * (vx + vy + c2));
+                ++nw;
+            }
+        }
+    }
+    return sum / (double)nw;
+}
+
 // Compress path, shared by report mode and -z. Report mode differs only in that
 // it writes no output file. T is the element type of the input array.
 
@@ -520,13 +586,28 @@ int run_compress(const Options& o, bool write_output) {
 
     double max_err   = 0.0;
     size_t nonfinite = 0;
+    double sum_sq    = 0.0;
+    double omin      = (double)h_in[0];
+    double omax      = (double)h_in[0];
     for (size_t i = 0; i < n; ++i) {
+        const double x = (double)h_in[i];
+        if (x < omin) omin = x;
+        if (x > omax) omax = x;
         const T v = h_out[i];
         if (!std::isfinite(v)) { ++nonfinite; continue; }
-        const double e = std::fabs((double)h_in[i] - (double)v);
+        const double e = std::fabs(x - (double)v);
         if (e > max_err) max_err = e;
+        sum_sq += e * e;
     }
     const bool ok = (nonfinite == 0) && (max_err <= 1.01 * (double)eb_abs);
+    const double range = omax - omin;
+    const double rmse  = std::sqrt(sum_sq / (double)n);
+    const double psnr  = (range > 0.0 && rmse > 0.0)
+                       ? 20.0 * std::log10(range / rmse)
+                       : std::numeric_limits<double>::infinity();
+    const double nrmse = (range > 0.0) ? rmse / range : 0.0;
+    double ssim = -1.0;
+    if (o.ssim) ssim = windowed_ssim(h_in, h_out, dims, ndims);
 
     uint64_t total_bytes = 0;
     if (write_file) {
@@ -576,6 +657,10 @@ int run_compress(const Options& o, bool write_output) {
     }
     std::printf("  compress   = %.2f ms (%.2f GB/s)\n", c_ms, gbps(in_bytes, c_ms));
     std::printf("  decompress = %.2f ms (%.2f GB/s)\n", d_ms, gbps(in_bytes, d_ms));
+    std::printf("  PSNR       = %.2f dB\n", psnr);
+    std::printf("  NRMSE      = %.6e\n", nrmse);
+    if (o.ssim)
+        std::printf("  SSIM       = %.6f\n", ssim);
     std::printf("  max_err    = %.6e  (bound %.6e)\n", max_err, (double)eb_abs);
     if (ok) {
         std::printf("  %sPASS%s\n", green(), reset());