Non-Local Means

Non-local means (NL-means) is an image denoising algorithm that estimates each pixel as a weighted average of other pixels, giving large weights to pixels whose surrounding patches resemble the patch around the pixel being denoised. Similar patches are searched for across a large region rather than only among immediate neighbors, hence the name. Introduced by Buades, Coll and Morel in 2005, it helped make image self-similarity a central idea of later denoising methods such as BM3D.

Problem

The input is an image corrupted by additive noise, usually white Gaussian noise of known standard deviation $\sigma$ (see Image Noise). Local filters such as the Gaussian filter average pixels because they are close, and so blur edges and texture. The bilateral filter also requires similar values, but compares single noisy pixel values, which works only when edges are much stronger than the noise. NL-means needs a more reliable test of which pixels share the same underlying intensity.

Inputs and Outputs

• Input: a noisy grayscale or color image, the noise standard deviation $\sigma$, a patch size, a search window size, and a filtering parameter $h$.
• Output: a denoised image of the same size.

Intuition

Natural images are redundant: a small patch, say $7 \times 7$ pixels, usually has many near copies in the same image. Along an edge, every patch centered on the edge resembles its neighbors along the edge direction. In a periodic texture, every patch resembles patches one period away. In a flat region, almost every patch looks alike.

If two patches look alike, their center pixels probably had the same clean value, and averaging them reduces the noise. Comparing whole patches rather than single values is what makes this test reliable: noise changes the distance between two patches by roughly the same amount for every pair, so the most similar patches in the noisy image tend to be the most similar in the clean one. A pixel on an edge is averaged with other pixels on the same side of the same edge, wherever they are, and not with the pixels across it.

Algorithm

For each pixel $p$:

1. Take the patch of size $(2f + 1) \times (2f + 1)$ centered on $p$.
2. For each pixel $q$ in a search window of size $(2r + 1) \times (2r + 1)$ around $p$, compute the distance between the patch at $p$ and the patch at $q$.
3. Convert each distance to a weight that decreases quickly with distance, controlled by $h$.
4. Output the weighted average of the values at all $q$, divided by the sum of the weights.

As with the bilateral filter, every output is computed from the noisy input in a single pass.

Mathematical Formulation

Let $v$ be the noisy image and $\mathcal{N}_p$ the patch around pixel $p$. Buades, Coll and Morel (2005) define the estimate

$$\mathrm{NL}[v](p) = \sum_{q} w(p, q)\, v(q), \qquad w(p, q) = \frac{1}{Z(p)} \exp\!\left(-\frac{\| v(\mathcal{N}_p) - v(\mathcal{N}_q) \|_{2,a}^2}{h^2}\right),$$

with the normalizing factor

$$Z(p) = \sum_{q} \exp\!\left(-\frac{\| v(\mathcal{N}_p) - v(\mathcal{N}_q) \|_{2,a}^2}{h^2}\right),$$

so the weights are between 0 and 1 and sum to one. The patch distance $\| \cdot \|_{2,a}^2$ is a squared Euclidean distance in which each pixel offset within the patch is weighted by a Gaussian of standard deviation $a$, so that pixels near the patch center count more. In principle the sum runs over the whole image; in practice it is restricted to the search window.

For additive white noise of variance $\sigma^2$, with the patch weights summing to one, the expected distance between two noisy patches is the distance between the clean patches plus a constant:

$$\mathbb{E}\,\| v(\mathcal{N}_p) - v(\mathcal{N}_q) \|_{2,a}^2 = \| u(\mathcal{N}_p) - u(\mathcal{N}_q) \|_{2,a}^2 + 2\sigma^2 ,$$

where $u$ is the clean image. This is why the patch distance preserves the order of similarity on average.

The reference implementation by the same authors in Image Processing On Line (2011) uses a box-weighted patch distance $d^2$, the mean over the patch pixels (and color channels) of the squared differences, and subtracts the expected noise contribution:

$$w(p, q) = \exp\!\left(-\frac{\max\!\left(d^2(p, q) - 2\sigma^2,\; 0\right)}{h^2}\right).$$

Patches whose distance is within what noise alone would produce receive weight 1. Because $w(p, p) = 1$ would always be the largest weight, the pixel’s own weight is set instead to the largest weight among the other pixels in the search window, which avoids over-weighting the noisy center value.

Parameters

• Patch size $(2f+1)^2$. Larger patches make the distance more robust to noise but find fewer close matches for small details. Buades, Coll and Morel (2005) recommended $7 \times 7$ or $9 \times 9$ for grayscale images and $5 \times 5$ or $3 \times 3$ for color images with little noise. The parameter table of the IPOL article (given for its patchwise implementation) grows the patch with the noise: for grayscale images, from $3 \times 3$ at $\sigma \le 15$ (on a 0–255 scale) to $11 \times 11$ at $\sigma$ up to 100. In the limit of a single-pixel patch, the method reduces to the Yaroslavsky neighborhood filter.
• Search window $(2r+1)^2$. Larger windows find more similar patches but cost more. Buades, Coll and Morel used $21 \times 21$; the IPOL table switches to $35 \times 35$ for $\sigma > 30$.
• Filtering parameter $h$. It sets how quickly weights decay with patch distance and should be proportional to $\sigma$. Too small, and only near-identical patches are averaged, leaving noise; too large, and dissimilar patches are averaged, blurring detail. The IPOL table sets $h = k\sigma$, with $k$ from 0.4 down to 0.3 for grayscale images, decreasing for larger patches, whose distances between pure-noise patches concentrate more tightly around $2\sigma^2$. Values of $h$ depend on how the distance is normalized, so they do not transfer between implementations.

Complexity

Computed directly, each of $N$ pixels compares its patch with every patch in the search window, for a cost of $O(N (2r+1)^2 (2f+1)^2)$. Buades, Coll and Morel report that with $7 \times 7$ patches and a $21 \times 21$ window, about $49 \times 441$ operations per pixel, a $512 \times 512$ image took about 30 seconds on a PC of the time. They also gave a multiscale version that searches only near the matches found in a half-resolution image.

The patch distance can be computed in constant time per pixel and offset. For a fixed offset between $p$ and $q$, the box-weighted patch distance is a sum of squared differences between the image and its shifted copy over a window. Darbon and colleagues (2008) computed that sum for every pixel at once from an integral image of the squared-difference image, replacing the Gaussian patch weighting with a constant one. The cost becomes $O(N (2r+1)^2)$, independent of patch size, and the loop over offsets vectorizes and parallelizes well. The implementation below uses the same idea with a box filter.

The patchwise variant, also described in the IPOL article, denoises whole patches with the weights computed for their centers and averages the overlapping estimates of each pixel, which raises the PSNR and reduces spurious noise oscillations near edges, though the authors note it does not preserve more detail.

Implementation

The example builds a $128 \times 128$ test image with a low-contrast checkerboard of $4 \times 4$-pixel squares (values 0.4 and 0.6), a smooth ramp, and a bright disk, and adds Gaussian noise with $\sigma = 0.08$, about 20 on a 0–255 scale. NL-means uses the IPOL table’s settings for that noise level: $5 \times 5$ patches, a $21 \times 21$ search window, and $h = 0.4\sigma$. It is compared with Gaussian filters and with bilateral filters (spatial $\sigma_s = 1.5$, range $\sigma_r$ 0.1 or 0.2), reporting PSNR over the whole image and over the textured left half.

import numpy as np
from scipy import ndimage


def nl_means(v, sigma, f=2, r=10, k=0.4):
    """Pixelwise non-local means (grayscale), in the IPOL formulation.

    f: patch radius (patch is (2f+1)^2), r: search radius, h = k * sigma.
    """
    h = k * sigma
    pad = np.pad(v, r + f, mode="reflect")
    H, W = v.shape
    num = np.zeros_like(v)
    den = np.zeros_like(v)
    wmax = np.zeros_like(v)
    center = pad[r:r + H + 2 * f, r:r + W + 2 * f]
    for dy in range(-r, r + 1):
        for dx in range(-r, r + 1):
            if dy == 0 and dx == 0:
                continue
            shifted = pad[r + dy:r + dy + H + 2 * f, r + dx:r + dx + W + 2 * f]
            # Mean squared patch difference for this shift, for every pixel:
            # one box filter of the squared-difference image.
            d2 = ndimage.uniform_filter((center - shifted) ** 2, 2 * f + 1)
            d2 = d2[f:f + H, f:f + W]
            w = np.exp(-np.maximum(d2 - 2 * sigma**2, 0) / h**2)
            num += w * shifted[f:f + H, f:f + W]
            den += w
            wmax = np.maximum(wmax, w)
    num += wmax * v            # the pixel itself gets the largest weight seen
    den += wmax
    return num / den


def bilateral(img, sigma_s, sigma_r, truncate=3.0):
    r = int(np.ceil(truncate * sigma_s))
    pad = np.pad(img, r, mode="reflect")
    num = np.zeros_like(img)
    den = np.zeros_like(img)
    H, W = img.shape
    for dy in range(-r, r + 1):
        for dx in range(-r, r + 1):
            nb = pad[r + dy:r + dy + H, r + dx:r + dx + W]
            w = np.exp(-(dx**2 + dy**2) / (2 * sigma_s**2)
                       - (nb - img) ** 2 / (2 * sigma_r**2))
            num += w * nb
            den += w
    return num / den


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


# Synthetic 128 x 128 test image in [0, 1]: a low-contrast checkerboard of small
# squares on the left, a smooth ramp on the right, and a bright disk.
y, x = np.mgrid[0:128, 0:128]
clean = np.where(((x // 4) + (y // 4)) % 2 == 0, 0.4, 0.6)
clean[:, 64:] = 0.2 + 0.6 * (x[:, 64:] - 64) / 63
clean[(x - 96) ** 2 + (y - 32) ** 2 < 15**2] = 0.95

rng = np.random.default_rng(0)
sigma = 0.08                                # about 20 on a 0-255 scale
noisy = clean + rng.normal(0, sigma, clean.shape)

results = {"noisy input": noisy}
for s in [0.5, 0.7, 1.0]:
    results[f"Gaussian, sigma {s}"] = ndimage.gaussian_filter(noisy, s)
for sr in [0.1, 0.2]:
    results[f"bilateral, 1.5 / {sr}"] = bilateral(noisy, 1.5, sr)
results["non-local means"] = nl_means(noisy, sigma)

tex = (slice(None), slice(0, 64))           # textured half
for name, out in results.items():
    print(f"{name:22s} PSNR {psnr(out, clean):5.2f} dB   "
          f"texture half {psnr(out[tex], clean[tex]):5.2f} dB")

Output:

noisy input            PSNR 21.97 dB   texture half 21.95 dB
Gaussian, sigma 0.5    PSNR 25.25 dB   texture half 24.89 dB
Gaussian, sigma 0.7    PSNR 26.31 dB   texture half 24.94 dB
Gaussian, sigma 1.0    PSNR 25.39 dB   texture half 23.41 dB
bilateral, 1.5 / 0.1   PSNR 25.66 dB   texture half 24.24 dB
bilateral, 1.5 / 0.2   PSNR 26.47 dB   texture half 24.07 dB
non-local means        PSNR 29.15 dB   texture half 27.44 dB

The checkerboard contrast, 0.2, is only 2.5 times the noise. The Gaussian filter gains at most about 3 dB on the texture before it starts erasing the squares, and the bilateral filter does no better there: single-pixel differences in the texture are comparable to the noise, so the range kernel cannot tell them apart. NL-means compares $5 \times 5$ patches, finds the many exact repetitions of the checkerboard within the search window, and gains about 5.5 dB on the texture and 7 dB overall. The test image favors it, since a repetitive pattern is its best case.

Properties and Behavior

• A neighborhood filter on patches. NL-means generalizes the Yaroslavsky neighborhood filter, a bilateral filter with a box spatial kernel, by replacing the difference of two pixel values with the distance between two patches, which makes the weights far more robust to noise.
• Consistency. Buades, Coll and Morel proved that, if the image is modeled as a stationary random process, the NL-means estimate converges as the image grows to the conditional expectation of the clean pixel value given its noisy neighborhood.
• Method noise. The same authors judged denoisers by their method noise, the difference $v - \mathrm{NL}[v]$ between the input and its denoised version. An ideal method removes only noise, so its method noise should look like white noise and contain no edges or texture. In their visual comparison with Gaussian smoothing, anisotropic filtering, total-variation methods, neighborhood filters, wavelet thresholding, and a DCT empirical Wiener filter, they judged the method noise of NL-means the most like white noise.

Limitations

• Rare patches. A corner, a small isolated detail, or any unique structure has few similar patches within the search window. Its weights then concentrate on the pixel itself and the noise largely remains; with a large $h$, the structure is instead averaged with dissimilar patches and blurred. Buades, Coll and Morel proposed a correction that mixes back more of the noisy value where the variance of the averaged pixels is much larger than $\sigma^2$.
• Cost. Even with constant-time patch distances, a $21 \times 21$ window means hundreds of operations per pixel.
• Noise model. The weights assume white noise with a known, uniform $\sigma$. Signal-dependent or spatially correlated noise makes distances unreliable unless the noise is first transformed or modeled.

Variants

• BM3D. Dabov and colleagues (2007) group the most similar patches into a 3D stack and filter it jointly by shrinkage in a 3D transform domain, in two stages: hard thresholding, then Wiener filtering guided by the first estimate. It became a standard baseline for later methods.
• Low-rank and learned models. Later methods treat groups of similar patches as approximately low rank, as in weighted nuclear norm minimization, or learn patch priors from data. Deep denoisers usually do not match patches explicitly, though attention layers compare features across the image in a related way.
• Color and video. For color images the patch distance is computed over all channels and the same weights average each channel; Buades, Coll and Morel also applied the method to films.