Structural Similarity Index

The structural similarity index (SSIM) compares a test image with a reference image of the same size, window by window, through three local statistics: the mean intensity, the contrast (standard deviation) and the correlation between the two patches. Each window gets a score of at most 1, reached only when the patches are identical; the scores form an SSIM map, and their mean is the single number reported in papers. Wang, Bovik, Sheikh and Simoncelli (2004) proposed it as an alternative to error measures such as the mean squared error and PSNR, on the premise that viewers judge images by their structure, and it is now reported next to PSNR in compression, restoration, super-resolution and view synthesis. Its value also depends on implementation choices that must be stated for two numbers to be comparable: the window, whether the image is downsampled first, the colour channels and the assumed dynamic range.

Formula

For two aligned, non-negative image patches $x$ and $y$ (in quality assessment, $x$ is the reference), the 2004 paper defines a luminance, a contrast and a structure comparison (eqs. 6, 9 and 10):

$$l(x, y) = \frac{2\mu_x\mu_y + C_1}{\mu_x^2 + \mu_y^2 + C_1}, \qquad c(x, y) = \frac{2\sigma_x\sigma_y + C_2}{\sigma_x^2 + \sigma_y^2 + C_2}, \qquad s(x, y) = \frac{\sigma_{xy} + C_3}{\sigma_x\sigma_y + C_3},$$

where $\mu$ is the mean, $\sigma$ the standard deviation and $\sigma_{xy}$ the covariance of the patches. Without the constants, $s$ is the correlation coefficient, the cosine of the angle between the mean-removed patches. The general index is $[l]^{\alpha}[c]^{\beta}[s]^{\gamma}$ (eq. 12); the paper sets $\alpha = \beta = \gamma = 1$ and $C_3 = C_2/2$, which cancels the $\sigma_x\sigma_y$ terms and gives the form everyone computes (eq. 13):

$$\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)} .$$

Constants. $C_1 = (K_1 L)^2$ and $C_2 = (K_2 L)^2$, with $L$ the dynamic range of the pixel values (255 for 8-bit grayscale) and $K_1 = 0.01$, $K_2 = 0.03$ (eq. 7 and §III-C). The authors call these values “somewhat arbitrary” and report that performance was fairly insensitive to them. For 8-bit data, $C_1 = 6.5025$ and $C_2 = 58.5225$; for data in $[0, 1]$ with $L = 1$, $C_1 = 10^{-4}$ and $C_2 = 9 \times 10^{-4}$, which gives the same SSIM, since all terms scale with $L^2$. The constants keep the ratios stable where $\mu_x^2 + \mu_y^2$ or $\sigma_x^2 + \sigma_y^2$ is near zero; without them the index reduces to the earlier universal quality index (UQI), which the authors found unstable in nearly flat regions.

Local statistics. The statistics are computed in an $11 \times 11$ circularly symmetric Gaussian window $w$ with a standard deviation of 1.5 pixels, normalized to sum to one (§III-C, eqs. 14–16):

$$\mu_x = \sum_i w_i x_i, \qquad \sigma_x^2 = \sum_i w_i (x_i - \mu_x)^2, \qquad \sigma_{xy} = \sum_i w_i (x_i - \mu_x)(y_i - \mu_y).$$

The window replaced the $8 \times 8$ square window of UQI, whose maps showed blocking artifacts; it is the sampled and truncated kernel of the Gaussian filter. The weighted variance has no $N-1$ correction (the paper’s unweighted eq. 4 has one).

Pooling. The window slides pixel by pixel; its score at each position forms the SSIM map, and the mean SSIM over the $M$ window positions is the image score (eq. 17):

$$\mathrm{MSSIM}(X, Y) = \frac{1}{M}\sum_{j=1}^{M} \mathrm{SSIM}(x_j, y_j).$$

The paper calls this MSSIM; it is what “SSIM” means in results tables. In the authors’ MATLAB code the windows are the “valid” positions, those that fit entirely inside the image, so an $H \times W$ image gives an $(H-10) \times (W-10)$ map.

MS-SSIM

Wang, Simoncelli and Bovik (2003) argued that the right scale for single-scale SSIM depends on viewing conditions, display resolution and distance, and proposed a multi-scale version. The images are low-pass filtered and downsampled by 2 repeatedly, as in an image pyramid; scale 1 is the original and scale $M$ the coarsest. Contrast and structure are compared at every scale and luminance only at the coarsest (eq. 7):

$$\mathrm{MS\text{-}SSIM}(x, y) = [l_M(x, y)]^{\alpha_M} \prod_{j=1}^{M} [c_j(x, y)]^{\beta_j}\,[s_j(x, y)]^{\gamma_j}.$$

With $\alpha_j = \beta_j = \gamma_j$ and the exponents normalized to sum to one, the authors calibrated five scales with a perceptual experiment: eight subjects chose, for ten $64 \times 64$ originals, synthesized images of equal perceived quality across scales. The resulting exponents are $\beta_1 = \gamma_1 = 0.0448$, $\beta_2 = \gamma_2 = 0.2856$, $\beta_3 = \gamma_3 = 0.3001$, $\beta_4 = \gamma_4 = 0.2363$ and $\alpha_5 = \beta_5 = \gamma_5 = 0.1333$ (§3.2), so the finest scale carries the least weight. The paper does not specify the low-pass filter; the authors’ msssim.m uses a $2 \times 2$ average, computes the mean of the combined contrast–structure map $(2\sigma_{xy} + C_2)/(\sigma_x^2 + \sigma_y^2 + C_2)$ at each scale and raises those means, not the per-window values, to the exponents.

Intuition

Three questions per window. Is the brightness the same ($l$), is the amount of variation the same ($c$), and do the two patches vary together ($s$)? The authors motivate the split by image formation: illumination changes the local mean and contrast but not the structure of the objects, so they define structure as what remains after both are removed (§III-B). Two perceptual properties are built in: with $C_1$ small, $l$ depends only on the relative luminance change, qualitatively consistent with Weber’s law (eq. 8), and $c$ is less sensitive to a given contrast change when the base contrast is high, consistent with contrast masking.

A locally normalized MSE. With $C_3 = C_2/2$, SSIM is the product of a luminance term and a combined contrast–structure term, and each of them is one minus a normalized squared error (Brunet, Vrscay & Wang, 2012, eqs. 4–9):

$$1 - l = \frac{(\mu_x - \mu_y)^2}{\mu_x^2 + \mu_y^2 + C_1}, \qquad 1 - cs = \frac{\sigma_x^2 + \sigma_y^2 - 2\sigma_{xy}}{\sigma_x^2 + \sigma_y^2 + C_2}.$$

The numerators add up to the local MSE in the window (a one-line identity: the squared error of the means plus the squared error of the mean-removed patches). SSIM therefore measures the same errors as the MSE, but divides them by the local signal energy, separately for the mean and for the variation. The same error costs little where the reference is bright or busy and a lot where it is dark or flat. Nilsson and Akenine-Möller (2020, §2) summarize work by Dosselmann and Yang (a 2008 technical report and the 2011 journal paper) showing a direct mathematical transform between SSIM with zero constants and a local MSE computed over the same window, and read such results as questioning the claim that SSIM is perception-based.

Interpretation

Range. For non-negative images, SSIM lies in $[-1, 1]$, and it equals 1 if and only if the patches are identical (Brunet et al., 2012, §II); the 2004 paper lists symmetry, boundedness by 1 and this unique maximum as design conditions (§III-B). Negative values come from the structure term, when the patches are anticorrelated. The mean over windows inherits the range, and it is 1 only for identical images.

Values. Higher is better, and the scale is far from linear in perceived quality: in the 2004 paper’s scatter plot (Fig. 8, our reading of the figure), compressed images with MSSIM below about 0.8 all have mean opinion scores below about 40 on a 100-point scale, and the rest of the subjective range is spread over 0.8 to 1. As with PSNR, there is no universal threshold for good; scores are comparable only between methods evaluated on the same images under the same convention.

Properties. SSIM is symmetric and differentiable. It is not a distance metric: $1 - \mathrm{SSIM}$ violates the triangle inequality. Brunet et al. (2012) proved that the square roots $\sqrt{1 - l}$ and $\sqrt{1 - cs}$, computed for a pair of signals, are normalized metrics (on the means and on the mean-removed signals respectively), and combined them into metrics on the whole signal (Theorems 2.6 and 2.9); these are modified indices, not SSIM itself.

Agreement with human scores. On the authors’ database of 344 JPEG and JPEG2000 images from 29 originals (the LIVE JPEG/JPEG2000 database, as the MS-SSIM paper calls it), the 2004 paper reports a Spearman rank-order correlation (SROCC) with mean opinion scores of 0.963 for MSSIM against 0.901 for PSNR, and a Pearson correlation after a logistic fit (CC) of 0.967 against 0.905 (Table I); it computed SSIM on the luminance, averaged over $2 \times 2$ blocks and downsampled by 2 (§IV-B). On the same database, MS-SSIM scored 0.969 CC and 0.966 SROCC, while single-scale SSIM varied with scale between 0.936 SROCC at the coarsest of five scales and 0.959 one scale down, with 0.945 at full resolution (Wang et al., 2003, Table 1). Databases with more distortion types give lower correlations: on the authors’ project page, ssim.m (with downsampling and grayscale conversion by MATLAB’s rgb2gray) reaches SROCC 0.9479 on LIVE (the table does not name the release; elsewhere the page describes the full database of 779 distorted images), 0.8756 on CSIQ and 0.7749 on TID2008.

When to Use It

SSIM is a full-reference index: it needs a pixel-aligned reference, like PSNR, and is usually reported with it.

• Compression, where the 2004 paper validated it and suggested it for rate–distortion optimization.
• Restoration: image denoising, deblurring and inpainting. Super-resolution benchmarks compute it on the luma channel with the border cropped, the conventions described in PSNR.
• View synthesis, where renderings of held-out poses are scored by PSNR, SSIM and LPIPS (see neural rendering and novel view synthesis).
• Training losses, as $1 - \mathrm{SSIM}$, discussed under Variants.

It fits when the errors to be judged are local loss of detail, noise and compression artifacts, and the reference is the image the output should match. It is a poor choice for misaligned outputs, for colour fidelity, for generated content without a unique correct answer and, used alone, for comparing methods whose differences are small relative to the convention effects shown below.

Failure Cases

Implementation conventions change the number. The same pair of images gives different SSIMs under different, equally common, conventions, and the example below shows differences far larger than, for instance, the 0.013 SSIM gap between NeRF and LLFF on real forward-facing scenes in NeRF’s Table 1 (see neural rendering). The main choices:

• Window. Gaussian $11 \times 11$, $\sigma = 1.5$, as in the paper, or a square uniform window. scikit-image 0.25.2 defaults to gaussian_weights=False and win_size=7, a $7 \times 7$ uniform window, with use_sample_covariance=True, the $N/(N-1)$ correction; its documentation states that matching Wang et al. requires gaussian_weights=True, sigma=1.5, use_sample_covariance=False and an explicit data_range.
• Downsampling before computing. The authors’ project page recommends, for images viewed at a typical distance (3 to 5 image heights), averaging $F \times F$ blocks and downsampling by $F = \max(1, \operatorname{round}(N/256))$, with $N$ the image height or width, before applying the single-scale code ssim_index.m. Their ssim.m does this automatically, using the smaller of the two dimensions; ssim_index.m and scikit-image’s function do not. A $512 \times 512$ image is therefore scored at $256 \times 256$ by ssim.m and at full resolution by scikit-image.
• Colour. The paper and the project page compute SSIM on one channel: luminance in the 2004 experiments, MATLAB’s rgb2gray for the page’s database results. scikit-image, given channel_axis, instead averages the per-channel SSIMs of an RGB image, and Zhao et al. (2017) report SSIM-type indices averaged over the three colour channels. Luma in the BT.601 studio range (16 to 235) with $L = 255$ is a third variant.
• Data range. $L$ must be the range of the data format. Passing $L = 255$ for images in $[0, 1]$ makes $C_1$ and $C_2$ dominate every window and drives SSIM to 1. In scikit-image 0.25.2, floating-point inputs without data_range raise an error, although the Notes section of the function’s documentation still describes an estimate of 2 for floats.
• Borders. The authors’ code averages only valid windows; implementations that pad the image and keep every pixel include windows that see the padding. scikit-image pads by reflection and then crops the window radius from the map, which matches the valid convention.

Insensitivity to some changes. A uniform brightness shift or a contrast stretch keeps the structure term near 1 and scores high, by design: the 2004 paper’s motivating figure (Fig. 2) shows images of a boat with MSE 210 where the mean shift scores 0.9900 and the contrast stretch 0.9168, against 0.6949 for JPEG compression. That suits the authors’ argument that the structural information is recoverable, but a viewer comparing two renderings may well see the shift. Zhao et al. (2017) found that networks trained with SSIM-based losses tended to shift brightness and dull colours for the same reason. Converting to grayscale first hides any colour change that preserves luma, and Nilsson and Akenine-Möller (2020, §4.2) show visibly different colours scoring above 0.99.

Dark and flat regions. Because the errors are normalized by $\mu_x^2 + \mu_y^2$ and $\sigma_x^2 + \sigma_y^2$, a small error in a dark or flat region can produce a very low score. Nilsson and Akenine-Möller (2020, Fig. 4) report MSSIM 0.61914 between a black image and one of grey level 2/255, which are hard to tell apart, but 0.99047 between a white image and one of level 222/255, which are clearly different. Pixel-scale checkerboards that are invisible at normal viewing can score near the component minima they derive, for example $c = 0.0036$ against a uniform grey (§3.1). In the other direction, SSIM maps stay high near strong edges even when noise of the same magnitude covers both sides (their Fig. 8); in our reading, the edge’s large variance masks the error.

Negative values and non-integer exponents. The structure term can be negative, and a negative number raised to a non-integer power, as in MS-SSIM’s exponents, is undefined in real arithmetic; Nilsson and Akenine-Möller (2020, §3.1) point out that implementations can then return NaN. Their mirrored-gradient example scores −0.82 at $16 \times 16$ pixels (Fig. 6).

Misalignment. Like the MSE, SSIM compares pixels at the same positions; a one-pixel shift wrecks the score in fine texture while leaving flat regions at 1. In the two-alternative perceptual tests of the LPIPS authors, SSIM agreed with human choices no more often than the $\ell_2$ distance (figures in novel view synthesis).

Variants

• MS-SSIM, the multi-scale form above, is common for compression and as a loss.
• Dissimilarities. $1 - \mathrm{SSIM}$ serves as a loss; Monodepth uses $(1 - \mathrm{SSIM})/2$ with a $3 \times 3$ box window, combined with an $L^1$ term (see monocular depth estimation), and 3D Gaussian splatting mixes $L^1$ with a D-SSIM term.
• SSIM as a training loss. Zhao et al. (2017) trained small convolutional networks for joint demosaicking and denoising with $\ell_2$, $\ell_1$, SSIM, MS-SSIM and a mixture $\alpha \cdot \mathcal{L}_{\text{MS-SSIM}} + (1 - \alpha) \cdot G_{\sigma_G^M} \cdot \mathcal{L}_{\ell_1}$ with $\alpha = 0.84$, the $\ell_1$ term weighted by the Gaussian of the coarsest scale (eq. 14). For the loss they set all MS-SSIM exponents to 1 and replaced the pyramid by Gaussians of $\sigma_G = 0.5, 1, 2, 4, 8$ on the full-resolution patch (§III-C). The SSIM and MS-SSIM networks did not clearly beat $\ell_1$, which the authors attribute to colour and brightness shifts; in their Table I the MS-SSIM-trained network even scores slightly lower MS-SSIM than the $\ell_1$-trained one (0.9741 against 0.9745), while the mixture scores best on every quality index (0.9757 MS-SSIM), and it was also best in their super-resolution and JPEG artifact removal tests. These results come from the arXiv v3 copy of the journal paper. Training on loss functions built from SSIM makes it harder to use SSIM as an independent evaluation of the same model.
• UQI, the special case $C_1 = C_2 = 0$, computed in $8 \times 8$ square windows (2004 paper, §III-B and §III-C).
• Weighted pooling. The 2004 paper notes that the map can be averaged with weights, for example from a region of interest or a model of fixations, but uses uniform weights.

Example

The example builds a synthetic $512 \times 512$ RGB image with 8-bit values in $[0, 255]$: colour ramps, two disks, a blue square, fine grey stripes, a band of smooth random texture and a near-black patch. It implements SSIM in NumPy and SciPy with the paper’s $11 \times 11$ Gaussian window, $\sigma = 1.5$, $L = 255$ and valid windows, on luma computed with MATLAB’s rgb2gray weights, and MS-SSIM as in the authors’ msssim.m. It then cross-checks against scikit-image 0.25.2 with the paper-matching flags; scores three distortions (Gaussian noise of $\sigma = 10$ per channel, a Gaussian blur and JPEG-like block DCT quantization with twice the ITU-T T.81 luminance table) under eight conventions; shows the data-range pitfall; and scores seven grayscale distortions tuned to the same MSE, with the mean of the SSIM map inside a flat, a textured and a dark region.

import numpy as np
from scipy.fft import dctn, idctn
from scipy.ndimage import correlate, gaussian_filter
from scipy.optimize import brentq
from scipy.signal import fftconvolve
from skimage.metrics import structural_similarity  # scikit-image 0.25.2, cross-check only

def window(kind="gauss", size=11, sigma=1.5):
    r = np.arange(size) - size // 2
    w = np.exp(-(r[:, None] ** 2 + r[None, :] ** 2) / (2 * sigma ** 2)) if kind == "gauss" \
        else np.ones((size, size))
    return w / w.sum()

def ssim(x, y, L=255.0, K1=0.01, K2=0.03, win=None, sample_cov=False, downsample=False,
         valid=True, full=False):
    """Mean SSIM of Wang et al. (2004), eq. 13 with weighted statistics (eqs. 14-16), over
    'valid' windows (valid=False: reflect-pad, keep all pixels). downsample=True applies
    the rule of the authors' ssim.m (the same result for F = 2 and even image sizes)."""
    x, y = np.asarray(x, float), np.asarray(y, float)
    if downsample:  # F = max(1, round(min(M, N) / 256)); F x F mean, keep every F-th pixel
        F = max(1, round(min(x.shape) / 256))
        x, y = [a if F == 1 else
                a[:a.shape[0] // F * F, :a.shape[1] // F * F]
                .reshape(a.shape[0] // F, F, a.shape[1] // F, F).mean(axis=(1, 3)) for a in (x, y)]
    w = window() if win is None else win
    f = (lambda a: fftconvolve(a, w[::-1, ::-1], mode="valid")) if valid else \
        (lambda a: correlate(a, w, mode="reflect"))
    mx, my = f(x), f(y)
    vx, vy, cxy = f(x * x) - mx ** 2, f(y * y) - my ** 2, f(x * y) - mx * my
    if sample_cov:  # N/(N-1) correction for an unweighted window of N pixels
        n = w.size
        vx, vy, cxy = vx * n / (n - 1), vy * n / (n - 1), cxy * n / (n - 1)
    C1, C2 = (K1 * L) ** 2, (K2 * L) ** 2
    m = (2 * mx * my + C1) * (2 * cxy + C2) / ((mx ** 2 + my ** 2 + C1) * (vx + vy + C2))
    return (m.mean(), m) if full else m.mean()

def ms_ssim(x, y, L=255.0, K1=0.01, K2=0.03,
            weights=(0.0448, 0.2856, 0.3001, 0.2363, 0.1333)):
    """MS-SSIM as in the authors' msssim.m: per scale the mean of c*s, luminance at the
    coarsest scale only, 2x2 averaging and downsampling by 2 between scales."""
    x, y, out = np.asarray(x, float), np.asarray(y, float), 1.0
    C1, C2 = (K1 * L) ** 2, (K2 * L) ** 2
    for j, wj in enumerate(weights):
        f = lambda a: fftconvolve(a, window(), mode="valid")
        mx, my = f(x), f(y)
        vx, vy, cxy = f(x * x) - mx ** 2, f(y * y) - my ** 2, f(x * y) - mx * my
        cs = (2 * cxy + C2) / (vx + vy + C2)
        if j < len(weights) - 1:
            out *= cs.mean() ** wj
            x, y = [a.reshape(a.shape[0] // 2, 2, a.shape[1] // 2, 2).mean(axis=(1, 3))
                    for a in (x, y)]
        else:
            out *= (((2 * mx * my + C1) / (mx ** 2 + my ** 2 + C1)) * cs).mean() ** wj
    return out

# A 512x512 RGB test image in [0, 255]: colour ramps, disks, a square, stripes, a smooth
# random texture, and a near-black patch.
rng = np.random.default_rng(0)
yy, xx = np.mgrid[0:512, 0:512].astype(float)
img = np.stack([40 + 160 * xx / 511, 50 + 140 * yy / 511, np.full_like(xx, 120)], -1)
img[(yy - 140) ** 2 + (xx - 150) ** 2 < 90 ** 2] = [210, 60, 50]
img[(yy - 380) ** 2 + (xx - 120) ** 2 < 60 ** 2] = [60, 170, 90]
img[300:460, 260:420] = [40, 70, 190]
img[40:220, 300:480] = 128 + 70 * np.sin(2 * np.pi * xx[40:220, 300:480] / 5)[..., None]
img[260:300, 20:500] = (128 + 60 * gaussian_filter(rng.standard_normal((40, 480)), 1.2)
                        / 0.25)[..., None].clip(0, 255)
img[470:510, 300:500] = 8                                   # near-black patch
ref = np.round(img).clip(0, 255)                           # 8-bit values, as floats
Y = lambda rgb: rgb @ np.array([0.2989, 0.5870, 0.1140])     # MATLAB rgb2gray (BT.601 weights)
Ys = lambda rgb: 16 + 219 * (rgb / 255) @ np.array([0.299, 0.587, 0.114])  # BT.601 studio

Q = np.array([[16, 11, 10, 16, 24, 40, 51, 61], [12, 12, 14, 19, 26, 58, 60, 55],
              [14, 13, 16, 24, 40, 57, 69, 56], [14, 17, 22, 29, 51, 87, 80, 62],
              [18, 22, 37, 56, 68, 109, 103, 77], [24, 35, 55, 64, 81, 104, 113, 92],
              [49, 64, 78, 87, 103, 121, 120, 101], [72, 92, 95, 98, 112, 100, 103, 99]])
def jpeg_like(a, scale):  # 8x8 block DCT quantized with scale * (ITU-T T.81 Table K.1)
    b = dctn((a - 128).reshape(64, 8, 64, 8), axes=(1, 3), norm="ortho")
    q = scale * Q[None, :, None, :]
    return idctn(q * np.round(b / q), axes=(1, 3), norm="ortho").reshape(512, 512) + 128

def u8(a):
    return np.round(a).clip(0, 255)
dist = {"Gaussian noise s=10": u8(ref + 10 * rng.standard_normal(ref.shape)),
        "Gaussian blur s=1.5": u8(gaussian_filter(ref, (1.5, 1.5, 0))),
        "block DCT, 2 x Q": u8(np.stack([jpeg_like(ref[..., c], 2) for c in range(3)], -1))}

# 1. Cross-check against scikit-image with the flags that match the paper.
d = dist["Gaussian noise s=10"]
mine = ssim(Y(ref), Y(d))
sk = structural_similarity(Y(ref), Y(d), gaussian_weights=True, sigma=1.5,
                           use_sample_covariance=False, data_range=255)
print(f"paper convention: NumPy {mine:.6f}  skimage {sk:.6f}  |diff| {abs(mine - sk):.1e}")

# 2. One distorted image, many conventions.
cols = ["paper", "same", "downsmp", "sk dflt", "unif7", "RGB", "Ystudio", "MS-SSIM"]
print(f"\n{'distortion':20s}" + "".join(f"{c:>8s}" for c in cols))
for name, d in dist.items():
    yr, yd = Y(ref), Y(d)
    v = [ssim(yr, yd), ssim(yr, yd, valid=False), ssim(yr, yd, downsample=True),
         structural_similarity(yr, yd, data_range=255),  # defaults: 7x7 uniform, N-1
         ssim(yr, yd, win=window("uniform", 7)),
         np.mean([ssim(ref[..., c], d[..., c]) for c in range(3)]),
         ssim(Ys(ref), Ys(d)), ms_ssim(yr, yd)]
    print(f"{name:20s}" + "".join(f"{s:8.4f}" for s in v))

# 3. Data range: the same pair as floats in [0, 1].
r1, d1 = Y(ref) / 255, Y(dist["Gaussian noise s=10"]) / 255
print(f"\n[0,1] data, L=1: {ssim(r1, d1, L=1):.4f}   L=2: {ssim(r1, d1, L=2):.4f}   "
      f"L=255: {ssim(r1, d1, L=255):.4f}")
try:
    structural_similarity(r1, d1)
except ValueError as e:
    print("skimage, float input without data_range: ValueError:", str(e)[:44] + "...")

# 4. Equal MSE: distortions of the grayscale image, clipped to [0, 255] and tuned to the
#    MSE of a one-pixel shift.
g = Y(ref)
mse = lambda a, b: np.mean((a - b) ** 2)
shifted = np.concatenate([g[:, :1], g[:, :-1]], axis=1)
t = mse(g, shifted)
n = rng.standard_normal(g.shape)
order, val = rng.permutation(g.size), rng.integers(0, 2, g.size) * 255.0
def impulses(fr):
    o = g.copy().ravel(); k = int(round(fr * g.size)); o[order[:k]] = val[:k]
    return o.reshape(g.shape)
gm = g.mean()
make = {"Gaussian noise": lambda a: g + a * n,
        "brightness shift": lambda a: g + a,
        "contrast stretch": lambda a: gm + a * (g - gm),
        "Gaussian blur": lambda a: gaussian_filter(g, a),
        "block DCT": lambda a: jpeg_like(g, a),
        "salt and pepper": impulses}
lo_hi = {"Gaussian noise": (1, 60), "brightness shift": (1, 60), "contrast stretch": (1, 3),
         "Gaussian blur": (0.1, 5), "block DCT": (0.1, 30), "salt and pepper": (1e-4, 0.1)}
eq = {k: np.clip(f(brentq(lambda a: mse(g, np.clip(f(a), 0, 255)) - t, *lo_hi[k], xtol=1e-9)),
                 0, 255) for k, f in make.items()}
eq["one-pixel shift"] = shifted
print(f"\n{'distortion':17s}   MSE  SSIM  MS-SSIM  SSIM flat  SSIM texture  SSIM dark")
for name, d in eq.items():
    s, m = ssim(g, d, full=True)  # map index = top-left corner of the 11x11 window;
    # windows inside the blue square (flat), the stripes (texture), the near-black patch
    print(f"{name:17s} {mse(g, d):5.1f} {s:5.3f} {ms_ssim(g, d):8.3f} "
          f"{m[300:440, 270:400].mean():10.3f} {m[45:205, 305:465].mean():13.3f} "
          f"{m[472:495, 305:485].mean():10.3f}")

# 5. 1 - SSIM is not a metric: the triangle inequality fails along a noise path.
x0, x1, x2 = g, g + 3 * n, g + 6 * n
D = lambda a, b: 1 - ssim(a, b)
print(f"\n1-SSIM: d(x0,x1) + d(x1,x2) = {D(x0, x1) + D(x1, x2):.4f} < d(x0,x2) = {D(x0, x2):.4f}")

Output (NumPy 2.3.5, SciPy 1.16.3, scikit-image 0.25.2):

paper convention: NumPy 0.688047  skimage 0.688047  |diff| 5.0e-15

distortion             paper    same downsmp sk dflt   unif7     RGB Ystudio MS-SSIM
Gaussian noise s=10   0.6880  0.6858  0.8949  0.6770  0.6807  0.5447  0.7400  0.9439
Gaussian blur s=1.5   0.8668  0.8707  0.8742  0.8693  0.8695  0.8623  0.8699  0.9150
block DCT, 2 x Q      0.9823  0.9828  0.9926  0.9791  0.9794  0.9763  0.9849  0.9949

[0,1] data, L=1: 0.6880   L=2: 0.8817   L=255: 1.0000
skimage, float input without data_range: ValueError: Since image dtype is floating point, you mus...

distortion          MSE  SSIM  MS-SSIM  SSIM flat  SSIM texture  SSIM dark
Gaussian noise    513.1 0.298    0.720      0.110         0.907      0.187
brightness shift  513.1 0.968    0.998      0.966         0.987      0.492
contrast stretch  513.1 0.928    0.949      0.923         0.877      0.092
Gaussian blur     513.1 0.773    0.842      0.999         0.025      0.990
block DCT         513.1 0.865    0.814      0.994         0.838      0.928
salt and pepper   513.2 0.473    0.753      0.340         0.910      0.473
one-pixel shift   513.1 0.875    0.904      1.000         0.316      1.000

1-SSIM: d(x0,x1) + d(x1,x2) = 0.1590 < d(x0,x2) = 0.2757

Cross-check. The NumPy function and scikit-image with gaussian_weights=True, sigma=1.5, use_sample_covariance=False and data_range=255 agree to rounding error, because scikit-image’s reflect-padded filtering followed by cropping the window radius equals the valid convolution.

Conventions. For the noisy image, the eight columns range from 0.54 to 0.94, a spread far larger than the gaps between methods in published tables such as the one cited above. Noise is the most convention-sensitive distortion because every convention that averages pixels reduces it: the authors’ $2 \times 2$ downsampling quarters the noise variance, and MS-SSIM puts most of its weight on the coarser scales, where the noise is averaged away. Averaging SSIM over RGB channels scores lower than luma, since each channel carries noise of standard deviation 10 while the luma weights combine the three independent noises into about $0.67 \times 10$ (our arithmetic, $\sqrt{0.2989^2 + 0.587^2 + 0.114^2}$). Studio-range luma compresses the signal to 219 levels while $C_2$ stays at its 255-level value, which raises the score. scikit-image’s defaults and the uniform $7 \times 7$ window are within about 0.01 of the paper here, and padding instead of valid windows changes the score by less than 0.004 at this size; with a smaller image the border windows are a larger share of the map (our reasoning). Blur and block quantization move less, by up to about 0.05 (MS-SSIM on the blur).

Data range. With $L = 255$ for data in $[0, 1]$, the constants are 65,025 times too large and every window scores essentially 1; the range of 2 described in the Notes of scikit-image’s documentation for floats inflates the noisy image’s score from 0.69 to 0.88.

Equal MSE. All seven grayscale distortions have the same MSE (salt and pepper misses slightly, since it corrupts whole pixels), and SSIM separates them, ranking the brightness shift and contrast stretch highest and Gaussian noise lowest; the mean squared error article shows why the MSE cannot. SSIM and MS-SSIM do not agree on the order: MS-SSIM ranks the blur above the block DCT, single-scale SSIM the reverse. The region columns show the normalization at work. Noise of the same strength scores far lower in the flat square than in the stripes, the contrast-masking effect, while blur leaves the flat square untouched and destroys the stripes. The brightness shift scores high in the mid-grey square but about half as much in the near-black patch, because the luminance term responds to the relative change. The contrast stretch clips the near-black patch, of luma about 8, to 0, so its luminance term is $C_1/(8^2 + C_1) \approx 0.092$, the printed value. The one-pixel shift is invisible to SSIM in flat regions and costly in the stripes.

Not a metric. Along a path of growing noise, $1 - \mathrm{SSIM}$ from the clean image to the strongly noisy one exceeds the sum of the two half steps, so the triangle inequality fails. The normalized form above explains why: for small errors $1 - \mathrm{SSIM}$ grows roughly with the square of the error amplitude, like the MSE, and squared distances do not satisfy the triangle inequality (our reasoning).

Related Metrics

• Mean squared error and PSNR measure the same pixel errors without local normalization; SSIM is usually reported alongside them.
• LPIPS compares deep-network features of the two images, with weights calibrated on human similarity judgments, and is the third standard metric in view synthesis.
• NIQE scores an image without any reference, for outputs that have no ground truth.
• Fréchet Inception Distance compares distributions of images rather than aligned pairs, for generative models.