BM3D

BM3D (block-matching and 3D filtering) is an image denoising algorithm for additive white Gaussian noise of known standard deviation. For each small block of the image it finds similar blocks nearby, stacks them into a three-dimensional array that its authors call a group, transforms the group with a separable 3D transform, shrinks the transform coefficients and transforms back. This joint processing, which Dabov, Foi, Katkovnik and Egiazarian (2007) named collaborative filtering, gives a filtered version of every block in the group. Because the blocks overlap and each block can join many groups, every pixel receives many estimates, and a weighted average, the aggregation, combines them. The whole procedure runs twice: first with hard thresholding to obtain a basic estimate, then with an empirical Wiener filter whose gains come from that basic estimate.

BM3D builds on the self-similarity that non-local means exploits, but where non-local means takes a weighted average of similar pixels, BM3D uses the similar blocks to make a sparse representation. Its authors reported the highest PSNR they knew of for grayscale and colour Gaussian denoising at the time (Dabov et al., 2007, Section VIII), and Image Denoising describes how it became one of the standard baselines for later methods, including learned ones.

Problem

The observation is $z(x) = y(x) + \eta(x)$ on the pixel grid $X$, where $y$ is the clean image and $\eta$ is independent zero-mean Gaussian noise of variance $\sigma^2$ (Dabov et al., 2007, Section III). The goal is an estimate of $y$ with small mean squared error, usually reported as PSNR.

Transform-domain denoisers assume that the clean signal is well approximated by a few basis functions: keep the few large coefficients, discard the many small ones, and most of the noise goes with them. Dabov et al. (Section I) point out the weakness of a fixed 2D transform: the discrete cosine transform represents sharp transitions poorly, wavelets do poorly on textures and smooth transitions, and no fixed basis is sparse for every kind of image content. Their answer is to make the representation sparser by adding a dimension along which the signal barely changes: a stack of blocks that already look alike.

Inputs and Outputs

• Input: a grayscale image corrupted by additive white Gaussian noise, and the noise standard deviation $\sigma$ (BM3D does not estimate it).
• Output: a denoised image of the same size. The intermediate basic estimate is also available.

Intuition

Suppose a group holds $n$ blocks that are nearly identical. Transformed one at a time by a 2D transform, each block needs about the same number $\alpha$ of significant coefficients, so the group costs about $n\alpha$. If a 1D transform with a constant (DC) basis vector is then applied across the stack, the shared content collects in the DC slice, and the whole group is described by about $\alpha$ significant coefficients instead of $n\alpha$ (Dabov et al., 2007, Section II-D, p. 4, where they call this a qualitative statement). The noise, which an orthonormal transform keeps white, stays spread evenly over all $n \times N_1^2$ coefficients when the blocks do not overlap. Hard thresholding therefore removes much more noise relative to the signal it removes, and the authors note that the sparsity increases with the number of grouped blocks.

Two kinds of correlation are at work. The 2D transform uses the correlation between neighbouring pixels within a block; the 1D transform uses the correlation between corresponding pixels of different blocks. Dabov et al. found the second more important: in their Table II (p. 8), replacing the 2D transform by one whose basis vectors other than DC are random lowered the PSNR by only 0.1–0.4 dB, while keeping only the DC along the stack, which amounts to averaging the grouped blocks, cost much more.

Grouping alone would be weaker. Averaging the blocks of a group is unbiased only when they are truly identical; real groups contain blocks that differ, and the average blurs those differences. The 3D transform keeps them: coefficients beyond the DC slice carry what distinguishes one block from another, and the large ones survive thresholding. This is what “collaborative” means in the paper: each block helps filter the others, while its own features are kept.

Algorithm

Both steps visit reference blocks of size $N_1 \times N_1$ on a grid with spacing $N_{\text{step}}$, and for each one do grouping, collaborative filtering and accumulation; the aggregation divides the accumulated sums at the end (Dabov et al., 2007, Fig. 3 and Sections III and IV).

Step 1: basic estimate.

1. Grouping by block matching. Within an $N_S \times N_S$ window around the reference block, find the blocks whose distance from the reference is at most $\tau^{\text{ht}}_{\text{match}}$. Keep at most $N_2^{\text{ht}}$ of them, nearest first; the reference block, at distance zero, is always first. Stack them into a 3D array.
2. Collaborative hard thresholding. Apply the 3D transform to the group, set to zero every coefficient whose magnitude is at most $\lambda_{3D}\sigma$, and invert the transform. The result is one estimate of every block in the group.
3. Accumulation. Add each block estimate, multiplied by the group’s weight, into a numerator image at the block’s position, and add the weight into a denominator image.

After all reference blocks, the basic estimate is the numerator divided by the denominator.

Step 2: final estimate.

1. Grouping in the basic estimate. Repeat block matching, now on the basic estimate, with threshold $\tau^{\text{wie}}_{\text{match}}$ and at most $N_2^{\text{wie}}$ blocks. Use the matched positions to form two groups: one from the basic estimate and one from the noisy image.
2. Collaborative Wiener filtering. Transform both groups. Use the energy of the basic group’s coefficients as the signal spectrum in a Wiener gain, multiply the noisy group’s coefficients by it, and invert.
3. Accumulation and aggregation as in step 1, with the step 2 weights.

The authors give two reasons for the second step (Section III): matching in the basic estimate is more reliable than matching in the noisy image, and Wiener filtering with the basic estimate as a pilot is more effective than hard thresholding the noisy spectrum.

Mathematical Formulation

The notation follows Dabov et al. (2007). $Z_x$ is the $N_1 \times N_1$ block of $z$ whose top-left corner is $x$, and $\mathbf{Z}_S$ is the 3D array of the blocks at the positions in $S$.

Block distance. With the clean image, the natural distance would be $\lVert Y_{x_R} - Y_x \rVert_2^2 / N_1^2$ (eq. 1). Only noisy blocks are available, and for two non-overlapping blocks the noisy distance

$$d^{\text{noisy}}(Z_{x_R}, Z_x) = \frac{\lVert Z_{x_R} - Z_x \rVert_2^2}{N_1^2}$$

has mean $d^{\text{ideal}} + 2\sigma^2$ and variance $8\sigma^4/N_1^2 + 8\sigma^2 d^{\text{ideal}}/N_1^2$ (eqs. 2 and 3, p. 5). The variance grows as $\sigma^4$, so at high noise or with small blocks the distances of good and bad candidates overlap and blocks are grouped wrongly. As a remedy the paper computes the distance between 2D-transformed blocks after hard-thresholding them at $\lambda_{2D}\sigma$ (eq. 4); its parameter table sets $\lambda_{2D} = 0$ for $\sigma \le 40$, which reduces eq. (4) to the plain distance for an orthonormal transform. The group of a reference block is

$$S^{\text{ht}}_{x_R} = \{\, x \in X : d(Z_{x_R}, Z_x) \le \tau^{\text{ht}}_{\text{match}} \,\} \qquad \text{(eq. 5)}.$$

Collaborative hard thresholding. With $\mathcal{T}_{3D}$ the normalized 3D transform and $\Upsilon$ hard thresholding at $\lambda_{3D}\sigma$,

$$\widehat{\mathbf{Y}}^{\text{ht}}_{S^{\text{ht}}_{x_R}} = \mathcal{T}_{3D}^{-1}\!\left( \Upsilon\!\left( \mathcal{T}_{3D}\!\left( \mathbf{Z}_{S^{\text{ht}}_{x_R}} \right) \right) \right) \qquad \text{(eq. 6)}.$$

$\mathcal{T}_{3D}$ is separable: a 2D transform on each block followed by a 1D transform $\mathcal{T}_{1D}$ along the stack (Section IV). The thresholds only make sense if the transforms are orthonormal, so that white noise of variance $\sigma^2$ stays white with variance $\sigma^2$ in every coefficient. Lebrun (2012, Section 4.1) shows how much this matters: with a DCT whose first coefficient is not normalized, the PSNR in his test was 0.1–1.5 dB lower at $\sigma = 10$–$40$ (and slightly higher at $\sigma = 2$ and 5).

Collaborative Wiener filtering. In step 2 the group is $S^{\text{wie}}_{x_R}$, the positions whose blocks in the basic estimate are within $\tau^{\text{wie}}_{\text{match}}$ of the reference (eq. 7). With $\widehat{\mathbf{Y}}^{\text{basic}}$ the group from the basic estimate, the Wiener coefficients are

$$\mathbf{W}_{S^{\text{wie}}_{x_R}} = \frac{\left|\mathcal{T}^{\text{wie}}_{3D}\big(\widehat{\mathbf{Y}}^{\text{basic}}_{S^{\text{wie}}_{x_R}}\big)\right|^2}{\left|\mathcal{T}^{\text{wie}}_{3D}\big(\widehat{\mathbf{Y}}^{\text{basic}}_{S^{\text{wie}}_{x_R}}\big)\right|^2 + \sigma^2} \qquad \text{(eq. 8)},$$

and the estimate of the group is $\mathcal{T}^{\text{wie}\,-1}_{3D}\big(\mathbf{W} \cdot \mathcal{T}^{\text{wie}}_{3D}(\mathbf{Z}_{S^{\text{wie}}_{x_R}})\big)$ (eq. 9). This is the pointwise Wiener gain $S_k/(S_k + \sigma^2)$ of the Wiener filter applied to each coefficient, with the signal power $S_k$ of that coefficient replaced by the squared coefficient of the basic estimate: an empirical, oracle-guided gain, computed separately for every group.

Aggregation weights. If the noise in a group were independent, the total variance of the step 1 estimate would be $\sigma^2 N_{\text{har}}$, where $N_{\text{har}}$ is the number of coefficients that survive thresholding, and that of the step 2 estimate $\sigma^2 \lVert \mathbf{W} \rVert_2^2$. The weights are the inverses:

$$w^{\text{ht}}_{x_R} = \begin{cases} \dfrac{1}{\sigma^2 N_{\text{har}}}, & N_{\text{har}} \ge 1, \\ 1, & \text{otherwise,} \end{cases} \qquad w^{\text{wie}}_{x_R} = \sigma^{-2} \lVert \mathbf{W}_{S^{\text{wie}}_{x_R}} \rVert_2^{-2} \qquad \text{(eqs. 10, 11)}.$$

A group whose shrinkage zeroed many coefficients, typically one of flat blocks, gets a large weight; a group across an edge, with many surviving coefficients, a small one. Grouped blocks may overlap, so their noise is not independent, and the authors call the weights only loosely inversely proportional to the variances (pp. 6–7). The final image is the weighted average of all block estimates covering each pixel (eq. 12):

$$\hat{y}(x) = \frac{\sum_{x_R} \sum_{x_m \in S_{x_R}} w_{x_R}\, \widehat{Y}^{x_R}_{x_m}(x)}{\sum_{x_R} \sum_{x_m \in S_{x_R}} w_{x_R}\, \chi_{x_m}(x)},$$

where $\chi_{x_m}$ is 1 on the support of the block at $x_m$ and 0 elsewhere, and each block estimate is zero outside its support. In the paper’s implementation the weights in eq. (12) also include an $N_1 \times N_1$ Kaiser window, to reduce border effects from transforms such as the 2D DCT (Section IV).

Parameters

Dabov et al. (2007) give two profiles in Table I (p. 8). The Normal profile has two columns:

Parameter	$\sigma \le 40$	$\sigma > 40$
Step 1 2D transform $\mathcal{T}^{\text{ht}}_{2D}$	2D Bior1.5	2D DCT
Block size $N_1^{\text{ht}}$	8	12
Maximum group size $N_2^{\text{ht}}$	16	16
Reference step $N^{\text{ht}}_{\text{step}}$	3	4
Search window $N_S^{\text{ht}}$	39	39
Distance prefilter $\lambda_{2D}$	0	2.0
Threshold factor $\lambda_{3D}$	2.7	2.8
Matching threshold $\tau^{\text{ht}}_{\text{match}}$	2500	5000
Step 2 2D transform $\mathcal{T}^{\text{wie}}_{2D}$	2D DCT	2D DCT
Block size $N_1^{\text{wie}}$	8	11
Maximum group size $N_2^{\text{wie}}$	32	32
Reference step $N^{\text{wie}}_{\text{step}}$	3	6
Search window $N_S^{\text{wie}}$	39	39
Matching threshold $\tau^{\text{wie}}_{\text{match}}$	400	3500
Kaiser window parameter $\beta$	2.0	2.0
1D transform $\mathcal{T}_{1D}$	1D Haar	1D Haar

Distances and thresholds are on the 0–255 scale. Because the distance is normalized by $N_1^2$, $\tau^{\text{ht}}_{\text{match}} = 2500$ is a mean squared difference of 2500 per pixel, about 50 grey levels root mean square, against the $2\sigma^2 = 1250$ that noise alone adds at $\sigma = 25$ (our arithmetic). The Haar transform needs a power of 2 along the stack, so the number of grouped blocks is cut to the largest power of 2 not exceeding the number of matches (p. 9). The order of blocks within a group hardly matters, the authors found, and they keep the order of increasing distance (p. 9). The Fast profile searches only near the previous matches and, at moderate noise ($\sigma < 35$), lost 0.05–0.2 dB against the Normal profile for a sixfold speed-up (Section VI-A).

Lebrun (2012) re-examined the parameters in his IPOL implementation and kept most of the paper’s values; his summary table (Section 3.7, p. 189) uses 8 × 8 blocks at all noise levels, groups of at most 16 and 32, $\lambda_{3D} = 2.7$, matching thresholds 2500 and 400 for $\sigma \le 40$ (5000 and 3500 above), and a reference step of 3 in both steps. He found that small changes of $\lambda_{3D}$ have a strong influence (Section 3.3): 2.5 or 3.2 instead of 2.7 lowered the PSNR on his test image at every $\sigma$ from 5 to 100, by up to 0.44 dB at $\sigma = 100$. The matching thresholds can cost several decibels if chosen wrongly, but a fixed pair works over a wide range of $\sigma$ (Section 3.4). The maximum group sizes had little influence (Section 3.2). His implementation uses a Walsh–Hadamard transform along the stack, where the paper uses a Haar decomposition; for his transform comparison (Table 1, p. 196) the best combination was Bior1.5 in step 1 and the 2D DCT in step 2, each with Walsh–Hadamard along the stack.

Complexity

With every parameter fixed, the cost is linear in the number of pixels (Section IV-A). The paper’s operation count per pixel has three terms: pre-computing the 2D transforms of the blocks in each search window, the exhaustive block matching in that window, and the 3D transforms of the groups; the last two shrink with the square of the reference step $N_{\text{step}}$. The paper reports about 4.1 s for a 256 × 256 image with the Normal profile ($\sigma \le 40$) on a 1.5 GHz Celeron M, and 0.7 s with the Fast profile (Table I); block matching took about half of the execution time (Section V). The reference blocks are independent, so the method parallelizes well, which Lebrun notes in his conclusions (Section 8).

Implementation

The script below implements both steps in NumPy and compares them with two step 1 variants and with two scikit-image denoisers. It uses NumPy 2.3.5, SciPy 1.16.3, scikit-image 0.25.2 and PyWavelets (pywt.__version__ 1.8.0), with the BLAS thread count set to one, runs in about 20 seconds, and printed identical output on two runs. The bm3d package that wraps the authors’ own software is not used: that software is distributed under a licence that binds on access, and this article neither installs nor reads it.

The implementation follows the Normal profile for $\sigma \le 40$ (Table I): blocks of 8 × 8, groups of at most 16 in step 1 and 32 in step 2, a reference step of 3, a 39 × 39 search window, $\lambda_{2D} = 0$, $\lambda_{3D} = 2.7$, matching thresholds 2500 and 400, a 2D DCT in step 2 and an orthonormal Haar decomposition along the stack, with the eq. (10) and (11) weights. It departs from the paper in five places, all fixed before the first run:

• D1. Step 1 uses the orthonormal 2D DCT rather than Bior1.5. In the paper’s Table II, the step 1 transform changes the basic estimate’s PSNR by at most 0.21 dB on Boats and 0.23 dB on Lena across the structured transforms, with the DCT within 0.1 dB of Bior1.5 on both.
• D2. There is no Kaiser window. Lebrun (2012, Section 4.4) found that it made no difference to the PSNR with a reference step of 3.
• D3. If every Wiener coefficient of a group is zero, eq. (11) is infinite; the code gives such a group weight 1. The paper has no rule for this case. The guard was added during planning, after a first run of the checks returned NaN on the noise-free astronaut crop, which contains blocks of exact zeros, and before the example’s plan was written.
• D4. The step 2 match uses $\le \tau^{\text{wie}}_{\text{match}}$, where eq. (7) has a strict inequality.
• D5. Choices the paper leaves open: reference blocks start at 0 every 3 pixels, plus the last row and column position so every pixel is covered; the search window is centred on the reference block’s top-left corner and clipped at the image border.

Block distances are computed between 2D DCT spectra, which for an orthonormal transform equals the pixel-domain distance (check C4 confirms it). Two variants of step 1 isolate the two halves of collaborative filtering. No grouping uses groups of one block, which is sliding-window DCT hard thresholding with the same threshold and weights. Group mean uses the same matched groups as step 1 but keeps only the DC slice along the stack, unthresholded, and discards the rest; with an orthonormal 2D transform this is the elementwise average of the grouped blocks, the “DC-only” 1D transform of Dabov et al.’s Table II and the “basic NL-means” of Lebrun’s Section 4.1.

The data are the top-left 256 × 256 crop of scikit-image’s bundled astronaut image, converted to grey levels 0–255 (flag stripes and stars, hair, background and the edge of the face), with Gaussian noise of $\sigma = 25$ for seeds 0, 1 and 2, neither clipped nor rounded. A synthetic 64 × 64 image made of one random 8 × 8 tile repeated eight times in each direction supplies exact copies for the checks. The two library baselines receive the true noise level. denoise_wavelet runs with its defaults (BayesShrink, soft thresholding, the db1 wavelet). denoise_nl_means uses its default patch size 7 and patch distance 11, with h = 0.8 $\sigma$, which its docstring gives as a reasonable starting point when sigma is supplied in fast mode; these settings were fixed before the run. PSNR is $10 \log_{10}(255^2/\text{MSE})$, as in Dabov et al. (p. 10).

Before the first run we fixed seven checks with their tolerances, from a planning run of the same checks on the same data: the hand-written DCT matches SciPy’s (C1) and the 3D transform is orthonormal (C2), both within $10^{-12}$; block matching on the periodic image groups only exact copies, with group sizes from 8 to 16 as the window geometry implies (C3, within $10^{-9}$); the matching distances agree with a direct pixel loop (C4, within $10^{-9}$); every variant and step 2 return a noise-free input unchanged when $\sigma = 10^{-6}$, within $10^{-4}$, a tolerance set from the largest change thresholding can make at that $\sigma$ (C5; the group mean is excluded on the astronaut crop, where averaging non-identical blocks is not the identity); group sizes are powers of 2 with the reference first (C6); and the noise standard deviation of each seed is within four standard errors, 0.2762, of 25 (C7).

We also stated hypotheses. On sparsity (H1), counting the clean-signal coefficients above the step 1 threshold $2.7\sigma = 67.5$, per block, in the groups that step 1’s matching on the noisy seed 0 image fills to 16 blocks, taking their first $n$ blocks: for matched groups the count falls strictly as the group grows from 1 to 2, 4, 8 and 16 blocks (H1a), and at every group size of 2 or more a control group of the reference plus randomly drawn blocks from the same window has a higher count than the matched group (H1b). On denoising, on every seed: the full step 1 beats both variants (H2a); step 2 improves on step 1 (H2b); the final estimate beats both library baselines (H2c), and so does step 1 alone (H2d); and every method improves on the noisy input (H2e).

import os
for v in ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS",
          "VECLIB_MAXIMUM_THREADS", "NUMEXPR_NUM_THREADS"):
    os.environ[v] = "1"

import numpy as np
import scipy
import scipy.fft
import skimage
import pywt
from skimage import color, data
from skimage.restoration import denoise_nl_means, denoise_wavelet

# ---------------- transforms ----------------
K1 = 8                                          # block size N1 (Table I)


def dct_matrix(n):
    """Orthonormal DCT-II matrix; rows are basis vectors."""
    k = np.arange(n)[:, None]
    i = np.arange(n)[None, :]
    D = np.cos(np.pi * (2 * i + 1) * k / (2 * n)) * np.sqrt(2.0 / n)
    D[0] /= np.sqrt(2.0)
    return D


def haar_matrix(n):
    """Orthonormal full dyadic Haar transform of length n (a power of 2)."""
    if n == 1:
        return np.ones((1, 1))
    h = haar_matrix(n // 2)
    top = np.kron(h, [1.0, 1.0])               # coarse part, recursively
    bottom = np.kron(np.eye(n // 2), [1.0, -1.0])  # finest differences
    return np.vstack([top, bottom]) / np.sqrt(2.0)


D8 = dct_matrix(K1)
T2D = np.kron(D8, D8)                           # 2D DCT of a row-major flattened 8x8 block
HAAR = {n: haar_matrix(n) for n in (1, 2, 4, 8, 16, 32)}


def block_spectra(img):
    """2D DCT spectra of every 8x8 block, indexed by its top-left corner."""
    win = np.lib.stride_tricks.sliding_window_view(img, (K1, K1))
    return win.reshape(win.shape[0], win.shape[1], K1 * K1) @ T2D.T


def ref_grid(n, step):
    g = list(range(0, n - K1 + 1, step))
    if g[-1] != n - K1:
        g.append(n - K1)                         # make sure the last row/column is covered
    return g


def match(spec, i, j, n_max, tau, ns=39):
    """Block matching: positions of up to n_max blocks within the ns x ns window
    whose distance ||.||^2 / N1^2 to the reference is <= tau, nearest first,
    cut to the largest power of 2. Returns positions and distances."""
    h, w = spec.shape[:2]
    r = ns // 2
    i0, i1 = max(i - r, 0), min(i + r, h - 1)
    j0, j1 = max(j - r, 0), min(j + r, w - 1)
    cand = spec[i0:i1 + 1, j0:j1 + 1].reshape(-1, K1 * K1)
    d = np.sum((cand - spec[i, j]) ** 2, axis=1) / K1**2
    ii, jj = np.meshgrid(np.arange(i0, i1 + 1), np.arange(j0, j1 + 1), indexing="ij")
    ii, jj = ii.ravel(), jj.ravel()
    d_sort = d.copy()
    d_sort[(ii == i) & (jj == j)] = -1.0         # the reference block comes first
    order = np.argsort(d_sort, kind="stable")
    order = order[d[order] <= tau][:n_max]
    n = 1 << (len(order).bit_length() - 1)       # largest power of 2 <= count
    order = order[:n]
    return ii[order], jj[order], d[order]


def aggregate(shape, pos, blocks, weight, num, den):
    for (a, b), blk in zip(pos, blocks):
        num[a:a + K1, b:b + K1] += weight * blk.reshape(K1, K1)
        den[a:a + K1, b:b + K1] += weight


def step1(z, sigma, mode="3d", n2=16, tau=2500.0, lam=2.7, step=3, stats=None):
    """BM3D step 1 (hard thresholding). mode: '3d' (BM3D), 'dc' (keep only the
    DC along the group, i.e. average the matched blocks), '2d' (no grouping)."""
    spec = block_spectra(z)
    num = np.zeros_like(z)
    den = np.zeros_like(z)
    for i in ref_grid(z.shape[0], step):
        for j in ref_grid(z.shape[1], step):
            if mode == "2d":
                pi, pj = np.array([i]), np.array([j])
            else:
                pi, pj, _ = match(spec, i, j, n2, tau)
            n = len(pi)
            if stats is not None:
                stats.append((i, j, pi, pj))
            g = HAAR[n] @ spec[pi, pj]           # 3D spectrum: 1D Haar along the stack
            if mode == "dc":
                g[1:] = 0.0                      # discard everything but the DC slice
            else:
                g[np.abs(g) <= lam * sigma] = 0.0
            n_har = np.count_nonzero(g)
            w = 1.0 / (sigma**2 * n_har) if n_har >= 1 else 1.0
            est = (HAAR[n].T @ g) @ T2D          # inverse 1D Haar, inverse 2D DCT
            aggregate(z.shape, zip(pi, pj), est, w, num, den)
    return num / den


def step2(z, basic, sigma, n2=32, tau=400.0, step=3):
    """BM3D step 2: matching in the basic estimate, empirical Wiener filtering."""
    spec_z = block_spectra(z)
    spec_b = block_spectra(basic)
    num = np.zeros_like(z)
    den = np.zeros_like(z)
    for i in ref_grid(z.shape[0], step):
        for j in ref_grid(z.shape[1], step):
            pi, pj, _ = match(spec_b, i, j, n2, tau)
            n = len(pi)
            gb = HAAR[n] @ spec_b[pi, pj]
            gz = HAAR[n] @ spec_z[pi, pj]
            wien = gb**2 / (gb**2 + sigma**2)
            nw = np.sum(wien**2)
            w = 1.0 / (sigma**2 * nw) if nw > 0 else 1.0   # guard: all-zero basic group
            est = (HAAR[n].T @ (wien * gz)) @ T2D
            aggregate(z.shape, zip(pi, pj), est, w, num, den)
    return num / den


def psnr(x, ref):
    return 10 * np.log10(255.0**2 / np.mean((x - ref) ** 2))


print(f"NumPy {np.__version__}, SciPy {scipy.__version__}, "
      f"scikit-image {skimage.__version__}, PyWavelets {pywt.__version__}")

# ---------------- data ----------------
SIGMA = 25.0
clean = 255.0 * color.rgb2gray(data.astronaut())[:256, :256]   # flag, hair, face
rng_p = np.random.default_rng(123)
tile = rng_p.uniform(0, 255, (K1, K1))
periodic = np.tile(tile, (8, 8))                                  # 64 x 64, period 8

# ---------------- A. checks on exact data ----------------
print("\nA. checks")
x = np.random.default_rng(7).standard_normal((K1, K1))
c1 = np.max(np.abs((T2D @ x.ravel()).reshape(K1, K1) - scipy.fft.dctn(x, norm="ortho")))
print(f"C1 2D DCT matrix vs scipy.fft.dctn(norm='ortho'): max diff {c1:.1e}")
c2 = 0.0
for n, H in HAAR.items():
    T3 = np.kron(H, T2D)
    c2 = max(c2, np.max(np.abs(T3 @ T3.T - np.eye(T3.shape[0]))))
print(f"C2 3D transform (2D DCT x 1D Haar, n = 1..32) orthonormal: max |T T' - I| {c2:.1e}")

# C3: block matching on the periodic image finds only exact copies
st = []
_ = step1(periodic, 1e-6, stats=st)
spec_p = block_spectra(periodic)
dmax, sizes = 0.0, []
for (i, j, pi, pj) in st:
    dmax = max(dmax, np.max(np.abs(spec_p[pi, pj] - spec_p[i, j])))
    sizes.append(len(pi))
print(f"C3 periodic image: max |grouped - reference| coefficient {dmax:.1e}, "
      f"group sizes {min(sizes)}-{max(sizes)}")

# C4: matching distances against a direct pixel loop (noisy crop, seed 0)
noisy0 = clean + SIGMA * np.random.default_rng(0).standard_normal(clean.shape)
spec0 = block_spectra(noisy0)
rr = np.random.default_rng(11)
c4 = 0.0
for _ in range(20):
    i, j = rr.integers(0, 249, 2)
    pi, pj, d = match(spec0, i, j, 16, 2500.0)
    for a, b, dd in zip(pi, pj, d):
        direct = np.sum((noisy0[a:a + K1, b:b + K1] - noisy0[i:i + K1, j:j + K1]) ** 2) / K1**2
        c4 = max(c4, abs(direct - dd))
print(f"C4 matching distance vs direct pixel loop (20 references): max diff {c4:.1e}")

# C5: noise-free inputs are returned unchanged (sigma = 1e-6)
for name, img in (("periodic", periodic), ("astronaut crop", clean)):
    out = {"3d": step1(img, 1e-6), "2d": step1(img, 1e-6, mode="2d")}
    if name == "periodic":
        out["dc"] = step1(img, 1e-6, mode="dc")
    out["step 2"] = step2(img, img, 1e-6)
    print(f"C5 {name:15s} noise-free, max |out - in|: " +
          ", ".join(f"{k} {np.max(np.abs(v - img)):.1e}" for k, v in out.items()))

# C6: group sizes on the noisy crop are powers of 2, reference first
st = []
_ = step1(noisy0, SIGMA, stats=st)
sz = np.array([len(p[2]) for p in st])
first = all(p[2][0] == p[0] and p[3][0] == p[1] for p in st)
pow2 = bool(np.all((sz & (sz - 1)) == 0))
print(f"C6 noisy crop, step 1: {len(st)} references, group sizes {sz.min()}-{sz.max()}, "
      f"powers of 2: {pow2}, reference first: {first}, size 16: {np.mean(sz == 16):.3f}")

# C7: the noise has the planned standard deviation
seeds = (0, 1, 2)
noisy = {s: clean + SIGMA * np.random.default_rng(s).standard_normal(clean.shape) for s in seeds}
print("C7 noise std per seed: " + ", ".join(f"{np.std(noisy[s] - clean):.4f}" for s in seeds))

# ---------------- B. sparsity of grouped against single blocks ----------------
print("\nB. clean-signal coefficients above 2.7 sigma, per block (seed 0 grouping)")
spec_c = block_spectra(clean)
full = [p for p in st if len(p[2]) == 16]
rs = np.random.default_rng(5)
thr = 2.7 * SIGMA
rows = {"matched": [], "random": []}
for n in (1, 2, 4, 8, 16):
    cm, cr = [], []
    for (i, j, pi, pj) in full:
        gm = HAAR[n] @ spec_c[pi[:n], pj[:n]]
        cm.append(np.count_nonzero(np.abs(gm) > thr) / n)
        i0, i1 = max(i - 19, 0), min(i + 19, 248)
        j0, j1 = max(j - 19, 0), min(j + 19, 248)
        ri = np.concatenate([[i], rs.integers(i0, i1 + 1, n - 1)])
        rj = np.concatenate([[j], rs.integers(j0, j1 + 1, n - 1)])
        gr = HAAR[n] @ spec_c[ri, rj]
        cr.append(np.count_nonzero(np.abs(gr) > thr) / n)
    rows["matched"].append(np.mean(cm))
    rows["random"].append(np.mean(cr))
print(f"  {len(full)} references with 16 matched blocks")
print("  group size n        1      2      4      8     16")
for k, v in rows.items():
    print(f"  {k:10s}     " + " ".join(f"{a:6.3f}" for a in v))

# ---------------- C. denoising ----------------
print("\nC. PSNR (dB), sigma = 25, astronaut crop 256 x 256")
methods = ["noisy", "wavelet (skimage)", "NL-means (skimage)", "step 1, no grouping",
           "step 1, group mean", "step 1 (basic)", "steps 1 + 2 (final)"]
res = {m: [] for m in methods}
for s in seeds:
    z = noisy[s]
    res["noisy"].append(psnr(z, clean))
    res["wavelet (skimage)"].append(psnr(255 * denoise_wavelet(z / 255, sigma=SIGMA / 255), clean))
    res["NL-means (skimage)"].append(psnr(255 * denoise_nl_means(
        z / 255, patch_size=7, patch_distance=11, h=0.8 * SIGMA / 255,
        sigma=SIGMA / 255, fast_mode=True), clean))
    res["step 1, no grouping"].append(psnr(step1(z, SIGMA, mode="2d"), clean))
    res["step 1, group mean"].append(psnr(step1(z, SIGMA, mode="dc"), clean))
    basic = step1(z, SIGMA)
    res["step 1 (basic)"].append(psnr(basic, clean))
    res["steps 1 + 2 (final)"].append(psnr(step2(z, basic, SIGMA), clean))
print("  method                  seed 0   seed 1   seed 2     mean")
for m in methods:
    v = res[m]
    print(f"  {m:22s} " + " ".join(f"{a:8.3f}" for a in v) + f" {np.mean(v):8.3f}")

Output:

NumPy 2.3.5, SciPy 1.16.3, scikit-image 0.25.2, PyWavelets 1.8.0

A. checks
C1 2D DCT matrix vs scipy.fft.dctn(norm='ortho'): max diff 4.1e-15
C2 3D transform (2D DCT x 1D Haar, n = 1..32) orthonormal: max |T T' - I| 2.7e-15
C3 periodic image: max |grouped - reference| coefficient 0.0e+00, group sizes 8-16
C4 matching distance vs direct pixel loop (20 references): max diff 3.0e-12
C5 periodic        noise-free, max |out - in|: 3d 8.2e-13, 2d 8.8e-13, dc 8.2e-13, step 2 9.4e-13
C5 astronaut crop  noise-free, max |out - in|: 3d 7.0e-09, 2d 9.1e-10, step 2 6.5e-10
C6 noisy crop, step 1: 7056 references, group sizes 1-16, powers of 2: True, reference first: True, size 16: 0.960
C7 noise std per seed: 24.9860, 24.8986, 24.9249

B. clean-signal coefficients above 2.7 sigma, per block (seed 0 grouping)
  6774 references with 16 matched blocks
  group size n        1      2      4      8     16
  matched         1.854  1.161  0.762  0.549  0.466
  random          1.854  1.814  1.801  1.771  1.779

C. PSNR (dB), sigma = 25, astronaut crop 256 x 256
  method                  seed 0   seed 1   seed 2     mean
  noisy                    20.177   20.207   20.198   20.194
  wavelet (skimage)        26.605   26.620   26.616   26.614
  NL-means (skimage)       29.190   29.107   29.243   29.180
  step 1, no grouping      29.171   29.052   29.064   29.096
  step 1, group mean       29.455   29.361   29.483   29.433
  step 1 (basic)           30.501   30.495   30.563   30.520
  steps 1 + 2 (final)      31.052   31.008   31.084   31.048

All seven checks held. The DCT, the 3D transform and the matching distances agree with their references to between $10^{-15}$ and $10^{-12}$, so C1, C2 and C4 held. On the periodic image every grouped block equals its reference exactly, with 8 to 16 blocks per group, so C3 held, and every variant returns both noise-free images to within $7 \times 10^{-9}$, far inside the $10^{-4}$ tolerance, so C5 held: the aggregation reproduces its input when nothing is removed, and the group mean is exact when the grouped blocks are identical. Every group size is a power of 2 with the reference first; 96% of the step 1 groups reach the full 16 blocks, and the smallest hold a single block because nothing else in the window is within the matching threshold; C6 held. The noise standard deviations, 24.8986 to 24.9860, are within 0.11 of 25, so C7 held.

Part B shows the sparsity gain directly, and H1a and H1b held. A single 8 × 8 block of the clean image has on average 1.854 DCT coefficients above 67.5. In the matched groups the count per block falls with every doubling, to 0.466 at 16 blocks: a group of 16 blocks needs about 7.5 significant coefficients, against about 30 for 16 blocks at the single-block rate (our arithmetic from the printed values). With random blocks from the same window the count stays near 1.8 per block at every size. In our reading, the DC slice of a random group still collects the blocks’ shared mean brightness, but their differences are as large as the blocks themselves and spread over the other slices; matching is what makes the stack nearly constant. Not hypothesized: the random row is not monotone, with 1.771 at 8 blocks and 1.779 at 16.

Part C held for every hypothesis on every seed. Step 1 alone reaches 30.50–30.56 dB. Without grouping, sliding-window DCT thresholding gives 29.05–29.17 dB; averaging the matched blocks gives 29.36–29.48 dB; transforming and thresholding the same groups adds 1.05–1.13 dB over averaging them. H2a held: the full step 1 is above both variants on every seed. The empirical Wiener step adds 0.51–0.55 dB on each seed, so H2b held, close to the “typically greater than 0.5 dB” that Dabov et al. read from their Table II (Section VII). Both BM3D estimates are above scikit-image’s non-local means (29.11–29.24 dB) and its BayesShrink wavelet denoiser with the default Haar wavelet (26.61–26.62 dB) on every seed, so H2c and H2d held, and every method improves on the noisy 20.18–20.21 dB, so H2e held. Not hypothesized: scikit-image’s non-local means is above the no-grouping variant on all three seeds, by 0.02 to 0.18 dB. The ranking is for this crop, this noise level and these settings of the baselines; Dabov et al. and Lebrun report their comparisons on their own test images.

Properties and Behavior

• Sparsity through grouping. The 1D transform across matched blocks, not the 2D transform within them, carries most of the gain. Dabov et al.’s Table II shows that the choice of 2D transform hardly matters, that a transform along the stack without a DC element (the DST) does worse, and that keeping only the DC along the stack, plain averaging, loses about 1.9 dB on Boats and 1.4 dB on Lena at $\sigma = 25$ against the Haar transform.
• Two steps. The second step uses the basic estimate twice, for better matching and as the signal spectrum of the Wiener gain. In Lebrun’s table (Section 4.2, p. 192) it raised the PSNR by 0.1 dB at $\sigma = 2$ and by about 1 dB at $\sigma = 40$ on his test image.
• Weights favour sparse groups. Groups with few surviving coefficients receive large weights, which reduces artefacts near edges (Lebrun, 2012, Section 2.2). In his tables the weighting changed the PSNR by at most about 0.4 dB in step 1 and very little in step 2, while being visually useful in step 1 (Sections 4.1, 4.2 and 8).
• Ordering does not matter. A random orthonormal basis along the stack, apart from its DC vector, performed like the structured transforms, which the authors read as evidence that the order of blocks in the group is irrelevant (p. 9).
• Near the oracle’s reach. Lebrun compared step 2 with an ideal Wiener filter that takes its gains from the clean image: the ideal filter was 2.3–3.3 dB better at $\sigma = 2$–$40$, which his text summarizes as “some 3 dB” (Section 4.2). NL-Bayes later recast both steps as Gaussian estimates for patch groups, as Wiener Filter describes.
• Reported results. Dabov et al. report grayscale PSNR on 12 standard images at $\sigma$ from 2 to 100 (Table III; for example 29.45 dB on Cameraman and 32.86 dB on House at $\sigma = 25$), better than every method they compared in Fig. 4. Lebrun compared his implementation with BLS-GSM, K-SVD, non-local means, DCT denoising and total variation denoising on six images at $\sigma = 2$–$40$. His text says BM3D gave the best results for every noise level and image, but his tables show K-SVD slightly ahead in three cells, Flowers and Traffic at $\sigma = 5$ and Traffic at $\sigma = 10$ (Tables 2 and 3).

Limitations

• Known, white, Gaussian noise. The thresholds and the Wiener gain assume i.i.d. Gaussian noise of known $\sigma$; the method neither estimates $\sigma$ nor handles signal-dependent or correlated noise as it stands. Dabov et al. note that other noise models need changes to how the coefficient variances are computed (Section VIII), and Mäkinen, Azzari and Foi (2020) argue that earlier collaborative filters of correlated noise approximated the transform-domain noise spectrum in ways that are inaccurate when the noise in one block is correlated with that in another. Real camera noise is signal-dependent, as the Poisson–Gaussian model in Image Noise describes.
• Grouping fails at high noise. Because the distance variance grows as $\sigma^4$, matching degrades at high noise; the $\sigma \le 40$ parameters showed a sharp drop in PSNR at about $\sigma = 40$, which the $\sigma > 40$ parameters, with prefiltered distances and larger blocks, avoid (Fig. 9 and footnote 2). Lebrun notes “many artifacts” at $\sigma > 30$ (Section 8).
• Texture and detail. Lebrun reports that the method often flattens micro-textures, such as skin, and suggests adding back a little noise to restore the grain (Section 8).
• Rare structures. In our reading, as with non-local means, a block with no similar blocks in its window forms a small group and gains little from collaboration; in the example some reference blocks were left with groups of one.
• Cost. It is slower and more complex than local transform thresholding (Lebrun, Section 8), although its cost is linear in the image size.

Variants

• Colour (C-BM3D). The image is converted to a luminance–chrominance space; grouping is done once on the luminance, which has the higher signal-to-noise ratio, and the same groups are used to filter all three channels (Dabov et al., 2007, Section V). In their Table VI this beat running BM3D independently on the opponent colour channels by 0.2–0.4 dB and on the RGB channels by 0.3–0.8 dB.
• Fast profile. Predictive search restricts matching to small neighbourhoods of the previous matches, with an exhaustive search every few blocks (Table I).
• BM3D-SAPCA. Lebrun (2012, Section 1) describes a later version by the same group that uses adaptive-shape neighbourhoods instead of square blocks and a PCA basis estimated from each group’s second-moment matrix.
• Volumetric data (BM4D). Maggioni, Katkovnik, Egiazarian and Foi (2013) stack cubes of voxels into four-dimensional groups and filter them with a 4D transform, for volumetric data with Gaussian or Rician noise, such as MRI, and also use the filter as a regularizer for reconstruction from incomplete measurements.
• Correlated noise. Mäkinen, Azzari and Foi (2020) compute the exact transform-domain noise variance of a group under stationary correlated noise and use it in shrinkage, patch matching and aggregation; the paper points to the authors’ BM3D software, distributed on their page and on PyPI; as of November 2026 the page offers version 4.0.3.
• Despeckling. SAR-BM3D is among the nonlocal despeckling filters for radar images that synthetic aperture radar lists.
• Restoration beyond denoising. Lebrun (2012, Section 1) notes that later work generalized BM3D into a restoration tool that includes deblurring; Image Deblurring covers that problem.
• Patch-group successors. NL-Bayes models each group as a Gaussian vector (see Wiener Filter), and weighted nuclear norm minimization treats a group as a low-rank matrix; Image Denoising summarizes weighted nuclear norm minimization and the learned denoisers, such as DnCNN, that later reported higher PSNR than BM3D.