Image Denoising: the Human–Machine Competition

Contents

• Introduction and overview
• 0. Review: probability, linear algebra and calculus
• 1. Image sampling and the discrete Fourier transform
• 2. A fast tutorial on Gaussian vectors
• 3. Noise
• 4. Multiscale DCT denoising
• 5. Image self-similarity and denoising
• 6. Bayesian patch-based methods
• 7. A crash course on deep learning
• 8. Deep and shallow CNN denoisers: a comparative study
• 9. Comparison of all denoising methods
• 10. Research directions and open problems
• Exercises
• Projects and assessment
• References

Introduction and overview

Audience: a student with a good M1 in mathematics, computer science or physics who knows elementary probability and Fourier analysis, has perhaps never trained a neural network, and wants to understand why a modern camera can photograph a dark room.

Empirical first: add Gaussian noise of standard deviation 20 to a 256 × 256 grey crop of the china.jpg image shipped with scikit-learn, and the PSNR falls to 22.1 dB; the best Gaussian blur recovers only 24.6 dB, while sliding-window DCT thresholding reaches 26.8 dB, non-local means 27.0 dB and a two-step NL-Bayes 27.8 dB (code in ch. 4–6). Why do these principles rank this way, and what does a trained network add?

The course studies one of the oldest problems of signal processing: separating a signal from its noise. It follows the history of the solutions — Fourier and Gaussian models, transform thresholding, self-similarity, Bayesian patch models — and then trains neural denoisers and compares both families on the same images. The common object is the image patch. Related MVA pages: Sub-Pixel Image Processing shares the sampling theory, Mathematical Models for Digital Photography the camera noise, and Deep Learning the networks.

What you will be able to do

• Explain how an image is sampled, what aliasing is, and compute with the 2D discrete Fourier and cosine transforms.
• Manipulate Gaussian vectors: affine images, conditioning, the norm of white noise, Stein’s lemma.
• Model the noise of a camera (Gaussian, Poisson, signal dependent), stabilise its variance and estimate its level from a single image.
• Derive the Wiener filter and the MMSE estimator, and state when a linear filter is optimal.
• Implement and analyse thresholding in an orthonormal basis, sliding-window DCT denoising and its multiscale version.
• Implement non-local means and explain the statistics of patch distances and the structure of the space of patches.
• Derive NL-Bayes and BM3D as empirical Wiener filters on groups of similar patches.
• Build, train and evaluate a convolutional denoiser in PyTorch, knowing why the L2 loss learns the conditional expectation.
• Prove that a sliding-window transform-thresholding denoiser is a shallow convolutional network, and read a bias-free network as an adaptive linear filter.
• Design a fair comparison of denoisers: metrics, paired statistics, noise mismatch, real noise.

Official objectives and syllabus (2026–2027)

The course description on master-mva.com states three objectives: explore the structure of images at the level of patches — small image extracts, processed both by neural networks and by recent image processing, whose currently understood size ranges from about 8 × 8 = 64 to 60 × 60 = 3600 pixels; apply this structure to a fundamental problem that requires understanding it, image denoising; and compare image processing designed by humans with image processing learned by deep neural networks.

The description recalls that the separation of noise and signal was already the key problem of Shannon’s Mathematical Theory of Communication, that denoising approaches identified since the 1970s group into a handful of “denoising principles”, and that since 2016 neural denoisers have started to outperform, slightly, the human-made algorithms, with inner workings not yet fully understood. Each session combines a course and online experiments on IPOL.

Every official topic is taught on this page. Its extensions beyond the syllabus are the review chapter 0, Stein’s unbiased risk estimate (ch. 2 and 4), Gaussian-mixture patch priors (ch. 6.4), and the research chapter 10 on self-supervised denoising, denoisers as priors and video.

Where this course fits in the programme

Official prerequisites: elementary probability, Fourier analysis and differential calculus; deep learning is introduced from scratch. It is a Modelling course of semester 1, with 5 ECTS and 30 h (20 h of lectures and 10 h of tutorials).

• Preparation — the refreshers Probability and Statistics (Gaussian vectors, conditional expectation), Linear Algebra, Analysis and Optimization (spectral theorem, gradients) and Scientific Python, NumPy and PyTorch (the code of every chapter).
• Related courses of semester 1 — Sub-Pixel Image Processing (sampling and interpolation), Mathematical Models for Digital Photography (camera pipeline and noise), The Machine Intelligence of Images, Deep Learning, Computational Statistics and Convex Optimization.
• Courses of semester 2 that build on it — Inverse Problems and Imaging, Deep Learning for Image Restoration and Synthesis, Generative Models for Images, Diffusion and Multiscale Models in AI (denoisers are the engine of diffusion models), Remote Sensing Data and Deep Learning for Medical Imaging.

Course outline

• Review: probability, linear algebra and calculus — covariance, the bias–variance split of the error, the conditional expectation as best estimator, orthonormal bases.
• Image sampling and the discrete Fourier transform — the Shannon–Whittaker theorem, aliasing, the 2D DFT and DCT, the spectrum of natural images.
• A fast tutorial on Gaussian vectors — definition, independence, invariance of white noise, conditioning and Stein’s lemma.
• Noise — Shannon’s channel, camera noise, variance stabilisation, noise estimation, the Wiener filter and the denoising principles.
• Multiscale DCT denoising — thresholding in an orthonormal basis, its oracle risk, sliding-window DCT and the multiscale fix of low-frequency noise.
• Image self-similarity and denoising — neighbourhood filters, non-local means, method noise and the space of patches.
• Bayesian patch-based methods — Gaussian patch models, NL-Bayes, BM3D, Gaussian mixtures and EPLL.
• A crash course on deep learning — networks, the L2 loss and the MMSE estimator, backpropagation, SGD, convolutions, residual learning.
• Deep and shallow CNN denoisers — MLP, DnCNN and FFDNet, the shallow network hidden in DCT denoising, bias-free networks.
• Comparison of all denoising methods — protocols, inherent bounds and what each family teaches the other.
• Research directions and open problems — self-supervision, denoisers as priors and diffusion, video and real noise (self-study).

Session plan

30 h in class = 10 sessions of 3 h. The volume is the Université Paris-Saclay programme table (20 h of lectures and 10 h of tutorials); master-mva.com announces 10 courses ending in the first week of December, each combining a course and online IPOL experiments. The plan follows the official order, one plan item per session (the crash course on deep learning takes two); the split of each session into 2 h of lecture and 1 h of exercises is a proposal.

Level 1 = foundations (reviews, definitions, first worked examples), 2 = core results and algorithms, 3 = advanced comparisons and recent research; the level never decreases.

Self-study, outside the 30 h (about 4 h in total): ch. 10 Research directions and open problems — self-supervised denoising, Tweedie’s formula and diffusion, plug-and-play priors, video and real noise (≈ 4 h, with Ex. 15).

The official assessment includes an intermediate written exam; on this plan it fits naturally after session 6, when the classical half of the course is complete.

How to study this module

• Read the chapters in order: each chapter uses only the earlier ones, and ch. 2 (Gaussian vectors) is used everywhere afterwards.
• After each session, do the exercises listed in its row, then run the code of the chapter and change one parameter (noise level, patch size, threshold) to see its effect.
• Reproduce the IPOL experiments named in the “IPOL experiment” notes on www.ipol.im: the official weekly report is about such online experiments and needs no programming.
• The Python blocks use NumPy, SciPy, scikit-learn (for its two bundled sample images, china.jpg and flower.jpg) and PyTorch on CPU; every block runs in less than half a minute on a laptop.
• Keep a table of PSNR values for the same noisy image across chapters: it becomes the comparison of ch. 9.

0. Review: probability, linear algebra and calculus

Every chapter of the course measures an error, averages noisy values or changes basis. This review fixes the notation and proves the three facts used everywhere: covariances transform simply, the error of an estimator splits into bias and variance, and the conditional expectation is the best possible estimator.

The following Monte Carlo experiment checks Example 0.4.

import numpy as np

rng = np.random.default_rng(0)
x = np.r_[np.zeros(10), 60 * np.ones(10)]                # a step of height 60
sigma, trials = 20.0, 20000
Y = x + sigma * rng.standard_normal((trials, x.size))  # noisy realisations
Xh = (np.roll(Y, 1, 1) + Y + np.roll(Y, -1, 1)) / 3    # 3-pixel moving average
bias2 = (Xh.mean(0) - x) ** 2
var = Xh.var(0)
print("squared bias, samples 8-11:", np.round(bias2[8:12], 1))
print("variance,     samples 8-11:", np.round(var[8:12], 1))

The output shows squared biases close to 0, 400, 400, 0 and variances close to 133 on the four samples around the step.

1. Image sampling and the discrete Fourier transform

A digital image is a finite array of samples of a continuous light field. This chapter says what sampling preserves and what it destroys (aliasing), and builds the two orthonormal bases — the discrete Fourier and cosine transforms — in which signal and noise will be separated.

1.1 Band-limited images and the sampling theorem

We use the Fourier transform \hat u(\xi) = \int u(x)\, e^{-i x \xi}\, dx on \mathbb{R}; the 2D statements follow by tensor products. A function is band-limited to [-\pi, \pi] when \hat u vanishes outside this interval. Sampling means keeping the values u(k), k \in \mathbb{Z}.

When \hat u does not vanish outside [-\pi, \pi], the sampled spectrum mixes several frequencies: this is aliasing. In a camera, the blur of the optics and of the sensor area acts as the low-pass filter that limits it.

1.2 The 2D discrete Fourier transform

The DFT treats the image as periodic. The jump between opposite borders creates artificial high frequencies — the bright cross visible in the spectrum of most photographs. Symmetrising the image before the transform removes this jump; this is the cosine transform of Section 1.4.

1.3 Trigonometric interpolation and zoom

For an N-periodic discrete signal with N odd, the DFT interpolate is the trigonometric polynomial

\bar u(x) = \frac{1}{N}\sum_{|p| \le (N-1)/2} \hat u(p)\, e^{2i\pi p x / N}.

It is the only trigonometric polynomial of degree at most (N-1)/2 that matches the samples. It gives exact sub-pixel translations — multiply \hat u(p) by e^{-2i\pi p \tau / N} — and zooms by zero-padding the spectrum. Its price is ringing (Gibbs oscillations) near sharp edges.

1.4 The discrete cosine transform

The DCT is the DFT of the image reflected across its borders. There is no artificial jump, so the energy of smooth patches concentrates in few low-frequency coefficients. This energy compaction explains its use in JPEG and in the denoisers of ch. 4.

1.5 The spectrum of natural images

Measured on large collections, the power spectrum of natural images decays roughly like |\xi|^{-2} (Field, 1987; Ruderman and Bialek, 1994). White noise has a flat spectrum, equal to \sigma^2 for the normalised DFT. The signal-to-noise ratio per frequency therefore decreases with |\xi|.

The code below builds the 8 \times 8 DCT, checks its orthonormality, and compares the energy compaction of image patches and of noise patches.

import numpy as np
from sklearn.datasets import load_sample_image

n = 8
j = np.arange(n)
D = np.sqrt(2 / n) * np.cos(np.pi * np.outer(j, j + 0.5) / n)   # rows = DCT-II atoms
D[0] /= np.sqrt(2)
print("orthonormal:", np.allclose(D @ D.T, np.eye(n)))

u = load_sample_image("china.jpg").mean(axis=2)
H, W = (s - s % n for s in u.shape)
P = u[:H, :W].reshape(H // n, n, W // n, n).transpose(0, 2, 1, 3).reshape(-1, n, n)
rng = np.random.default_rng(0)
N = 20 * rng.standard_normal(P.shape)                         # white noise, sigma = 20

def energy_in_top6(patches):
    """Mean fraction of the energy carried by the 6 largest of the 64 DCT coefficients."""
    C = np.einsum("ki,pij,lj->pkl", D, patches, D).reshape(len(patches), -1) ** 2
    C.sort(axis=1)
    return (C[:, -6:].sum(1) / C.sum(1)).mean()

print(f"image patches: {energy_in_top6(P - P.mean((1, 2), keepdims=True)):.2f}")
print(f"noise patches: {energy_in_top6(N):.2f}")
print(f"noise variance per DCT coefficient: {np.einsum('ki,pij,lj->pkl', D, N, D).var():.1f}")

With their mean removed, image patches put 69 % of their energy in their 6 largest coefficients, against 41 % for noise. The noise variance per coefficient stays 400: an orthonormal transform keeps white noise white (Theorem 2.4).

2. A fast tutorial on Gaussian vectors

White Gaussian noise is the reference model of the course, and Gaussian models of patches are the core of ch. 6. This chapter collects what we need: the definition, independence, the invariance of white noise under orthonormal transforms, the size of a noise vector, conditioning, and Stein’s lemma.

2.1 Definition and first properties

2.2 White noise in an orthonormal basis

This is why every transform method of ch. 4–6 — DCT, wavelets, the PCA of a group of patches — can treat its coefficients one by one with the same noise level \sigma. For the normalised DFT of real white noise, every coefficient \hat n(p) satisfies \mathbb{E}|\hat n(p)|^2 = \sigma^2, with the conjugate symmetry \hat n(-p) = \overline{\hat n(p)}.

2.3 The norm of a noise vector

2.4 Conditioning: the Gaussian denoiser

The estimator keeps the directions where the signal variance dominates the noise and kills the others. All the patch methods of ch. 6 apply this rule with estimated \lambda_i.

import numpy as np

rng = np.random.default_rng(0)
lam, sigma2, n = np.array([100.0, 25.0, 1.0]), 25.0, 200_000
U, _ = np.linalg.qr(rng.standard_normal((3, 3)))          # a random orthonormal basis
C = U @ np.diag(lam) @ U.T                                # prior covariance of X
X = rng.multivariate_normal(np.zeros(3), C, size=n)
Y = X + np.sqrt(sigma2) * rng.standard_normal((n, 3))
K = C @ np.linalg.inv(C + sigma2 * np.eye(3))             # Wiener gain (Theorem 2.7)
print("shrinkage factors:", np.round(np.linalg.eigvalsh(K), 4))
print(f"MSE of Y:        {((Y - X) ** 2).sum(1).mean():.2f}  (theory 75)")
print(f"MSE of K Y:      {((Y @ K.T - X) ** 2).sum(1).mean():.2f}  (theory {np.sum(lam * sigma2 / (lam + sigma2)):.2f})")

The simulation returns the factors 0.0385, 0.5, 0.8 and errors of 74.87 and 33.53, against 75 and 33.46 in theory.

2.5 Stein’s lemma and an unbiased estimate of the risk

The right-hand side involves only Y: it estimates the error of a denoiser without knowing the clean image. Chapter 4 uses it to choose a threshold, and ch. 8 reads \operatorname{div} f as the trace of a network Jacobian.

3. Noise

Before removing noise we must say what it is. This chapter starts from Shannon’s formulation of the problem, follows the noise through a camera, shows how to make it Gaussian and how to measure its level, and ends with the optimal linear denoiser and the map of the classical denoising principles.

3.1 Signal plus noise: Shannon’s problem

Shannon’s Mathematical Theory of Communication (1948) treats a message corrupted by a noisy channel. Its Gaussian case gives the amount of information that noise leaves in each sample.

For an image whose grey levels have standard deviation 50, observed with \sigma = 20, the bound is \frac12 \log_2(1 + 6.25) = 1.43 bit per pixel. A denoiser cannot create this missing information. It can only exploit the redundancy of images — smoothness, sparsity, self-similarity — exactly as Shannon’s codes exploit redundancy to correct errors. Each denoising principle of Section 3.6 is a hypothesis on this redundancy.

3.2 Noise in a camera

A pixel counts photons during the exposure. Three sources dominate.

• Photon (shot) noise. The count follows a Poisson law \mathcal{P}(\lambda), with mean and variance \lambda. Its signal-to-noise ratio is \sqrt{\lambda}: 10 for 100 photons, 100 for 10 000. A dark scene collects few photons and is noisy.
• Read-out and thermal noise. The electronics add a signal-independent, nearly Gaussian noise.
• Quantisation. Rounding to a step \Delta adds an error that is close to uniform.

The camera pipeline then demosaics, white-balances, applies a tone curve and compresses. After these steps the noise is correlated between neighbouring pixels and colour channels, and its variance depends on intensity and frequency. Removing such noise blindly, from a single JPEG image, is the problem addressed by the Noise Clinic of Lebrun, Colom and Morel (2015).

3.3 Variance stabilisation

The same computation for the affine noise curve gives the generalised Anscombe transform f(z) = \frac{2}{a}\sqrt{a z + \frac38 a^2 + b}. The usual pipeline stabilises, denoises as if the noise were Gaussian with \sigma = 1, and inverts; at low counts, the inverse must be unbiased rather than algebraic (Mäkitalo and Foi, 2011).

3.4 Estimating the noise level from one image

Take for d the finest diagonal Haar coefficient,

d(k,l) = \tfrac12\big(\tilde u(k,l) - \tilde u(k+1,l) - \tilde u(k,l+1) + \tilde u(k+1,l+1)\big).

Its noise part is \mathcal{N}(0, \sigma^2) by Proposition 0.2, and its signal part vanishes on smooth regions. Edges and textures occupy a minority of pixels, which the median ignores — up to a point.

More accurate methods select the blocks of lowest variance and measure their high-frequency DCT energy, intensity by intensity, which yields a whole noise curve (Colom and Buades, 2013).

import numpy as np
from sklearn.datasets import load_sample_image

rng = np.random.default_rng(0)
# 1) Anscombe: the variance of 2 sqrt(k + 3/8) for Poisson counts k
for lam in [1, 4, 10, 50, 200]:
    k = rng.poisson(lam, 200_000)
    print(f"lambda = {lam:3d}: Var(k) = {k.var():7.2f}, Var(Anscombe) = {(2 * np.sqrt(k + 3 / 8)).var():.3f}")

# 2) MAD estimate of sigma from the finest diagonal Haar coefficients
for name in ["flower.jpg", "china.jpg"]:
    u = load_sample_image(name).mean(axis=2)
    for sigma in [5, 20, 40]:
        f = u + sigma * rng.standard_normal(u.shape)
        d = (f[:-1, :-1] - f[1:, :-1] - f[:-1, 1:] + f[1:, 1:]) / 2   # noise part ~ N(0, sigma^2)
        print(f"{name:10s} sigma = {sigma:2d}: estimate {np.median(np.abs(d)) / 0.6745:5.2f}")

The Anscombe variance is 0.72 for \lambda = 1 and within 0.4 % of 1 from \lambda = 4 on. On the smooth flower.jpg the MAD estimates are 5.39, 20.23 and 40.23. On the textured china.jpg they are 8.32, 23.26 and 42.47: texture leaks into the high-pass coefficients and biases the estimate upwards, most visibly at low noise.

3.5 Linear filtering, method noise and the Wiener filter

A linear filter w applied to white noise gives a noise of variance \sigma^2 \|w\|_2^2 per pixel (Proposition 0.2). A k \times k box filter divides the variance by k^2. What it removes from the image is measured by the method noise.

The Gaussian filter removes the Laplacian of the image: edges, corners and textures. No choice of h escapes the bias–variance dilemma of Example 0.4.

The Wiener filter is optimal among all estimators when the image is Gaussian and stationary — which a photograph is not. One spectrum cannot describe both a flat sky and a textured roof: the global filter blurs the textures or leaves noise in the sky. Restricting the Gaussian model to small regions, the patches, is the step taken in ch. 4–6.

3.6 A handful of denoising principles

The course description groups the methods found since the 1970s into a few principles. Each one is a hypothesis on the redundancy of images.

The survey by Lebrun, Colom, Buades and Morel (2012) analyses the classical methods as combinations of a few ingredients shared across principles. Three of them return in every later chapter: the aggregation of overlapping estimates, a second pass guided by a first oracle estimate, and the choice of a patch model.

4. Multiscale DCT denoising

If the clean image is sparse in an orthonormal basis — few large coefficients — and the noise is spread evenly over all coefficients (Theorem 2.4), then keeping the large coefficients and killing the small ones separates the two. This chapter makes the idea quantitative, applies it to sliding DCT patches, and repairs its main defect with a multiscale scheme.

4.1 Thresholding in an orthonormal basis

In an orthonormal basis, the observation becomes y_k = x_k + n_k, k = 1, \dots, d, with independent n_k \sim \mathcal{N}(0, \sigma^2) and unknown coefficients x_k of the clean image.

The oracle risk counts the coefficients above the noise level. It is small exactly when the image is sparse in the basis: the choice of the basis is the hypothesis on the image.

With \tau_d, all pure-noise coefficients die with high probability, but slowly: for d = 64, \tau_d = 2.88\sigma and the bound is 0.28.

Up to a factor 2\log d, a simple rule matches the keep-or-kill oracle, and Donoho and Johnstone show that this factor cannot be improved asymptotically, uniformly over all signals. In practice the universal threshold is too high for the mean square error, and \tau is tuned on the data.

Minimising SURE over \tau chooses the threshold from the noisy data alone (the SureShrink rule of Donoho and Johnstone, 1995). The code applies it to the global DCT of a 256 \times 256 crop of flower.jpg with \sigma = 20.

import numpy as np
from scipy.fft import dctn
from sklearn.datasets import load_sample_image

rng = np.random.default_rng(0)
u = load_sample_image("flower.jpg").mean(axis=2)[:256, :256]
sigma = 20.0
x = dctn(u, norm="ortho").ravel()                        # clean coefficients
y = x + sigma * rng.standard_normal(x.size)              # noisy coefficients (Theorem 2.4)
d = y.size

def sure(t):    # Proposition 4.5, divided by d
    return (np.minimum(y**2, t**2).sum() - d * sigma**2 + 2 * sigma**2 * (np.abs(y) > t).sum()) / d

def true_mse(t):
    return ((np.sign(y) * np.maximum(np.abs(y) - t, 0) - x) ** 2).mean()

grid = sigma * np.array([0.5, 1.0, 1.5, 2.0, 2.5, 3.0, np.sqrt(2 * np.log(d))])
for t in grid:
    print(f"tau = {t / sigma:4.2f} sigma: SURE {sure(t):6.1f}   true MSE {true_mse(t):6.1f}")

SURE follows the true error within 3 units on the whole grid. Both are minimal at \tau = 1.5\sigma (MSE 70.2, i.e. 29.7 dB), while the universal threshold 4.71\sigma gives 131.6.

4.2 Sliding-window DCT denoising

A global basis adapts poorly to an image whose content changes from place to place. The DCT of small patches adapts to local content, and processing every patch position removes the dependence on a grid.

A pure-noise coefficient survives the threshold 3\sigma with probability 0.0027, about 0.17 coefficient per 8 \times 8 patch.

When the errors of overlapping patches are weakly correlated, the gain approaches a division by K. Averaging over all translations is the cycle spinning of Coifman and Donoho (1995).

4.3 Why one scale is not enough

The same holds for the few lowest frequencies, which survive in most patches. After aggregation these errors are smooth over about s pixels: faint blotches, visible in flat regions although small in MSE. For \sigma = 40 and 8 \times 8 patches, the DC error of each patch estimate has standard deviation 5 grey levels. Aggregation averages strongly correlated errors and only brings it down to about 3.3.

4.4 The multiscale DCT

At the coarse scale, the slow variations of the image become fast, and the patches of the same size see them at a lower noise level. This motivates the following scheme, in the spirit of the Multi-Scale DCT Denoising of Pierazzo, Morel and Facciolo (2017); the code is a simplified version.

The fraction r trades the cleaner low frequencies of the coarse scale against the band-limiting artefacts of the half-size image; the code uses r = 1/4.

import numpy as np
from scipy.fft import dctn, idctn
from sklearn.datasets import load_sample_image

def dct_denoise(f, sigma, s=8, thr=3.0):
    """Sliding-window DCT: hard thresholding at thr * sigma, uniform aggregation."""
    P = np.lib.stride_tricks.sliding_window_view(f, (s, s))   # every s x s patch
    C = dctn(P, axes=(-2, -1), norm="ortho")
    C[np.abs(C) < thr * sigma] = 0.0
    Q = idctn(C, axes=(-2, -1), norm="ortho")
    out, cnt = np.zeros_like(f), np.zeros_like(f)
    h, w = Q.shape[:2]
    for i in range(s):                                       # put every estimate back
        for j in range(s):
            out[i:i + h, j:j + w] += Q[:, :, i, j]
            cnt[i:i + h, j:j + w] += 1
    return out / cnt

def two_scale(f, sigma, r=0.25):
    """Coarse scale = low-frequency quarter of the global DCT, with noise sigma / 2."""
    F = dctn(f, norm="ortho")
    n = f.shape[0] // 2
    coarse = idctn(F[:n, :n], norm="ortho") / 2
    c_hat = dct_denoise(coarse, sigma / 2)
    G = dctn(dct_denoise(f, sigma), norm="ortho")
    k = int(r * n)
    G[:k, :k] = dctn(2 * c_hat, norm="ortho")[:k, :k]        # replace the lowest frequencies
    return idctn(G, norm="ortho")

rng = np.random.default_rng(0)
crops = {"china.jpg": (slice(100, 356), slice(200, 456)), "flower.jpg": (slice(0, 256), slice(0, 256))}
for name, sigma in [("china.jpg", 20.0), ("china.jpg", 40.0), ("flower.jpg", 20.0), ("flower.jpg", 40.0)]:
    u = load_sample_image(name).mean(axis=2)[crops[name]]
    f = u + sigma * rng.standard_normal(u.shape)
    psnr = lambda a: 10 * np.log10(255**2 / np.mean((a - u) ** 2))
    print(f"{name:10s} sigma {sigma:2.0f}: noisy {psnr(f):5.2f} | DCT {psnr(dct_denoise(f, sigma)):5.2f}"
          f" | two scales {psnr(two_scale(f, sigma)):5.2f} dB")

The local DCT beats the global SURE-tuned thresholding by 6 dB on the flower crop (35.62 against 29.7 dB). The second scale gains 0.8 dB where it matters — a smooth image at high noise — and almost nothing on the textured crop, as Proposition 4.8 predicts.

5. Image self-similarity and denoising

Transform thresholding assumes that each patch is sparse in a fixed basis. Self-similarity makes a different bet: most patches of an image have near copies elsewhere in the same image, and averaging these copies removes the noise without blurring. This chapter builds non-local means on this bet and studies the space of patches that makes it work.

5.1 From local to non-local averaging

The comparison of two single noisy grey levels is unreliable. The expected value of |\tilde u(i) - \tilde u(j)|^2 is |u(i) - u(j)|^2 + 2\sigma^2, with fluctuations of the same order as 2\sigma^2. At \sigma = 20, two pixels that differ by 20 grey levels cannot be told apart. Comparing whole patches fixes this, thanks to the concentration of Proposition 2.5.

5.2 Non-local means

The subtraction of 2\sigma^2 and the parameter h come from the statistics of the patch distance.

The identity above — a noisy value is an unbiased target for a predictor that does not see it — reappears in the self-supervised denoisers of ch. 10.

The code implements Definition 5.2 with 5 \times 5 patches, a 21 \times 21 window and h = 0.7\sigma, on the crop of ch. 4.

import numpy as np
from scipy.ndimage import uniform_filter, gaussian_filter
from sklearn.datasets import load_sample_image

rng = np.random.default_rng(0)
u = load_sample_image("china.jpg").mean(axis=2)[100:356, 200:456]
sigma = 20.0
f = u + sigma * rng.standard_normal(u.shape)
psnr = lambda a: 10 * np.log10(255**2 / np.mean((a - u) ** 2))

def nl_means(f, sigma, half_patch=2, half_win=10, h=0.7):
    """Pixelwise NL-means (Definition 5.2), vectorised over the offsets of the window."""
    r = half_win
    g = np.pad(f, r, mode="reflect")
    num, den = np.zeros_like(f), np.zeros_like(f)
    for dx in range(-r, r + 1):
        for dy in range(-r, r + 1):
            shifted = g[r + dx: r + dx + f.shape[0], r + dy: r + dy + f.shape[1]]
            # mean squared distance between the patches centred at i and at i + (dx, dy)
            d2 = uniform_filter((f - shifted) ** 2, size=2 * half_patch + 1, mode="reflect")
            w = np.exp(-np.maximum(d2 - 2 * sigma**2, 0.0) / (h * sigma) ** 2)
            num += w * shifted
            den += w
    return num / den

best = max((psnr(gaussian_filter(f, s)), s) for s in [0.3, 0.4, 0.5, 0.6, 0.7, 1.0])
print(f"noisy {psnr(f):.2f} dB | best Gaussian {best[0]:.2f} dB (h = {best[1]})"
      f" | NL-means 5x5 / 21x21 {psnr(nl_means(f, sigma)):.2f} dB")

NL-means reaches 27.01 dB against 24.59 dB for the best Gaussian filter and 26.78 dB for the sliding DCT of ch. 4. The weights adapt to the image: on a straight edge, the similar patches lie along the edge, and the average does not cross it.

5.3 Method noise and parameters

The method noise of NL-means (Definition 3.8) is close to white noise on textured and periodic images, where many exact copies exist. Its parameters act as follows.

• h too small: only the central patch has a large weight, and noise remains.
• h too large: dissimilar patches enter the average, and details blur.
• Window size: a larger window finds more copies but also more false matches, and costs more.
• Rare patches — corners, ends of lines, unique details — have few copies. They are barely denoised: this is the rare patch effect, and the main residual defect of the method.

5.4 The space of patches

The objectives of the course single out patches from 8 \times 8 to 60 \times 60 pixels. Two facts govern their use: larger patches are compared more reliably, and larger patches have fewer copies.

For \sigma = 20 and k = 3, an 8 \times 8 patch separates differences of root mean square \bar\Delta \approx 20.6 grey levels, and a 60 \times 60 patch separates \bar\Delta \approx 7.5.

But a 60 \times 60 region rarely has a copy in the same image. Every patch method chooses a size between these two constraints, and deep networks face the same trade-off with their receptive field (ch. 7.5).

The global ensemble of patches is far from filling its space. Natural 3 \times 3 patches of high contrast concentrate near a low-dimensional manifold of edges (Lee, Pedersen and Mumford, 2003). Topological data analysis shows that a dense core of them has the topology of a Klein bottle (Carlsson, Ishkhanov, de Silva and Zomorodian, 2008). Small patches recur within the same image, at the same and at other scales (Zontak and Irani, 2011).

The PCA of all the patches of an image gives a first, linear, picture of this structure. The number of principal variances above \sigma^2 is the number of directions that a Wiener filter keeps (Corollary 2.8).

import numpy as np
from sklearn.datasets import load_sample_image

u = load_sample_image("china.jpg").mean(axis=2)
for s in [4, 8, 16]:
    P = np.lib.stride_tricks.sliding_window_view(u, (s, s))[::2, ::2].reshape(-1, s * s)
    lam = np.linalg.eigvalsh(np.cov(P, rowvar=False))[::-1]       # patch PCA spectrum
    above = [int((lam > sig**2).sum()) for sig in (5, 20, 40)]
    print(f"{s:2d}x{s:<2d} (d = {s*s:3d}): top eigenvalues {np.round(lam[:3]).astype(int)},"
          f" eigenvalues above sigma^2 for sigma = 5, 20, 40: {above}")

At \sigma = 40 a global Gaussian model keeps only 7 of the 64 directions of an 8 \times 8 patch, the same for every patch. Groups of similar patches have a much lower dimension and their own principal axes: ch. 6 builds one Gaussian model per group.

6. Bayesian patch-based methods

NL-means averages similar patches with scalar weights. A group of similar patches carries more information: its mean and its covariance. This chapter models each group as a Gaussian vector and applies the Wiener estimator of ch. 2 to it. The result, NL-Bayes, and its sibling BM3D are among the best classical denoisers.

6.1 A Gaussian model for groups of similar patches

Collect K noisy patches \tilde P_1, \dots, \tilde P_K similar to a reference patch, with d pixels each. Model the clean patches as independent draws of \mathcal{N}(m, C) and the noise as white.

The estimator is a PCA of the group followed by a Wiener shrinkage: directions where the group varies more than the noise are kept, the others are attenuated. Unlike the DCT, the basis is adapted to each group.

When K is not much larger than d, the sample eigenvalues are noisy. For pure noise, as d, K \to \infty with d/K \to \gamma < 1, the eigenvalues of \tilde C spread over the interval of Marchenko and Pastur (1967):

\big[\sigma^2 (1 - \sqrt\gamma)^2,\ \sigma^2 (1 + \sqrt\gamma)^2\big].

6.2 NL-Bayes

The full algorithm adds a test for flat groups — if the group variance is close to \sigma^2, every patch is replaced by the mean — takes the step-2 group mean from the basic estimate, and works on colour patches. The code below keeps the two steps.

import numpy as np
from sklearn.datasets import load_sample_image

rng = np.random.default_rng(0)
u = load_sample_image("china.jpg").mean(axis=2)[100:356, 200:456]
sigma = 20.0
f = u + sigma * rng.standard_normal(u.shape)
psnr = lambda a: 10 * np.log10(255**2 / np.mean((a - u) ** 2))

def nl_bayes_step(f, guide, sigma, s=5, K=60, win=10, stride=4, oracle=False):
    """One step of a simplified NL-Bayes: group, Gaussian (Wiener) estimate, aggregate."""
    H, W = f.shape
    view = lambda a: np.lib.stride_tricks.sliding_window_view(a, (s, s)).reshape(H - s + 1, W - s + 1, -1)
    Pf, Pg = view(f), view(guide)
    off = (np.arange(s)[:, None] * W + np.arange(s)[None, :]).ravel()   # pixel offsets in a patch
    out, cnt = np.zeros(H * W), np.zeros(H * W)
    rows = sorted(set(range(0, H - s + 1, stride)) | {H - s})
    cols = sorted(set(range(0, W - s + 1, stride)) | {W - s})
    for i in rows:
        for j in cols:
            i0, i1 = max(i - win, 0), min(i + win, H - s) + 1
            j0, j1 = max(j - win, 0), min(j + win, W - s) + 1
            d2 = ((Pg[i0:i1, j0:j1] - Pg[i, j]) ** 2).sum(-1).ravel()
            idx = np.argpartition(d2, K)[:K]                  # the K most similar patches
            ii, jj = np.unravel_index(idx, (i1 - i0, j1 - j0))
            ii, jj = ii + i0, jj + j0
            Y, m = Pf[ii, jj], Pf[ii, jj].mean(0)               # noisy group (K x d) and its mean
            if oracle:   # step 2: covariance of the basic estimate
                C = np.cov(Pg[ii, jj], rowvar=False)
            else:        # step 1: Cov(noisy) - sigma^2 I, eigenvalues clipped at 0
                lam, U = np.linalg.eigh(np.cov(Y, rowvar=False))
                C = (U * np.maximum(lam - sigma**2, 0)) @ U.T
            Wf = C @ np.linalg.inv(C + sigma**2 * np.eye(s * s))  # Wiener matrix
            est = m + (Y - m) @ Wf.T
            pix = ((ii * W + jj)[:, None] + off[None, :]).ravel()
            np.add.at(out, pix, est.ravel())                  # aggregate every estimated patch
            np.add.at(cnt, pix, 1.0)
    return (out / cnt).reshape(H, W)

basic = nl_bayes_step(f, f, sigma)
final = nl_bayes_step(f, basic, sigma, oracle=True)
print(f"noisy {psnr(f):.2f} dB | step 1 {psnr(basic):.2f} dB | step 2 {psnr(final):.2f} dB")

The first step reaches 27.33 dB on the crop of ch. 4 (\sigma = 20), and the oracle step 27.75 dB, above NL-means (27.01 dB) and the sliding DCT (26.78 dB). The second step gains 0.4 dB because the covariance of the basic estimate is no longer inflated by noise eigenvalues (Example 6.3).

6.3 BM3D as a collaborative Wiener filter

Step 2 is the oracle attenuation of Proposition 4.2, with the unknown coefficients x_k replaced by the pilot estimate \hat b_k. NL-Bayes and BM3D share the same skeleton — group, estimate jointly, aggregate, repeat with an oracle. They differ in the basis: fixed for BM3D, the PCA of the group for NL-Bayes. The collaborative filtering of BM3D thus exploits both sparsity (ch. 4) and self-similarity (ch. 5).

6.4 Learned patch priors: Gaussian mixtures and EPLL

NL-Bayes learns one Gaussian per group from the noisy image itself: an internal prior. One can instead learn a prior on patches from a database of clean images: an external prior.

The MMSE estimate is a mixture of Wiener filters, each selected softly by how well its component explains the noisy patch. Keeping only the most likely component gives a piecewise linear estimator (Yu, Sapiro and Mallat, 2012).

The EPLL method of Zoran and Weiss (2011) learns a Gaussian mixture with about two hundred components on 8 \times 8 patches of natural images. It then computes a whole image \hat u by approximately solving

\min_u\ \frac{\lambda}{2}\|u - \tilde u\|^2 - \sum_i \log p(R_i u),

where R_i extracts the patch at i. The solver alternates a patch step (the Wiener filter of the most likely component) and an image step (a closed-form aggregation), while the weight of the patch term grows.

A Gaussian mixture learned on a database is a small statistical model trained on examples. The networks of ch. 7–8 take the same path with millions of parameters.

7. A crash course on deep learning

Since 2016, the best denoisers have been neural networks trained on examples. This chapter introduces them from scratch, as the official course does. Part I (ch. 7.1–7.4) covers networks, losses, gradients and optimisation; part II (ch. 7.5–7.7) covers convolutional networks and the training of a denoiser.

7.1 Neurons, layers and networks

Continuous functions on an interval are uniform limits of piecewise linear interpolants, so one hidden layer approximates any of them. This is the 1D case of the universal approximation theorem (Cybenko, 1989, for sigmoids; Leshno, Lin, Pinkus and Schocken, 1993, for every non-polynomial activation). The theorem says nothing about the size needed, nor about learning the weights from data.

7.2 Learning: the L2 loss estimates the MMSE denoiser

A trained network is therefore an approximation of the MMSE denoiser for the distribution of its training pairs: both the images and the noise model. This is the precise meaning of the course description’s remark that neural algorithms “learn jointly noise and image structure”. A change of either distribution at test time breaks the guarantee (ch. 8.5).

7.3 Backpropagation

Applied from the last layer to the first, these formulas compute all the gradients for about twice the cost of one forward pass. This is backpropagation (Rumelhart, Hinton and Williams, 1986); libraries such as PyTorch build it automatically.

7.4 Gradient descent and its stochastic version

The loss of a network is not convex, so this guarantees a small gradient, not a global minimum. In practice the gradient is computed on a random mini-batch of training pairs: an unbiased but noisy estimate (stochastic gradient descent, SGD). The Adam method (Kingma and Ba, 2015) rescales each coordinate by running averages of the gradient g_t and of its square:

m_t = \beta_1 m_{t-1} + (1 - \beta_1) g_t, \qquad v_t = \beta_2 v_{t-1} + (1 - \beta_2) g_t^2, \qquad \theta_{t+1} = \theta_t - \eta\, \frac{\hat m_t}{\sqrt{\hat v_t} + \epsilon},

where \hat m_t and \hat v_t are the averages corrected for their initialisation at zero.

7.5 Convolutional layers and receptive fields

A network with L = 17 layers sees 35 \times 35 pixels, inside the range of patch sizes named by the course objectives. A convolutional denoiser is thus a patch-based method: each output pixel is a non-linear function of one patch of the input, and the same function is applied everywhere. One 3 \times 3 layer from 64 to 64 channels has 64 \times (9 \times 64 + 1) = 36\,928 parameters.

7.6 Residual learning and normalisation

A residual denoiser predicts the noise rather than the image:

\hat x = y - R_\theta(y).

The identity, which already scores the PSNR of the input, is then the network R_\theta = 0, and training only learns a correction. Residual connections inside deep networks (He, Zhang, Ren and Sun, 2016) play the same role between layers and ease the optimisation of deep stacks.

Batch normalisation (Ioffe and Szegedy, 2015) standardises each channel over a mini-batch,

\hat z = \gamma\, \frac{z - \mu_B}{\sqrt{\sigma_B^2 + \epsilon}} + \beta,

with learned \gamma and \beta; at test time running averages replace \mu_B and \sigma_B^2. It speeds up training; ch. 8.4 shows that its additive constants matter for generalisation across noise levels.

7.7 Training a denoiser in practice

The standard recipe has five ingredients.

• Data: random crops of clean images, flipped and rotated, with synthetic noise drawn afresh at every step.
• Noise: a fixed \sigma (non-blind) or a range of levels (blind).
• Loss and optimiser: mean square error and Adam, with a decreasing learning rate.
• Validation: PSNR on held-out images, never on training crops.
• Budget: real denoisers train on hundreds of images for hours on a GPU.

The code trains a four-layer residual CNN with 19 105 parameters for 300 Adam steps on the parts of china.jpg outside the test crop, then tests it on the crop of ch. 4–6.

import numpy as np, torch, torch.nn as nn
from sklearn.datasets import load_sample_image

torch.manual_seed(0); torch.set_num_threads(1)     # small convolutions: one thread is fastest
rng = np.random.default_rng(0)
sigma = 20.0
img = load_sample_image("china.jpg").mean(axis=2)
test = img[100:356, 200:456] / 255                                 # the crop used since ch. 4
train = torch.tensor(np.concatenate([img[:, :200], img[:, 456:]], axis=1) / 255, dtype=torch.float32)

class SmallDnCNN(nn.Module):
    """Residual CNN x_hat = y - R(y): four 3x3 convolutions, receptive field 9 x 9."""
    def __init__(self, depth=4, ch=32):
        super().__init__()
        layers = [nn.Conv2d(1, ch, 3, padding=1), nn.ReLU()]
        for _ in range(depth - 2):
            layers += [nn.Conv2d(ch, ch, 3, padding=1), nn.ReLU()]
        layers += [nn.Conv2d(ch, 1, 3, padding=1)]
        self.R = nn.Sequential(*layers)

    def forward(self, y):
        return y - self.R(y)

net = SmallDnCNN()
opt = torch.optim.Adam(net.parameters(), lr=2e-3)
H, W = train.shape
for step in range(300):
    i = torch.randint(0, H - 32, (16,)); j = torch.randint(0, W - 32, (16,))
    x = torch.stack([train[a:a + 32, b:b + 32] for a, b in zip(i, j)])[:, None]
    y = x + sigma / 255 * torch.randn_like(x)                      # clean / noisy pairs
    loss = ((net(y) - x) ** 2).mean()                              # L2 loss: learns E[X | Y]
    opt.zero_grad(); loss.backward(); opt.step()

y_test = test + sigma / 255 * rng.standard_normal(test.shape)
with torch.no_grad():
    x_hat = net(torch.tensor(y_test, dtype=torch.float32)[None, None])[0, 0].numpy()
psnr = lambda a: 10 * np.log10(1 / np.mean((a - test) ** 2))
n_par = sum(p.numel() for p in net.parameters())
print(f"noisy {psnr(y_test):.2f} dB | CNN {psnr(x_hat):.2f} dB | {n_par} parameters")

After about twenty seconds of CPU training, the network reaches 27.00 dB: as good as NL-means (27.01 dB), better than the sliding DCT (26.78 dB), below NL-Bayes (27.75 dB). Trained the same way on flower.jpg instead — a smooth image with little texture — it reaches only about 25.4 dB on the same test crop. The network has learned the image statistics of its training set, as Theorem 7.4 predicts.

8. Deep and shallow CNN denoisers: a comparative study

This chapter compares learned denoisers with each other and with the classical ones. It shows that the classical transform denoisers are shallow convolutional networks with hand-chosen weights, and that a deep network without biases is an adaptive linear filter. These two facts make the human–machine comparison precise.

8.1 The first competitive neural denoiser: a plain MLP

Burger, Schuler and Harmeling (2012) trained a plain multilayer perceptron to map noisy patches to clean patches. The training set contained a very large number of patch pairs drawn from natural images. The network matched BM3D at several noise levels, which showed that a generic learned function could rival the best hand-designed principle — given enough data and computation. Their experiments also showed that larger input patches and larger training sets kept improving the results.

8.2 DnCNN and FFDNet

With 17 layers, the receptive field is 35 \times 35 pixels (Proposition 7.9). DnCNN was among the first deep networks to clearly outperform BM3D on Gaussian noise, by a margin of a fraction of a dB — the “slight” advantage mentioned in the course description.

FFDNet (Zhang, Zuo and Zhang, 2018) changes two things. It takes a noise-level map as an extra input, so one network handles a range of levels and even spatially varying noise. It also rearranges the image into four half-resolution sub-images before the convolutions, which enlarges the receptive field and cuts the cost. This is a multiscale idea, as in ch. 4.4.

8.3 Shallow networks hidden in the classical algorithms

Soft thresholding fits the same frame with ReLUs, since

\eta_S(t) = \rho(t - \tau) - \rho(-t - \tau) .

A sliding-DCT soft-thresholding denoiser is therefore a ReLU network with 2 s^2 hidden channels: filters \pm c_k, all with bias -\tau, and fixed weights.

The comparison becomes concrete. The human algorithm fixes the filters (DCT atoms), the non-linearity (a threshold at 3\sigma) and the depth (two layers). Learning these filters and shrinkage functions, layer after layer, gives the shrinkage fields of Schmidt and Roth (2014) and the trainable reaction–diffusion of Chen and Pock (2017): shallow learned networks. Deep CNNs replace the explicit transform by twenty learned layers.

The code below measures the step from human to learned design. It starts from the soft-threshold version of Proposition 8.2, with 8 \times 8 DCT filters and a threshold 1.5\sigma, then trains its filters and its 64 thresholds with the data and the budget of ch. 7.7.

import numpy as np, torch, torch.nn as nn, torch.nn.functional as F
from scipy.fft import idctn
from sklearn.datasets import load_sample_image

torch.manual_seed(0); torch.set_num_threads(1)
rng = np.random.default_rng(0)
sigma, s = 20.0 / 255, 8
img = load_sample_image("china.jpg").mean(axis=2) / 255
test = img[100:356, 200:456]                                         # the crop of ch. 4-7
train = torch.tensor(np.concatenate([img[:, :200], img[:, 456:]], axis=1), dtype=torch.float32)
atoms = idctn(np.eye(s * s).reshape(s * s, s, s), axes=(-2, -1), norm="ortho")[:, None]

class ShallowDCT(nn.Module):
    """Proposition 8.2 with soft thresholds: analysis filters, one threshold per channel, synthesis filters."""
    def __init__(self):
        super().__init__()
        self.wa = nn.Parameter(torch.tensor(atoms, dtype=torch.float32))   # layer 1 starts as the DCT
        self.ws = nn.Parameter(torch.tensor(atoms, dtype=torch.float32))   # layer 2 starts as its inverse
        tau = torch.full((1, s * s, 1, 1), 1.5 * sigma); tau[0, 0] = 0      # the DC channel is kept
        self.tau = nn.Parameter(tau)

    def forward(self, y):
        a = F.conv2d(y, self.wa)
        a = F.relu(a - self.tau) - F.relu(-a - self.tau)                  # soft threshold with ReLUs
        cnt = F.conv_transpose2d(torch.ones_like(a[:, :1]), torch.ones(1, 1, s, s))
        return F.conv_transpose2d(a, self.ws) / cnt                       # uniform aggregation

net = ShallowDCT()
y_test = torch.tensor(test + sigma * rng.standard_normal(test.shape), dtype=torch.float32)[None, None]
psnr = lambda a: 10 * np.log10(1 / np.mean((a - test) ** 2))
with torch.no_grad():
    before = psnr(net(y_test)[0, 0].numpy())
opt = torch.optim.Adam(net.parameters(), lr=1e-3)
H, W = train.shape
for step in range(300):                                               # same budget as ch. 7.7
    i = torch.randint(0, H - 32, (16,)); j = torch.randint(0, W - 32, (16,))
    x = torch.stack([train[a:a + 32, b:b + 32] for a, b in zip(i, j)])[:, None]
    loss = ((net(x + sigma * torch.randn_like(x)) - x) ** 2).mean()
    opt.zero_grad(); loss.backward(); opt.step()
with torch.no_grad():
    after = psnr(net(y_test)[0, 0].numpy())
print(f"fixed DCT, soft 1.5 sigma: {before:.2f} dB | learned two-layer network: {after:.2f} dB"
      f" | {sum(p.numel() for p in net.parameters())} parameters")

On the same noisy crop, the fixed soft-threshold DCT scores 25.96 dB and the trained two-layer network 27.07 dB, with 8256 parameters. That is 0.3 dB above hard-threshold DCT (26.78 dB) and level with the four-layer CNN of ch. 7.7 (27.00 dB). At this small budget, a shallow network shaped by the human design matches a generic deeper one. Depth pays off with more data and training, as in DnCNN.

NL-means does not fit this frame. Its weights depend on the content of the image, through patch distances. Such content-dependent averaging is an attention mechanism; non-local and attention layers bring it back into networks (ch. 9.4).

The next code checks Proposition 8.2 numerically, and the two identities of Proposition 8.4 below.

import numpy as np, torch, torch.nn as nn, torch.nn.functional as F
from scipy.fft import dctn, idctn
from sklearn.datasets import load_sample_image

torch.manual_seed(0); rng = np.random.default_rng(0)
u = load_sample_image("china.jpg").mean(axis=2)[100:228, 200:328]      # 128 x 128 crop
sigma, s = 20.0, 8
f = u + sigma * rng.standard_normal(u.shape)

# 1) Algorithm 4.6 written with NumPy
P = np.lib.stride_tricks.sliding_window_view(f, (s, s))
C = dctn(P, axes=(-2, -1), norm="ortho"); C[np.abs(C) < 3 * sigma] = 0
Q = idctn(C, axes=(-2, -1), norm="ortho")
out, cnt = np.zeros_like(f), np.zeros_like(f)
for i in range(s):
    for j in range(s):
        out[i:i + Q.shape[0], j:j + Q.shape[1]] += Q[:, :, i, j]
        cnt[i:i + Q.shape[0], j:j + Q.shape[1]] += 1
numpy_result = out / cnt

# 2) The same algorithm as a two-layer CNN (Proposition 8.2)
atoms = idctn(np.eye(s * s).reshape(s * s, s, s), axes=(-2, -1), norm="ortho")  # 64 DCT atoms
w = torch.tensor(atoms[:, None], dtype=torch.float64)                # 64 filters of size 8 x 8
y = torch.tensor(f[None, None])
alpha = F.conv2d(y, w)                                               # layer 1: DCT of every patch
alpha = alpha * (alpha.abs() >= 3 * sigma)                           # hard-threshold non-linearity
z = F.conv_transpose2d(alpha, w)                                     # layer 2: inverse DCT + sum
count = F.conv_transpose2d(torch.ones_like(alpha[:, :1]), torch.ones(1, 1, s, s, dtype=torch.float64))
cnn_result = (z / count)[0, 0].numpy()
print("max |NumPy - CNN| =", np.abs(numpy_result - cnn_result).max())

# 3) A bias-free ReLU CNN is positively homogeneous and equals its Jacobian times y (Prop. 8.4)
net = nn.Sequential(nn.Conv2d(1, 16, 3, padding=1, bias=False), nn.ReLU(),
                    nn.Conv2d(16, 16, 3, padding=1, bias=False), nn.ReLU(),
                    nn.Conv2d(16, 1, 3, padding=1, bias=False)).double()
y = torch.randn(1, 1, 32, 32, dtype=torch.float64)
fy, jvp = torch.autograd.functional.jvp(net, y, y)                # f(y) and J_f(y) y
print("f(3.7 y) - 3.7 f(y):", (net(3.7 * y) - 3.7 * fy).abs().max().item())
print("f(y) - J_f(y) y:    ", (fy - jvp).abs().max().item())

The two implementations agree to 2 \times 10^{-13}, and the identities of Proposition 8.4 hold to machine precision.

8.4 Bias-free networks are adaptive linear filters

Read row i of f(y) = J_f(y) y:

\hat x(i) = \sum_j J_f(y)_{ij}\, y(j) .

The network computes a weighted average of the noisy pixels, with weights chosen by the image itself — exactly the form of NL-means, with learned instead of hand-designed weights. Mohan, Kadkhodaie, Simoncelli and Fernandez-Granda (2020) removed the biases of DnCNN-type networks. They observed that bias-free networks generalise to noise levels far outside their training range, where ordinary networks fail, and used the rows of the Jacobian to visualise the adaptive filters.

The Jacobian also gives SURE (Theorem 2.11) for a network: \operatorname{div} f(y) = \mathrm{tr}\, J_f(y). It is the effective number of degrees of freedom the denoiser keeps, which can be estimated with random probes \mathbb{E}[v^\top J_f v] for v \sim \mathcal{N}(0, I).

8.5 Generalisation: noise level and real noise

By Theorem 7.4, a network trained at one noise level approximates \mathbb{E}[X \mid X + N_\sigma] for that \sigma only. Three remedies exist:

• train blindly over a range of levels, as the blind DnCNN;
• give the level as input, as FFDNet;
• constrain the architecture, as the bias-free networks of Proposition 8.4.

Real camera noise is neither white nor Gaussian (ch. 3.2). Plötz and Roth (2017) built a benchmark of real photographs with nearly noise-free references. On it, several methods that excelled on synthetic Gaussian noise lost their advantage over BM3D. Real-noise datasets (Abdelhamed, Lin and Brown, 2018) and the synthesis of realistic raw noise by inverting the camera pipeline (Brooks et al., 2019) answer this mismatch.

9. Comparison of all denoising methods

Every chapter so far produced a number for the same noisy crop. This chapter turns such numbers into a fair comparison, asks how far any method can go, and draws the lessons of the human–machine competition.

9.1 How to compare denoisers

A sound protocol fixes the test images, the noise realisations and the metrics before tuning anything. Its main rules are the following.

• Tune parameters on images that are not in the test set.
• Use the same noise realisation for all methods, so that differences are paired.
• Report per-image results, not only the mean, and test several noise levels.
• Look at the method noise \tilde u - D(\tilde u) and at details: PSNR does not see textures that are smoothed away or invented.

Pairing removes the variability between images, which is much larger than the differences between methods. The code compares three methods on 18 crops of the two scikit-learn images at \sigma = 25.

import numpy as np
from scipy.fft import dctn, idctn
from scipy.ndimage import uniform_filter, gaussian_filter
from sklearn.datasets import load_sample_image

def dct_denoise(f, sigma, s=8):                      # Algorithm 4.6
    C = dctn(np.lib.stride_tricks.sliding_window_view(f, (s, s)), axes=(-2, -1), norm="ortho")
    C[np.abs(C) < 3 * sigma] = 0
    Q = idctn(C, axes=(-2, -1), norm="ortho")
    out, cnt = np.zeros_like(f), np.zeros_like(f)
    for i in range(s):
        for j in range(s):
            out[i:i + Q.shape[0], j:j + Q.shape[1]] += Q[:, :, i, j]
            cnt[i:i + Q.shape[0], j:j + Q.shape[1]] += 1
    return out / cnt

def nl_means(f, sigma, hp=2, r=10, h=0.7):           # Definition 5.2
    g = np.pad(f, r, mode="reflect"); num = np.zeros_like(f); den = np.zeros_like(f)
    for dx in range(-r, r + 1):
        for dy in range(-r, r + 1):
            sh = g[r + dx: r + dx + f.shape[0], r + dy: r + dy + f.shape[1]]
            d2 = uniform_filter((f - sh) ** 2, size=2 * hp + 1, mode="reflect")
            w = np.exp(-np.maximum(d2 - 2 * sigma**2, 0) / (h * sigma) ** 2)
            num += w * sh; den += w
    return num / den

rng = np.random.default_rng(1)
images = [load_sample_image(n).mean(axis=2) for n in ("china.jpg", "flower.jpg")]
crops = [img[i:i + 128, j:j + 128] for img in images for i in (0, 150, 290) for j in (0, 250, 500)]
sigma, rows = 25.0, []
for u in crops:
    f = u + sigma * rng.standard_normal(u.shape)
    p = lambda a: 10 * np.log10(255**2 / np.mean((a - u) ** 2))
    rows.append([p(gaussian_filter(f, 0.8)), p(dct_denoise(f, sigma)), p(nl_means(f, sigma))])
R = np.array(rows)
print("mean PSNR  Gaussian / DCT / NL-means:", np.round(R.mean(0), 2))
diff = R[:, 2] - R[:, 1]                             # paired differences NL-means - DCT
se = diff.std(ddof=1) / np.sqrt(len(diff))
print(f"NL-means - DCT: mean {diff.mean():+.2f} dB, standard error {se:.2f}, t = {diff.mean() / se:.1f},"
      f" NL-means better on {(diff > 0).sum()} of {len(diff)} crops")

The mean PSNR is 26.27 dB for a fixed Gaussian filter, 31.43 dB for the sliding DCT and 31.86 dB for NL-means. NL-means wins on 15 of the 18 crops, by 0.43 dB on average with a standard error of 0.12 dB (t = 3.5). The 18 crops come from two images only, so they are not independent draws: the conclusion holds for these images, not for natural images in general.

Standard test sets are Set12, Kodak24 and BSD68, a subset of the Berkeley segmentation dataset (Martin, Fowlkes, Tal and Malik, 2001). Real-noise benchmarks are discussed in ch. 8.5.

9.2 Inherent bounds

How much better can any denoiser get? Restrict the denoiser of pixel i to a window W of the noisy image around i.

Two studies made such bounds quantitative. Chatterjee and Milanfar (2010) derived lower bounds on the error of patch-based denoisers and concluded that, for many images, some room for improvement remained. Levin and Nadler (2011) estimated the MMSE of denoising from a small window with a very large collection of natural patches. For small windows, the best methods of the time were already close to the bound: further progress requires larger supports — which deep networks provide through their receptive field.

9.3 The landscape

The classical methods carry their hypothesis openly and need no training data; their parameters are few and readable. The learned ones absorb the image statistics of their training set and win when test images resemble it, and lose that edge when the noise or the content changes.

9.4 What each side teaches the other

• Networks rediscovered the classical ingredients. A convolutional layer is a bank of patch transforms (Proposition 8.2), its non-linearity a shrinkage, the sum of overlapping outputs an aggregation, the downsampling of FFDNet or of encoder–decoder networks such as U-Net (Ronneberger, Fischer and Brox, 2015) a multiscale scheme.
• Classical methods explain what networks compute. A bias-free denoiser is an adaptive weighted average (Proposition 8.4), the object NL-means made explicit.
• Hybrids combine both. Non-local operations inside networks bring self-similarity back (Lefkimmiatis, 2017). Learned priors replace the hand-made Gaussian models of ch. 6.
• Denoisers became modules. A good denoiser serves as the prior of other inverse problems and as the engine of generative models (ch. 10).

10. Research directions and open problems

This chapter is for self-study. It shows where the tools of the course lead in current research: training without clean images, denoisers used as priors and as generative models, video, real noise, and the questions that remain open.

10.1 Learning to denoise without clean images

Clean references are often unavailable: in microscopy, astronomy or medical imaging, every image is noisy. Two observations make supervised training possible anyway.

Case 1 is Noise2Noise (Lehtinen et al., 2018): train on pairs of independent noisy shots of the same scene. Case 2 underlies Noise2Void (Krull, Buchholz and Jug, 2019) and Noise2Self (Batson and Royer, 2019), which train on single noisy images with masked pixels. It is the identity already met in Theorem 5.5: a noisy value is an unbiased target for a predictor that does not see it. Theorem 2.11 gives a third route — minimising SURE needs no clean image when the noise is Gaussian with known \sigma.

10.2 Denoisers as priors: Tweedie’s formula and diffusion

An MMSE denoiser, and hence a well-trained network (Theorem 7.4), gives the score \nabla \log p_Y of the noisy image distribution. This is the engine of diffusion models, which generate images by denoising pure noise in many small steps (Ho, Jain and Abbeel, 2020; Song et al., 2021). Kadkhodaie and Simoncelli (2021) used a bias-free denoiser in the same way to sample from the prior implicit in it and to solve linear inverse problems.

The code checks Proposition 10.2 on a one-dimensional Gaussian mixture, where both sides can be computed.

import numpy as np

rng = np.random.default_rng(0)
sigma = 0.5
pi, mu, s = np.array([0.3, 0.7]), np.array([-2.0, 1.0]), np.array([0.3, 0.8])  # 1D Gaussian mixture prior

def p_y(y):                     # density of Y = X + N, again a Gaussian mixture
    v = s**2 + sigma**2
    return (pi * np.exp(-(y[:, None] - mu) ** 2 / (2 * v)) / np.sqrt(2 * np.pi * v)).sum(1)

y = np.linspace(-3, 3, 7)
eps = 1e-5
score = (np.log(p_y(y + eps)) - np.log(p_y(y - eps))) / (2 * eps)   # d/dy log p_Y
tweedie = y + sigma**2 * score                                     # Proposition 10.2

# Monte Carlo posterior mean: sample (X, Y), average X over a small bin around each y
k = rng.choice(2, size=4_000_000, p=pi)
X = mu[k] + s[k] * rng.standard_normal(k.size)
Y = X + sigma * rng.standard_normal(k.size)
mc = np.array([X[np.abs(Y - t) < 0.01].mean() for t in y])
for t, a, b in zip(y, tweedie, mc):
    print(f"y = {t:+.1f}: Tweedie {a:+.3f}   Monte Carlo {b:+.3f}")

The formula and the Monte Carlo posterior mean agree to about 0.004 at all seven points.

Denoisers are also used as priors without a probabilistic model. Plug-and-play methods split an inverse problem — deblurring, super-resolution, tomography — and replace the proximal step of the prior by a denoiser (Venkatakrishnan, Bouman and Wohlberg, 2013). Regularisation by denoising defines an explicit regulariser from a denoiser (Romano, Elad and Milanfar, 2017). The Inverse Problems and Imaging and Diffusion and Multiscale Models in AI courses of semester 2 develop both directions.

10.3 Video and burst denoising

Consecutive video frames show the same scene, so space-time patches have many more copies than spatial ones. VBM4D extends BM3D to space-time volumes (Maggioni, Boracchi, Foi and Egiazarian, 2012). Arias and Morel (2018) extend the empirical Bayesian estimation of ch. 6 to space-time patches. On the learned side, Ehret, Davy, Morel, Facciolo and Arias (2019) fine-tune a pre-trained network on the video to denoise itself, with frame-to-frame training in the spirit of Noise2Noise. FastDVDnet (Tassano, Delon and Veit, 2020) denoises video without explicit motion estimation.

10.4 Real noise and raw images

Most of this course assumes white Gaussian noise of known level. Real images break all three assumptions (ch. 3.2). Current work models the noise of the raw sensor (Foi, Trimeche, Katkovnik and Egiazarian, 2008), synthesises realistic training data by inverting the camera pipeline (Brooks et al., 2019), collects real noisy–clean pairs (Plötz and Roth, 2017; Abdelhamed, Lin and Brown, 2018), and denoises blindly from a single image (Lebrun, Colom and Morel, 2015).

10.5 Open problems

• Understanding learned denoisers. The course description notes that their inner workings are not fully understood. Jacobian analyses of bias-free networks (Proposition 8.4) are one entry point; a theory that predicts what a network learns from a given training set is missing.
• Bounds for large supports. The bounds of ch. 9.2 are tight for small windows. How close current networks are to the MMSE for large receptive fields is unknown.
• Guarantees under distribution shift. Theorem 7.4 holds for the training distribution only. Certifying a denoiser on new noise or new content remains open.
• Hallucination. A denoiser trained on natural images may invent plausible details. In medical and scientific imaging this is a failure, and quantifying it is an active subject.
• Architectures. Transformer-based restoration networks such as SwinIR (Liang et al., 2021) and Restormer (Zamir et al., 2022) bring attention, a learned form of non-local averaging, into the competition.

Exercises

The exercises are sorted by level: ★ applies the course, ★★ combines results or proves a short statement, ★★★ is open-ended. Exercises 5, 9, 11, 14 and 15 are computational.

1. ★ Aliasing by hand. A camera without optical blur samples the stripe pattern \cos(2\pi(0.6\, k + 0.2\, l)) on the integer grid. Which pattern is recorded, and in which direction do its stripes run? Which band-limiting condition of Theorem 1.2 is violated?

   Solution: for integers k, l, \cos(2\pi(0.6k + 0.2l)) = \cos(2\pi(-0.4k + 0.2l)), because 0.6 = 1 - 0.4. The recorded frequency is (-0.4, 0.2) cycle per pixel instead of (0.6, 0.2): a coarser pattern, with stripes in a different direction.

   The frequency 0.6 cycle per pixel, i.e. 1.2\pi radians per pixel, lies outside [-\pi, \pi], so the periodisation of Proposition 1.1 folds it.

2. ★ A 2 × 2 DFT. Compute the DFT of the image u = \begin{pmatrix} 1 & 2 \\ 3 & 4 \end{pmatrix} and check the Parseval identity of Proposition 1.5.

   Solution: for M = N = 2, e^{-2i\pi pk/2} = (-1)^{pk}. Hence \hat u(0,0) = 10, \hat u(1,0) = (1 + 2) - (3 + 4) = -4, \hat u(0,1) = (1 + 3) - (2 + 4) = -2 and \hat u(1,1) = 1 - 2 - 3 + 4 = 0.

   Parseval: \sum |u|^2 = 30 and \frac14(100 + 16 + 4 + 0) = 30.

3. ★ The energy of a noise patch. For \sigma = 10 and 7 \times 7 patches, give the mean and the standard deviation of \|N\|^2. Estimate the probability that \|N\|^2 exceeds 1.5 times its mean.

   Solution: d = 49, so the mean is 49 \times 100 = 4900 and the standard deviation is \sigma^2\sqrt{2d} = 100\sqrt{98} \approx 990 (Proposition 2.5).

   The threshold 7350 is 2450/990 \approx 2.47 standard deviations above the mean: the normal approximation gives 0.007. The exact \chi^2_{49} tail at 73.5 is 0.013 — the \chi^2 law is skewed to the right.

4. ★ Box filters and PSNR. An image is corrupted by noise with \sigma = 30. Compute the PSNR of the noisy image, then of a 3 \times 3 and a 5 \times 5 box filter on a flat region. What happens next to a vertical step of height 90?

   Solution: noisy: 10\log_{10}(65025/900) = 18.59 dB. A k \times k box divides the variance by k^2: 3 \times 3 gives MSE 100, i.e. 28.13 dB; 5 \times 5 gives MSE 36, i.e. 32.57 dB.

   Next to the step, a 3 \times 3 box mixes one column of the other side: the bias is 90/3 = 30, a squared bias of 900 — the error of the noisy image itself. The 5 \times 5 box spreads such errors over four columns (Example 0.4).

5. ★ Anscombe at low counts (Python). For Poisson counts with \lambda \in \{0.5, 1, 2, 5, 20\}, compute \mathbb{E}[2\sqrt{k + 3/8}] with scipy.stats.poisson.pmf, then apply the algebraic inverse a \mapsto (a/2)^2 - 3/8 to this mean. Compare with \lambda and conclude.

   Hints and expected result: the means are 1.742, 2.187, 2.928, 4.527 and 8.972. The algebraic inverse returns 0.383, 0.821, 1.769, 4.749 and 19.75 instead of 0.5, 1, 2, 5 and 20.

   Even a perfect denoiser in the stabilised domain, which returns the mean, would bias the result downwards by about 20 % at low counts. This is why an exact unbiased inverse is used (Mäkitalo and Foi, 2011).

6. ★★ Gaussian conditioning in two dimensions. Let X \sim \mathcal{N}(0, C) with C = \begin{pmatrix} 4 & 2 \\ 2 & 2 \end{pmatrix}, and Y = X + N with \sigma^2 = 1. Compute the gain K, the posterior covariance and the MMSE. Then prove that C - KC = \sigma^2 K in general.

   Solution: C + I = \begin{pmatrix} 5 & 2 \\ 2 & 3 \end{pmatrix} has determinant 11, so

   K = C(C + I)^{-1} = \frac{1}{11}\begin{pmatrix} 8 & 2 \\ 2 & 6 \end{pmatrix}.

   A direct product gives the posterior covariance

   C - KC = \frac{1}{11}\begin{pmatrix} 8 & 2 \\ 2 & 6 \end{pmatrix},

   so the MMSE is 14/11 \approx 1.27, against 2 for Y itself. In general, since C and (C + \sigma^2 I)^{-1} commute,

   C - KC = C(C + \sigma^2 I)^{-1}\big((C + \sigma^2 I) - C\big) = \sigma^2 K .

7. ★★ Two oracles. For x = (10, 3, 1, 0.5) and \sigma = 2, compute the keep-or-kill risk R_P and the attenuation risk R_A of Proposition 4.2. Then show that the lower bound R_A = R_P/2 is attained for a single coefficient with x_k^2 = \sigma^2.

   Solution: R_P = 4 + 4 + 1 + 0.25 = 9.25, and

   R_A = \frac{400}{104} + \frac{36}{13} + \frac{4}{5} + \frac{1}{4.25} = 3.85 + 2.77 + 0.80 + 0.24 = 7.65,

   a ratio of 0.83.

   For one coefficient with x_k^2 = \sigma^2, R_P = \sigma^2 and R_A = \sigma^4/(2\sigma^2) = \sigma^2/2: the lower bound of Proposition 4.2 is attained.

8. ★★ How often does NL-means reject a true copy? Take 7 \times 7 patches, \sigma = 30 and h = 0.7\sigma. For two noisy versions of the same patch, give the law of \delta^2 and the probability that their weight falls below e^{-1}.

   Solution: the difference of the two noisy patches is \mathcal{N}(0, 2\sigma^2 I_{49}), so \delta^2 = \frac{2\sigma^2}{49}\chi^2_{49}, with mean 2\sigma^2 = 1800 and standard deviation 2\sqrt2 \times 900/7 \approx 364 (Proposition 5.3).

   The weight is below e^{-1} when \delta^2 - 2\sigma^2 > h^2 = 441, i.e. \chi^2_{49} > 2241 \times 49/1800 = 61.0. The exact tail is 0.117, the normal approximation 0.113: about one true copy in nine is down-weighted.

9. ★★ SURE inside each patch (Python). In Algorithm 4.6, keep the DC coefficient, replace the hard threshold of the 63 others by a soft threshold, and choose \tau separately in each patch by minimising SURE (Proposition 4.5) over \tau \in \sigma \cdot \{0.5, 0.75, \dots, 3.5\}. Test on the flower.jpg crop of ch. 4 at \sigma = 20 against the hard threshold 3\sigma and a fixed soft threshold 1.5\sigma. Why does SURE not apply to hard thresholding?

   Hints and expected result: about 35.5 dB for the hard threshold, 34.6 dB for the per-patch SURE soft threshold (median \tau = 2\sigma), 33.8 dB for the fixed soft threshold.

   SURE improves soft thresholding, but soft thresholding shrinks the large coefficients and blurs, while aggregation already reduces the variance of hard thresholding. Hard thresholding is discontinuous, so Theorem 2.11 does not apply: its divergence has Dirac masses at \pm\tau.

10. ★★ The empirical Wiener filter of NL-Bayes. (a) Compute the clipped factor of Theorem 6.2 for \mu_i \le \sigma^2, for \mu_i = 1.2\sigma^2 and for \mu_i = 2\sigma^2. (b) Why is \tilde C singular when the group has K \le d patches? (c) Why does the second step of Algorithm 6.4 use the basic estimate for the covariance but the noisy patches as input?

    Solution: (a) The factor (1 - \sigma^2/\mu_i)_+ equals 0, then 1 - 1/1.2 = 1/6, then 1/2.

    2. The centred patches \tilde P_k - \tilde m sum to zero, so they span at most K - 1 dimensions, and \tilde C has rank at most K - 1 < d.

    3. The covariance of the basic estimate is no longer inflated by noise eigenvalues (Example 6.3), but the basic estimate is biased (it is already smoothed). The Wiener formula needs an unbiased observation X + N with known noise, which only the noisy patches provide.

11. ★★ Backpropagation through a convolution (Python). For the cross-correlation without padding z(p) = \sum_q w(q)\, a(p + q), written z = w \star a, and a loss with gradient g = \partial\ell/\partial z, prove that \partial \ell / \partial w(q) = \sum_p g(p)\, a(p + q), i.e. the cross-correlation of a with g. Check it with PyTorch.

    Solution: z(p) = \sum_q w(q)\, a(p + q), so \partial z(p)/\partial w(q) = a(p + q) and the chain rule gives the formula, as in Proposition 7.5.

    import torch, torch.nn.functional as F
    
    torch.manual_seed(0)
    a = torch.randn(1, 1, 6, 6, dtype=torch.float64)                    # input channel
    w = torch.randn(1, 1, 3, 3, dtype=torch.float64, requires_grad=True)
    z = F.conv2d(a, w)                       # 4 x 4 output: z(p) = sum_q w(q) a(p + q)
    g = torch.randn_like(z)                  # upstream gradient dl/dz
    (z * g).sum().backward()                 # loss l = <g, z>, so dl/dz = g
    manual = F.conv2d(a, g)                  # dl/dw(q) = sum_p g(p) a(p + q)
    print("autograd = cross-correlation of a with g:", torch.allclose(w.grad, manual))

    The script prints True.

12. ★★ Receptive field and parameter count. A DnCNN-type network has 17 layers of 3 \times 3 convolutions with biases: 1 \to 64, then fifteen 64 \to 64, then 64 \to 1. Count its parameters, ignoring batch normalisation, and give its receptive field. How many layers would a 60 \times 60 receptive field need? How does the answer change if the image is first rearranged into four half-resolution sub-images, as FFDNet does (ch. 8.2)?

    Solution: 64 \times 10 = 640, plus 15 \times 64 \times 577 = 553\,920, plus 577: in total 555\,137 parameters. The receptive field is 35 \times 35 (Proposition 7.9).

    A 61 \times 61 field needs L = 30 layers. After a 2 \times downsampling into four sub-images, each layer covers twice as many original pixels: about 15 layers suffice, plus the rearrangement.

13. ★★ Homogeneous classical denoisers. (a) Prove \eta_S(t) = \rho(t - \tau) - \rho(-t - \tau). (b) Show that the sliding DCT denoiser with a fixed \tau is not positively homogeneous. (c) Show that it becomes homogeneous of degree 1 if \tau = 3\hat\sigma(\tilde u) with the MAD estimator of Proposition 3.7, and relate this to Proposition 8.4.

    Solution: (a) Check the three cases t > \tau, |t| \le \tau and t < -\tau.

    2. Take \alpha > 0 so small that every DCT coefficient of every patch of \alpha\tilde u is below \tau. Then D(\alpha \tilde u) = 0, while \alpha D(\tilde u) \neq 0 as soon as D(\tilde u) \neq 0.

    3. The MAD estimate scales with the image: multiplying \tilde u by \alpha > 0 multiplies every high-pass coefficient, hence their median and \hat\sigma, by \alpha.

    Then every coefficient and the threshold scale together, so the kept set is unchanged and the output scales by \alpha. Like a bias-free network, this denoiser has no fixed intensity scale and adapts to the noise level.

14. ★★★ Noise-level mismatch (Python). Train the network of ch. 7.7 at \sigma = 20, once with biases and once with bias=False in every convolution. Test both on the crop of ch. 4 at \sigma = 5, 10, 20, 40, with the same noise realisations, and compare with the PSNR of the noisy input. Explain the results with Theorem 7.4 and Proposition 8.4.

    Hints and expected result: the noisy inputs score 34.2, 28.1, 22.1 and 16.1 dB. Call torch.manual_seed(0) before each training and draw each test noise with np.random.default_rng(0). With 300 steps as in ch. 7.7, the network with biases gives about 29.8, 29.3, 27.0 and 19.9 dB. At \sigma = 5 it is worse than its input: it removes noise of the training level from an image that has much less.

    The bias-free network gives about 30.8, 29.8, 26.8 and 20.9 dB. Over other seeds the values move by up to 0.5 dB at \sigma = 5, but the pattern stays. The bias-free network loses at most 0.25 dB at the training level and gains about 1 dB at \sigma = 40, but does not fix the failure at low noise.

    Proposition 8.4 only makes the network consistent under contrast scaling: denoising at \sigma = 5 is like denoising at \sigma = 20 an image of four times the contrast, outside the training distribution. Retrain both networks with \sigma drawn uniformly in [0, 25] for each patch, as Mohan et al. (2020) train over a range starting at zero, and compare again.

    Theorem 7.4 explains the failure: the trained function approximates \mathbb{E}[X \mid X + N_{20}], not the denoiser for other levels.

15. ★★★ Choosing a parameter without clean images (Python). Take two independent noisy shots y_1, y_2 of the flower.jpg crop at \sigma = 20. For a Gaussian filter of width h \in \{0.5, 1, \dots, 3\}, compare the Noise2Noise loss \frac{1}{|\Omega|}\|G_h * y_1 - y_2\|^2 with \mathrm{MSE}(G_h * y_1) + \sigma^2. Which h does each criterion select? Then design a blind-spot version of NL-means and explain which pixels must be excluded from the patch distances.

    Hints and expected result: the two columns differ by less than 5 at every h and are both minimal at h = 1.5 (Noise2Noise loss 440.8, true MSE plus \sigma^2 437.4), as Proposition 10.1 predicts.

    For NL-means, the value at i must not depend on \tilde u(i). Exclude j = i from the average, and remove the pixel i from both patches in every distance \delta^2(i, j).

    The pixel i appears in \tilde P_i at the centre, and in \tilde P_j whenever j is within the patch radius of i.

Projects and assessment

The official assessment (master-mva.com) has three parts: simple mathematical exercises and an experimental report on online image-processing experiments, both delivered each week and requiring no programming, plus Colab notebook exercises that are not evaluated; an intermediate written exam; and a final oral examination, written if the number of students is too large. The weights are not displayed. The official assessment contains no project. The subjects below are proposals for practice, for a research-oriented extension of the course, or for preparing the oral examination; the exam proposal mirrors the intermediate written exam.

Solo projects

1. [Solo · ≈ 20 h] Paper study: NL-Bayes. Reproduce the two-step algorithm of Lebrun, Buades and Morel (2013, SIAM Journal on Imaging Sciences, and its IPOL implementation) on grey images.
   • Implement both steps, the flat-area test and the aggregation; vary the patch size and K.
   • Data: five public test images at \sigma = 10, 20, 40, with fixed noise realisations.
   • Reference: the C++ source code published with the IPOL article, run locally on the same noisy images.
   • Success criterion: within 0.3 dB of the reference implementation, and a measured gain of the second step.
   • Deliverables: code, a 4–6 page report, a 10-minute talk. Chapters 2, 5 and 6.
2. [Solo · ≈ 20 h] Multiscale DCT and low-frequency noise. Measure the residual noise of Algorithm 4.6 and how Algorithm 4.10 removes it.
   • Estimate the power spectrum of the error \hat u - u on flat images, then with two and three scales and several recomposition fractions r.
   • Compare with the analysis of Pierazzo, Morel and Facciolo (2017).
   • Success criterion: a figure showing the low-frequency excess predicted by Proposition 4.8 and its removal, with PSNR gains on smooth images.
   • Deliverables: code, a 4–6 page report, a 10-minute talk. Chapters 1, 2 and 4.
3. [Solo · ≈ 20 h] Estimating a noise curve. Recover the parameters (a, b) of the Poisson–Gaussian model of Definition 3.5 from a single image.
   • Generate raw-like images with known (a, b) from the scikit-learn sample images; implement the MAD estimator per intensity bin and a low-variance-block DCT estimator.
   • Success criterion: relative error below 10 % on a and b for at least three settings.
   • Deliverables: code, a 4–6 page report, a 10-minute talk. Chapters 2 and 3.

Group projects

1. [Group of 2–3 · ≈ 30 h each] Paper study: bias-free denoisers. Reproduce the main experiment of Mohan, Kadkhodaie, Simoncelli and Fernandez-Granda (2020): networks trained on a limited range of noise levels, tested far outside it.
   • Data: a public set of natural images such as the Berkeley segmentation dataset; a free GPU tier suffices for DnCNN-size networks.
   • Split: one member builds the data pipeline and trains the networks with biases; one trains the bias-free versions; one evaluates with paired statistics and computes Jacobian rows (Proposition 8.4).
   • Protocol: train both networks on a limited range such as \sigma \in [0, 10] and test from \sigma = 5 to 100 on held-out images.
   • Success criterion: a figure of PSNR against test \sigma in which the two networks agree within 0.3 dB inside the training range and the bias-free one leads by several dB far outside it, the qualitative finding of the paper.
   • Deliverables: code, a 6–8 page report, a 20-minute talk. Chapters 7, 8 and 9.
2. [Group of 2–3 · ≈ 30 h each] The human–machine benchmark. Compare sliding DCT, multiscale DCT, NL-means, NL-Bayes, BM3D (a public implementation) and a DnCNN-type network trained by the group.
   • Protocol of ch. 9.1: Kodak24 or Set12, \sigma = 15, 25, 50, paired differences with standard errors, SSIM, method noise, run time.
   • Two stress tests: a noise level not seen in training, and correlated noise obtained by blurring white noise.
   • Split: classical methods; network training; evaluation and statistics.
   • Success criterion: a complete table with standard errors and an explanation of each ranking change by a principle of the course.
   • Deliverables: code, a 6–8 page report, a 20-minute talk. Chapters 4–9.
3. [Group of 2–3 · ≈ 30 h each] Denoising without clean images. Compare Noise2Noise and a blind-spot method (Proposition 10.1) with supervised training and NL-Bayes on synthetic fluorescence-like images with Poisson–Gaussian noise.
   • Use the generalised Anscombe transform of ch. 3.3 for the classical baseline.
   • Split: data simulation and variance stabilisation; self-supervised training; evaluation.
   • Success criterion: Noise2Noise within 0.5 dB of supervised training at equal architecture, and the larger gap of the blind-spot method measured and explained (it never sees the centre pixel).
   • Deliverables: code, a 6–8 page report, a 20-minute talk. Chapters 3, 6, 7 and 10.

Written exam (proposal)

Two hours, one handwritten A4 sheet allowed, no calculator needed; four exercises.

1. Question (4 points). State the Shannon–Whittaker theorem and explain aliasing with a numerical example. Why are image patches transformed with the DCT rather than the DFT? Expected: Theorem 1.2 with the band-limiting hypothesis; the periodisation of Proposition 1.1 and a folded frequency such as 0.7 \to 0.3 cycle per pixel; the border jump of the periodic DFT against the symmetric extension of the DCT (Proposition 1.7) and the resulting energy compaction.

2. Question (6 points). A patch model has principal variances \lambda = (9, 4, 1) and the noise has \sigma^2 = 4. Give the MMSE estimator in the principal basis, its error, the error of the noisy observation, and the risks of the two oracles of Proposition 4.2 for a patch whose coefficients are x = (3, 2, 1).

   Expected: shrinkage factors 9/13, 1/2 and 1/5 (Corollary 2.8); MMSE 36/13 + 2 + 0.8 \approx 5.57, against 12 for the observation.

   For x, R_P = 4 + 4 + 1 = 9 and R_A = 36/13 + 2 + 0.8 \approx 5.57, the same value as the MMSE because x_k^2 = \lambda_k.

3. Question (5 points). Prove Proposition 5.3 on the distance between noisy patches. Deduce why NL-means subtracts 2\sigma^2, and how the reliability of the comparison depends on the patch size.

   Expected: the decomposition \Delta + M with M \sim \mathcal{N}(0, 2\sigma^2 I) and the mean and variance computed from the \chi^2 moments. Then the bias 2\sigma^2 per pixel, and the fluctuation 2\sqrt2\sigma^2/\sqrt d that shrinks with d (Proposition 5.6).

4. Question (5 points). Show that sliding-window DCT hard thresholding is a two-layer convolutional network, and that a ReLU network without biases satisfies f(y) = J_f(y)\, y. Use both facts to compare the DCT denoiser, NL-means and DnCNN. Expected: the proof of Proposition 8.2 (filters = atoms, threshold non-linearity, transposed convolution, aggregation); homogeneity and Euler’s identity (Proposition 8.4); fixed filters versus content-adaptive weights versus learned adaptive filters.

Grading rubric (proposal)

References

Official material of the course (master-mva.com, 2026–2027)

• Lecture notes in English with exercises, provided to the students of the course (“Material provided”).
• IPOL — Image Processing On Line, www.ipol.im: the online workshops used in every session.

IPOL articles with online demos used in the notes

• A. Buades, B. Coll, J.-M. Morel, “Non-Local Means Denoising”, Image Processing On Line, 2011.
• G. Yu, G. Sapiro, “DCT Image Denoising: a Simple and Effective Image Denoising Algorithm”, Image Processing On Line, 2011.
• M. Lebrun, “An Analysis and Implementation of the BM3D Image Denoising Method”, Image Processing On Line, 2012.
• M. Lebrun, A. Buades, J.-M. Morel, “Implementation of the ‘Non-Local Bayes’ (NL-Bayes) Image Denoising Algorithm”, Image Processing On Line, 2013.
• M. Colom, A. Buades, “Analysis and Extension of the Ponomarenko et al. Method, Estimating a Noise Curve from a Single Image”, Image Processing On Line, 2013.
• M. Lebrun, M. Colom, J.-M. Morel, “The Noise Clinic: a Blind Image Denoising Algorithm”, Image Processing On Line, 2015.
• N. Pierazzo, J.-M. Morel, G. Facciolo, “Multi-Scale DCT Denoising”, Image Processing On Line, 2017.

Textbooks and surveys

• S. Mallat, A Wavelet Tour of Signal Processing: The Sparse Way, 3rd ed., Academic Press, 2009 — Fourier, sampling, thresholding estimators (ch. 1, 4).
• T. M. Cover, J. A. Thomas, Elements of Information Theory, 2nd ed., Wiley, 2006 — the Gaussian channel (ch. 3).
• I. Goodfellow, Y. Bengio, A. Courville, Deep Learning, MIT Press, 2016 — ch. 7.
• M. Lebrun, M. Colom, A. Buades, J.-M. Morel, “Secrets of image denoising cuisine”, Acta Numerica 21, 2012.
• A. Buades, B. Coll, J.-M. Morel, “A review of image denoising algorithms, with a new one”, Multiscale Modeling and Simulation 4(2), 2005.
• P. Milanfar, “A tour of modern image filtering”, IEEE Signal Processing Magazine 30(1), 2013.

Sampling, noise and statistics (ch. 1–3)

• C. E. Shannon, “A mathematical theory of communication”, Bell System Technical Journal 27, 1948; “Communication in the presence of noise”, Proceedings of the IRE 37(1), 1949.
• N. Wiener, Extrapolation, Interpolation, and Smoothing of Stationary Time Series, MIT Press and Wiley, 1949.
• D. J. Field, “Relations between the statistics of natural images and the response properties of cortical cells”, Journal of the Optical Society of America A 4(12), 1987; D. L. Ruderman, W. Bialek, “Statistics of natural images: scaling in the woods”, Physical Review Letters 73, 1994.
• C. M. Stein, “Estimation of the mean of a multivariate normal distribution”, Annals of Statistics 9(6), 1981.
• F. J. Anscombe, “The transformation of Poisson, binomial and negative-binomial data”, Biometrika 35, 1948; M. Mäkitalo, A. Foi, “Optimal inversion of the Anscombe transformation in low-count Poisson image denoising”, IEEE Transactions on Image Processing 20(1), 2011.
• A. Foi, M. Trimeche, V. Katkovnik, K. Egiazarian, “Practical Poissonian-Gaussian noise modeling and fitting for single-image raw-data”, IEEE Transactions on Image Processing 17(10), 2008.
• V. A. Marchenko, L. A. Pastur, “Distribution of eigenvalues for some sets of random matrices”, Mathematics of the USSR-Sbornik 1(4), 1967.

Classical denoising principles (ch. 3–6)

• L. P. Yaroslavsky, Digital Picture Processing: An Introduction, Springer, 1985; C. Tomasi, R. Manduchi, “Bilateral filtering for gray and color images”, ICCV 1998.
• P. Perona, J. Malik, “Scale-space and edge detection using anisotropic diffusion”, IEEE Transactions on Pattern Analysis and Machine Intelligence 12(7), 1990; L. Rudin, S. Osher, E. Fatemi, “Nonlinear total variation based noise removal algorithms”, Physica D 60, 1992.
• D. L. Donoho, I. M. Johnstone, “Ideal spatial adaptation by wavelet shrinkage”, Biometrika 81(3), 1994; “Adapting to unknown smoothness via wavelet shrinkage”, Journal of the American Statistical Association 90, 1995.
• R. R. Coifman, D. L. Donoho, “Translation-invariant de-noising”, in Wavelets and Statistics, Lecture Notes in Statistics 103, Springer, 1995.
• A. Buades, B. Coll, J.-M. Morel, “A non-local algorithm for image denoising”, CVPR 2005.
• K. Dabov, A. Foi, V. Katkovnik, K. Egiazarian, “Image denoising by sparse 3-D transform-domain collaborative filtering”, IEEE Transactions on Image Processing 16(8), 2007.
• M. Lebrun, A. Buades, J.-M. Morel, “A nonlocal Bayesian image denoising algorithm”, SIAM Journal on Imaging Sciences 6(3), 2013.
• D. Zoran, Y. Weiss, “From learning models of natural image patches to whole image restoration”, ICCV 2011; G. Yu, G. Sapiro, S. Mallat, “Solving inverse problems with piecewise linear estimators: from Gaussian mixture models to structured sparsity”, IEEE Transactions on Image Processing 21(5), 2012.
• A. B. Lee, K. S. Pedersen, D. Mumford, “The nonlinear statistics of high-contrast patches in natural images”, International Journal of Computer Vision 54, 2003; G. Carlsson, T. Ishkhanov, V. de Silva, A. Zomorodian, “On the local behavior of spaces of natural images”, International Journal of Computer Vision 76, 2008; M. Zontak, M. Irani, “Internal statistics of a single natural image”, CVPR 2011.

Neural networks and learned denoisers (ch. 7–9)

• G. Cybenko, “Approximation by superpositions of a sigmoidal function”, Mathematics of Control, Signals and Systems 2, 1989; M. Leshno, V. Ya. Lin, A. Pinkus, S. Schocken, “Multilayer feedforward networks with a nonpolynomial activation function can approximate any function”, Neural Networks 6(6), 1993.
• D. E. Rumelhart, G. E. Hinton, R. J. Williams, “Learning representations by back-propagating errors”, Nature 323, 1986; D. P. Kingma, J. Ba, “Adam: a method for stochastic optimization”, ICLR 2015.
• K. He, X. Zhang, S. Ren, J. Sun, “Deep residual learning for image recognition”, CVPR 2016; S. Ioffe, C. Szegedy, “Batch normalization”, ICML 2015; O. Ronneberger, P. Fischer, T. Brox, “U-Net: convolutional networks for biomedical image segmentation”, MICCAI 2015.
• H. C. Burger, C. J. Schuler, S. Harmeling, “Image denoising: can plain neural networks compete with BM3D?”, CVPR 2012.
• K. Zhang, W. Zuo, Y. Chen, D. Meng, L. Zhang, “Beyond a Gaussian denoiser: residual learning of deep CNN for image denoising”, IEEE Transactions on Image Processing 26(7), 2017; K. Zhang, W. Zuo, L. Zhang, “FFDNet: toward a fast and flexible solution for CNN-based image denoising”, IEEE Transactions on Image Processing 27(9), 2018.
• U. Schmidt, S. Roth, “Shrinkage fields for effective image restoration”, CVPR 2014; Y. Chen, T. Pock, “Trainable nonlinear reaction diffusion: a flexible framework for fast and effective image restoration”, IEEE Transactions on Pattern Analysis and Machine Intelligence 39(6), 2017.
• S. Mohan, Z. Kadkhodaie, E. P. Simoncelli, C. Fernandez-Granda, “Robust and interpretable blind image denoising via bias-free convolutional neural networks”, ICLR 2020.
• S. Lefkimmiatis, “Non-local color image denoising with convolutional neural networks”, CVPR 2017.
• Z. Wang, A. C. Bovik, H. R. Sheikh, E. P. Simoncelli, “Image quality assessment: from error visibility to structural similarity”, IEEE Transactions on Image Processing 13(4), 2004; D. Martin, C. Fowlkes, D. Tal, J. Malik, “A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics”, ICCV 2001.
• P. Chatterjee, P. Milanfar, “Is denoising dead?”, IEEE Transactions on Image Processing 19(4), 2010; A. Levin, B. Nadler, “Natural image denoising: optimality and inherent bounds”, CVPR 2011.
• T. Plötz, S. Roth, “Benchmarking denoising algorithms with real photographs”, CVPR 2017; A. Abdelhamed, S. Lin, M. S. Brown, “A high-quality denoising dataset for smartphone cameras”, CVPR 2018; T. Brooks, B. Mildenhall, T. Xue, J. Chen, D. Sharlet, J. T. Barron, “Unprocessing images for learned raw denoising”, CVPR 2019.

Research directions (ch. 10)

• J. Lehtinen, J. Munkberg, J. Hasselgren, S. Laine, T. Karras, M. Aittala, T. Aila, “Noise2Noise: learning image restoration without clean data”, ICML 2018; A. Krull, T.-O. Buchholz, F. Jug, “Noise2Void — learning denoising from single noisy images”, CVPR 2019; J. Batson, L. Royer, “Noise2Self: blind denoising by self-supervision”, ICML 2019.
• K. Miyasawa, “An empirical Bayes estimator of the mean of a normal population”, Bulletin of the International Statistical Institute 38, 1961; B. Efron, “Tweedie’s formula and selection bias”, Journal of the American Statistical Association 106, 2011.
• J. Ho, A. Jain, P. Abbeel, “Denoising diffusion probabilistic models”, NeurIPS 2020; Y. Song, J. Sohl-Dickstein, D. P. Kingma, A. Kumar, S. Ermon, B. Poole, “Score-based generative modeling through stochastic differential equations”, ICLR 2021; Z. Kadkhodaie, E. P. Simoncelli, “Stochastic solutions for linear inverse problems using the prior implicit in a denoiser”, NeurIPS 2021.
• S. V. Venkatakrishnan, C. A. Bouman, B. Wohlberg, “Plug-and-play priors for model based reconstruction”, IEEE GlobalSIP 2013; Y. Romano, M. Elad, P. Milanfar, “The little engine that could: regularization by denoising (RED)”, SIAM Journal on Imaging Sciences 10(4), 2017.
• M. Maggioni, G. Boracchi, A. Foi, K. Egiazarian, “Video denoising, deblocking, and enhancement through separable 4-D nonlocal spatiotemporal transforms”, IEEE Transactions on Image Processing 21(9), 2012; P. Arias, J.-M. Morel, “Video denoising via empirical Bayesian estimation of space-time patches”, Journal of Mathematical Imaging and Vision 60(1), 2018.
• T. Ehret, A. Davy, J.-M. Morel, G. Facciolo, P. Arias, “Model-blind video denoising via frame-to-frame training”, CVPR 2019; M. Tassano, J. Delon, T. Veit, “FastDVDnet: towards real-time deep video denoising without flow estimation”, CVPR 2020.
• J. Liang, J. Cao, G. Sun, K. Zhang, L. Van Gool, R. Timofte, “SwinIR: image restoration using Swin Transformer”, ICCV Workshops 2021; S. W. Zamir, A. Arora, S. Khan, M. Hayat, F. S. Khan, M.-H. Yang, “Restormer: efficient transformer for high-resolution image restoration”, CVPR 2022.