Test Images, Datasets & Quality Metrics Tooling

"Forty-seven years of being blurred, sharpened, compressed, and denoised, and they finally let me retire. My one request: let the astronaut take it from here. She is used to noise."

A Graciously Retired Test Image

The previous section (Section 8.2) made pipelines fast; this one makes their evaluation trustworthy. Every experiment in Part I needed three ingredients we mostly took for granted: an input image, a degraded-versus-clean pair, and a number that says how well a method did. Here we treat each ingredient as infrastructure: where the standard images live, which datasets benchmark the restoration tasks of Chapter 7, and how to compute the quality metrics introduced in Chapter 1 so that a reviewer, or your own CI pipeline, gets the same number you did.

1. The Standard Test Images Beginner

A handful of photographs have been blurred and denoised millions of times because everyone agreed to use the same ones. The historical home is the USC-SIPI database (maintained since 1977), source of the mandrill, the peppers, and the house; the Kodak suite adds 24 film scans that remain the default smoke test for compression and denoising. For code, the most convenient gallery ships inside scikit-image itself: skimage.data bundles the classic cameraman, an astronaut portrait, coins on a textured background, text, and more, each one load-by-function-call. Code 8.3.1 inventories it.

# Inventory the built-in test gallery: load each classic image by
# function call and report its shape and dtype, surfacing the mixed
# uint8/float64 conventions that Section 8.1 warned about.
from skimage import data

for name in ["camera", "astronaut", "coins", "text", "checkerboard",
             "page", "moon", "shepp_logan_phantom"]:
    img = getattr(data, name)()          # each image is a function call
    print(f"{name:20s} {img.shape!s:16s} {img.dtype}")

camera               (512, 512)       uint8
astronaut            (512, 512, 3)    uint8
coins                (303, 384)       uint8
text                 (172, 448)       uint8
checkerboard         (200, 200)       uint8
page                 (191, 384)       uint8
moon                 (512, 512)       uint8
shepp_logan_phantom  (400, 400)       float64

Synthetic targets complement the photographs because their ground truth is mathematical. The Shepp-Logan phantom (printed by Code 8.3.1) has known ellipse boundaries, which is why CT reconstruction papers have used it since 1974. Checkerboards and ramps expose interpolation and quantization artifacts from Chapter 1. The most instructive of all is the zone plate, $I(r) = \tfrac{1}{2}\left(1 + \cos(\alpha r^2)\right)$, whose local frequency grows linearly with radius (the rings pack tighter the farther out you go, because differentiating the $r^2$ phase gives a frequency proportional to $r$): a single image that sweeps every spatial frequency, so aliasing from poor resampling (Chapter 5) and the band behavior of filters (Chapter 4) appear as visible rings exactly where theory predicts. Code 8.3.2 builds one, then manufactures the (clean, noisy) pair every denoising experiment needs.

# Build a zone plate (a radial chirp with known frequency content)
# and add calibrated Gaussian noise to it, manufacturing a synthetic
# (clean, noisy) restoration pair with exact mathematical ground truth.
import numpy as np

def zone_plate(n=512, alpha=0.4):
    """Radial chirp: local frequency grows linearly with radius."""
    y, x = np.mgrid[-n // 2:n // 2, -n // 2:n // 2].astype(np.float64)
    r2 = (x**2 + y**2) / n                # normalized squared radius
    return 0.5 * (1.0 + np.cos(alpha * r2))

clean = zone_plate()                      # float64 in [0, 1], exact ground truth
rng = np.random.default_rng(7)
noisy = np.clip(clean + rng.normal(0.0, 25 / 255, clean.shape), 0.0, 1.0)

2. Benchmark Datasets for Restoration Intermediate

When a method must be compared against published results, single images give way to standard datasets, and the choice is rarely innocent: a denoiser that wins on synthetic Gaussian noise can lose badly on real sensor noise, so the dataset you pick quietly decides which method looks best. Table 8.3.1 lists the ones that dominate the literature for the tasks of Chapter 7; the names recur in virtually every paper you will read after this part of the book.

Two hygiene rules attach to every row. First, respect the published splits: DIV2K's validation images, for example, appear inside other datasets' test sets, and accidental train-test contamination has invalidated more than one results table. Second, record the exact degradation recipe (noise sigma and dtype scale, downsampling kernel for super-resolution, JPEG quality) alongside your scores; "PSNR 31.2 on BSD68" is meaningless without "sigma = 25, grayscale, values in [0, 255]". A quick histogram audit of any downloaded dataset, with the tools from Chapter 2, catches the most common surprises: 16-bit files masquerading as 8-bit, and already-compressed "clean" references.

3. Metric Tooling: PSNR, SSIM & the Default-Parameter Trap Intermediate

Chapter 1 defined the two workhorse full-reference metrics: peak signal-to-noise ratio, $$\mathrm{PSNR} = 10 \log_{10} \frac{L^2}{\mathrm{MSE}},$$ with $L$ the dynamic range, and the structural similarity index $$\mathrm{SSIM}(x, y) = \frac{(2\mu_x \mu_y + C_1)(2\sigma_{xy} + C_2)}{(\mu_x^2 + \mu_y^2 + C_1)(\sigma_x^2 + \sigma_y^2 + C_2)},$$ computed over local windows and averaged. The tooling lives in skimage.metrics, and it hides a reproducibility trap: the original SSIM publication specifies an 11×11 Gaussian-weighted window with $\sigma = 1.5$ and population statistics, while scikit-image defaults to a 7×7 uniform window with sample statistics. (Population statistics divide a variance by $N$; sample statistics divide by $N-1$, a small bias correction that shifts the score slightly.) Both are legitimate SSIMs; they are just different numbers, as the illustration below dramatizes with two judges scoring the same painting through different windows. Code 8.3.3 computes both and makes the second trap, data_range, explicit.

# Score the same image pair two ways: PSNR plus SSIM under both the
# scikit-image defaults and the original Wang et al. (2004) protocol,
# making the data_range and window-choice traps explicit.
from skimage.metrics import peak_signal_noise_ratio, structural_similarity

# Floats demand an explicit range; guessing it is the classic silent bug.
psnr = peak_signal_noise_ratio(clean, noisy, data_range=1.0)

ssim_default = structural_similarity(clean, noisy, data_range=1.0)

ssim_paper = structural_similarity(      # the Wang et al. (2004) protocol
    clean, noisy, data_range=1.0,
    gaussian_weights=True, sigma=1.5,    # 11x11 Gaussian window
    use_sample_covariance=False)         # population (1/N) statistics

print(f"PSNR              : {psnr:.2f} dB")
print(f"SSIM (defaults)   : {ssim_default:.4f}")
print(f"SSIM (paper flags): {ssim_paper:.4f}")

PSNR              : 20.17 dB
SSIM (defaults)   : 0.4664
SSIM (paper flags): 0.3923

4. Beyond Pixels: Perceptual and No-Reference Metrics Advanced

PSNR and SSIM compare pixels and local statistics, and they famously prefer over-smoothed results that humans dislike. LPIPS (Learned Perceptual Image Patch Similarity; Zhang et al., 2018) scores the distance between deep network features of the two images instead, and correlates far better with human judgments on restoration outputs. When no clean reference exists at all (real-world enhancement, in-the-wild quality monitoring), no-reference metrics such as NIQE, BRISQUE, and the learned MUSIQ and CLIP-IQA estimate quality from the image alone. Figure 8.3.1 organizes the families into the decision you actually face.

In practice you rarely implement any of these. The pyiqa toolbox wraps more than forty metrics, classical and learned, behind one interface, as Code 8.3.4 shows.

# Score one image pair with four metrics through a single interface:
# pyiqa exposes each metric's direction (lower_better) and reference
# requirement (metric_mode), so the harness picks the call form safely.
import pyiqa, torch

device = "cuda" if torch.cuda.is_available() else "cpu"
metrics = {name: pyiqa.create_metric(name, device=device)
           for name in ["psnr", "ssim", "lpips", "niqe"]}

# Tensors in NCHW, float, [0, 1]; pyiqa also accepts file paths directly.
ref = torch.from_numpy(clean).float()[None, None]
deg = torch.from_numpy(noisy).float()[None, None]

for name, m in metrics.items():
    score = m(deg) if m.metric_mode == "NR" else m(deg, ref)
    print(f"{name:6s}: {float(score):8.4f}   (lower is better: {m.lower_better})")

The three threads of this chapter, choosing a library deliberately (Section 8.1), measuring and comparing speed honestly (Section 8.2), and scoring quality reproducibly (this section), only become a single skill when you wire them into one harness. Put them together in the hands-on lab below, which turns this chapter into a runnable denoising benchmark studio you can keep and reuse.

pip install opencv-python scikit-image scipy numpy

import numpy as np
from skimage import data, img_as_float

def make_pairs(sigma=25 / 255, seed=7):
    rng = np.random.default_rng(seed)
    images = {
        "camera": data.camera(),
        "coins": data.coins(),
        "moon": data.moon(),
    }
    pairs = {}
    for name, raw in images.items():
        clean = img_as_float(raw)              # uint8 -> float64 in [0, 1]
        # TODO: add Gaussian noise of standard deviation `sigma`, then clip to [0, 1]
        # Hint: rng.normal(0.0, sigma, clean.shape), then np.clip(...)
        noisy = ...
        pairs[name] = (clean, noisy)
    return pairs

pairs = make_pairs()
print({k: v[0].shape for k, v in pairs.items()})

import cv2
from scipy import ndimage as ndi
from skimage.restoration import denoise_bilateral

def denoise_opencv_nlm(img):
    u8 = (np.clip(img, 0, 1) * 255).astype(np.uint8)   # OpenCV wants uint8
    out = cv2.fastNlMeansDenoising(u8, h=12)
    return out.astype(np.float64) / 255.0

def denoise_scipy_gaussian(img):
    # TODO: return a Gaussian-blurred copy using scipy.ndimage (sigma about 1.0)
    # Hint: ndi.gaussian_filter keeps float input float; no conversion needed
    return ...

def denoise_skimage_bilateral(img):
    # skimage works natively in float [0, 1]
    return denoise_bilateral(img, sigma_color=0.1, sigma_spatial=3)

denoisers = {
    "opencv_nlm": denoise_opencv_nlm,
    "scipy_gaussian": denoise_scipy_gaussian,
    "skimage_bilateral": denoise_skimage_bilateral,
}

import time

def timed(fn, img, repeats=3):
    fn(img)                                  # warm-up: ignore the first call
    best = float("inf")
    for _ in range(repeats):
        t0 = time.perf_counter()
        out = fn(img)
        best = min(best, time.perf_counter() - t0)
    # TODO: return both the denoised image `out` and the best time in milliseconds
    return ...

from skimage.metrics import peak_signal_noise_ratio, structural_similarity

def score(clean, restored):
    psnr = peak_signal_noise_ratio(clean, restored, data_range=1.0)
    ssim_default = structural_similarity(clean, restored, data_range=1.0)
    # TODO: compute SSIM under the Wang et al. (2004) protocol
    # Hint: gaussian_weights=True, sigma=1.5, use_sample_covariance=False
    ssim_paper = ...
    return psnr, ssim_default, ssim_paper

rows = []
for img_name, (clean, noisy) in pairs.items():
    for method_name, fn in denoisers.items():
        restored, ms = timed(fn, noisy)
        psnr, ssim_d, ssim_p = score(clean, restored)
        rows.append({
            "image": img_name, "method": method_name,
            "psnr": psnr, "ssim_default": ssim_d,
            "ssim_paper": ssim_p, "time_ms": ms,
        })
print(f"collected {len(rows)} runs")           # expect 3 images x 3 methods = 9

import csv, statistics, skimage, cv2 as _cv

# Leaderboard: mean paper-SSIM and mean time per method
by_method = {}
for r in rows:
    by_method.setdefault(r["method"], []).append(r)
print(f"{'method':20s} {'SSIM(paper)':>12s} {'time_ms':>10s}")
for method, recs in sorted(
        by_method.items(),
        key=lambda kv: statistics.mean(r['ssim_paper'] for r in kv[1]),
        reverse=True):
    s = statistics.mean(r["ssim_paper"] for r in recs)
    t = statistics.mean(r["time_ms"] for r in recs)
    print(f"{method:20s} {s:12.4f} {t:10.1f}")

# TODO: write all `rows` to results.csv with csv.DictWriter, prefixing
# a comment line that records skimage.__version__ and _cv.__version__

method               SSIM(paper)    time_ms
skimage_bilateral         0.8021       38.4
opencv_nlm                0.7864       21.7
scipy_gaussian            0.6190        2.1

import csv, statistics, time
import numpy as np
import cv2
from scipy import ndimage as ndi
import skimage
from skimage import data, img_as_float
from skimage.restoration import denoise_bilateral
from skimage.metrics import peak_signal_noise_ratio, structural_similarity

# Step 1: reproducible (clean, noisy) pairs
def make_pairs(sigma=25 / 255, seed=7):
    rng = np.random.default_rng(seed)
    images = {"camera": data.camera(), "coins": data.coins(), "moon": data.moon()}
    pairs = {}
    for name, raw in images.items():
        clean = img_as_float(raw)
        noisy = np.clip(clean + rng.normal(0.0, sigma, clean.shape), 0.0, 1.0)
        pairs[name] = (clean, noisy)
    return pairs

# Step 2: three denoisers, one signature (float [0, 1] in and out)
def denoise_opencv_nlm(img):
    u8 = (np.clip(img, 0, 1) * 255).astype(np.uint8)
    return cv2.fastNlMeansDenoising(u8, h=12).astype(np.float64) / 255.0

def denoise_scipy_gaussian(img):
    return ndi.gaussian_filter(img, sigma=1.0)

def denoise_skimage_bilateral(img):
    return denoise_bilateral(img, sigma_color=0.1, sigma_spatial=3)

denoisers = {
    "opencv_nlm": denoise_opencv_nlm,
    "scipy_gaussian": denoise_scipy_gaussian,
    "skimage_bilateral": denoise_skimage_bilateral,
}

# Step 3: honest timing (warm-up, then best of N)
def timed(fn, img, repeats=3):
    fn(img)
    best = float("inf")
    for _ in range(repeats):
        t0 = time.perf_counter()
        out = fn(img)
        best = min(best, time.perf_counter() - t0)
    return out, best * 1000.0

# Step 4: PSNR plus two SSIM protocols
def score(clean, restored):
    psnr = peak_signal_noise_ratio(clean, restored, data_range=1.0)
    ssim_default = structural_similarity(clean, restored, data_range=1.0)
    ssim_paper = structural_similarity(
        clean, restored, data_range=1.0,
        gaussian_weights=True, sigma=1.5, use_sample_covariance=False)
    return psnr, ssim_default, ssim_paper

# Step 5: run the grid
pairs = make_pairs()
rows = []
for img_name, (clean, noisy) in pairs.items():
    for method_name, fn in denoisers.items():
        restored, ms = timed(fn, noisy)
        psnr, ssim_d, ssim_p = score(clean, restored)
        rows.append({
            "image": img_name, "method": method_name,
            "psnr": round(psnr, 4), "ssim_default": round(ssim_d, 4),
            "ssim_paper": round(ssim_p, 4), "time_ms": round(ms, 1),
        })

# Step 6: leaderboard + reproducible CSV
by_method = {}
for r in rows:
    by_method.setdefault(r["method"], []).append(r)
print(f"{'method':20s} {'SSIM(paper)':>12s} {'time_ms':>10s}")
for method, recs in sorted(
        by_method.items(),
        key=lambda kv: statistics.mean(r['ssim_paper'] for r in kv[1]),
        reverse=True):
    s = statistics.mean(r["ssim_paper"] for r in recs)
    t = statistics.mean(r["time_ms"] for r in recs)
    print(f"{method:20s} {s:12.4f} {t:10.1f}")

with open("results.csv", "w", newline="") as f:
    f.write(f"# skimage={skimage.__version__} opencv={cv2.__version__}\n")
    writer = csv.DictWriter(f, fieldnames=list(rows[0].keys()))
    writer.writeheader()
    writer.writerows(rows)
print("wrote results.csv")