"They handed me a blurry, noisy, half-missing photograph and asked for the truth. I gave them the most probable truth, which is the only kind I stock."
A Bayesian Restoration Filter
Chapter Overview
Every photograph you have ever taken was damaged before you saw it. Photons arrived at the sensor randomly, so every pixel carries shot noise. The lens spread each point of light into a small smear, and if your hand moved, the smear grew a tail. The scene's brightness range exceeded what the sensor could record, so the highlights clipped or the shadows drowned. And somewhere between sensor and screen, compression quietly threw information away. Image restoration is the discipline of undoing this damage, and image enhancement is its pragmatic cousin: making the picture look better even when "the original" is not strictly recoverable.
What makes restoration intellectually special is that it is the book's first true inverse problem. Everything in Chapter 2 through Chapter 6 ran forward: take an image, apply an operation, get a result. Restoration runs the movie backward. We observe a degraded image $g$, we model how it was produced from an ideal image $f$ (typically blur plus noise, $g = h \ast f + n$), and we try to invert the model. The catch, and the recurring theme of this chapter, is that the inversion is ill-posed: many different clean images explain the same dirty one, and naive inversion amplifies noise catastrophically. Every method here resolves the ambiguity the same way, by adding a prior: an assumption about what natural images look like. Smoothness, self-similarity, sparse gradients, each classical algorithm is a hand-written belief about images. The three-word recipe in the callout below, model, prior, invert, is the spine that every section reuses; keep it in mind as the chapter changes nouns.
The journey runs from damage to remedy. Section 7.1 builds the vocabulary of degradation itself: where noise comes from physically, why its statistics are Poisson at heart, and how realistic degradation pipelines are assembled. Section 7.2 attacks noise, climbing from the Gaussian blur of Chapter 3 to the patch-based elegance of non-local means and BM3D. Section 7.3 takes on blur with frequency-domain deconvolution, where the Fourier machinery of Chapter 4 earns its keep: the Wiener filter and Richardson-Lucy iteration, including the algorithm that rescued the Hubble Space Telescope. Section 7.4 fills holes: inpainting by smoothness, by fast marching, and by copying texture. Section 7.5 asks how much resolution can honestly be recovered, from multi-frame fusion to the eve of deep super-resolution. Section 7.6 closes with dynamic range: merging bracketed exposures into high dynamic range radiance and tone mapping the result back onto an eight-bit screen.
This chapter is also a load-bearing wall for the rest of the book. The noise models of Section 7.1 return as the noise schedules of diffusion models in Chapter 33; denoising itself becomes the training objective of denoising autoencoders in Chapter 31 and, astonishingly, the generative engine of diffusion. Inpainting reappears as generative editing in Chapter 35, and super-resolution returns both at the edge in Chapter 28 and inside text-to-image upscalers in Chapter 34. When deep models eventually beat every method in this chapter, they did it by learning the priors we are about to write by hand. Knowing the hand-written versions is what lets you understand, debug, and distrust the learned ones.
Restoration is inference, not filtering: you model how the damage happened, you state what clean images look like, and the algorithm follows from those two commitments. Every method in this chapter, from the Wiener filter to exemplar inpainting, is the same Bayesian sentence with different nouns: find the image that best explains the observation while remaining plausible. Deep restoration networks did not change the sentence; they replaced the hand-written plausibility term with a learned one.
Every section below is the same three-step recipe wearing different clothes. Model the damage as a forward process (the degradation $g = h \ast f + n$ and its relatives). Name the prior: state what a clean image looks like (smooth, self-similar, sparse, well-exposed). Invert, almost always by iterating the forward model rather than dividing it out. If you can recite "model, prior, invert" you can reconstruct the shape of every algorithm in the chapter, and explain why a Wiener filter, an inpainter, and a diffusion restorer are cousins rather than strangers.
Learning Objectives
- Write down the degradation model $g = h \ast f + n$, identify the physical source of each term, and synthesize realistic noise (Gaussian, Poisson, salt-and-pepper, speckle) and full degradation pipelines in code.
- Denoise images with Gaussian, median, bilateral, non-local means, and BM3D filtering, and explain the prior each method encodes and where it fails.
- Deblur images by inverse filtering, the Wiener filter, and Richardson-Lucy iteration, and explain why regularization is unavoidable.
- Fill missing regions with PDE-based, fast-marching, and exemplar-based inpainting, and choose the right method from hole geometry.
- Reconstruct higher-resolution images from multiple aliased frames and explain why single-image interpolation cannot add information.
- Merge bracketed exposures into HDR radiance maps, tone map them for display, and decide when exposure fusion should replace the radiometric pipeline entirely.
Prerequisites
This chapter leans on most of Part I. From Chapter 1: Digital Image Fundamentals you need the sensor pipeline (where photons become numbers) and the PSNR and SSIM quality metrics, which we use to score every restoration. Chapter 3: Spatial Filtering & Convolution supplies convolution itself, plus the Gaussian, median, and bilateral filters that Section 7.2 builds upon. Chapter 4: The Frequency Domain & Multi-Scale Analysis is essential for Section 7.3 (deconvolution lives in the Fourier domain) and useful for Sections 7.5 and 7.6 (aliasing and pyramid blending). Chapter 5: Geometric Transformations & Image Warping provides interpolation and registration, which multi-frame super-resolution and HDR alignment both require. Chapter 2: Point Operations, Histograms & Thresholding helps with tone mapping, and Chapter 6: Morphology, Binary Images & Shape with preparing inpainting masks. If your environment needs setting up, start at Chapter 0.
Chapter Roadmap
- 7.1 Noise Models & Degradation Pipelines Where damage comes from: shot noise, read noise, blur, and compression, and how to model, synthesize, and measure each one before attempting any cure.
- 7.2 Classical Denoising: From Gaussian to Non-Local Means Denoising as the art of choosing what to average: local filters, edge-aware filters, patch self-similarity, and BM3D, the decade-long classical champion.
- 7.3 Deblurring & Deconvolution: Wiener & Richardson-Lucy Undoing convolution in the Fourier domain: why naive inversion explodes, how Wiener regularizes it, and the iterative method that fixed Hubble's eyesight.
- 7.4 Inpainting: Filling the Holes Reconstructing missing pixels by continuing structure, marching from the boundary, or copying texture, and knowing which trick the hole demands.
- 7.5 Classical Super-Resolution Why bicubic upscaling is not super-resolution, how sub-pixel shifts across frames buy genuine detail, and the back-projection loop that fuses them.
- 7.6 HDR Imaging & Tone Mapping Capturing more light range than the sensor allows and squeezing it onto a display: radiance recovery, tone-mapping operators, and exposure fusion.
The most famous deblurring job in history was performed on a telescope, not by one. When the Hubble Space Telescope launched in 1990 with a mirror ground 2.2 micrometers too flat, astronomers spent three years deconvolving its images with the Richardson-Lucy algorithm (Section 7.3) while waiting for the repair mission. The algorithm dates from 1972, the telescope from 1990, and the combination produced publishable science from a broken billion-dollar instrument. Software patience beat hardware perfection.
Hands-On Lab: Build a Restoration Pipeline and Scorecard
Objective
Build a single script that takes one clean image, runs it forward through a realistic degradation (blur plus Poisson and Gaussian noise), then runs it backward through denoising and deconvolution, and finally scores every stage with PSNR and SSIM. The artifact you keep is a labeled comparison panel and a printed scorecard that proves, in numbers, which step of the model, prior, invert recipe recovered the most quality.
What You'll Practice
- Synthesizing a realistic degradation pipeline, $g = h \ast f + n$, with a known blur kernel and mixed Poisson and Gaussian noise (Section 7.1).
- Denoising with non-local means and reasoning about the self-similarity prior it encodes (Section 7.2).
- Deblurring with Richardson-Lucy iteration and watching why the iteration count is a regularizer (Section 7.3).
- Measuring restoration quality with PSNR and SSIM rather than trusting your eyes (Chapter 1).
- Comparing your from-scratch pipeline against a one-call library equivalent, the book's Right Tool pattern.
Setup
Two libraries cover every step. Install with:
pip install scikit-image numpyThe lab uses the built-in skimage.data.camera() image, so no download is needed; swap in your own grayscale photo if you prefer.
Steps
Step 1: Load a clean reference and a known blur kernel
Restoration is an inverse problem, so you need ground truth to score against. Load a clean image as float in the range zero to one, and build a small Gaussian point-spread function that you will both apply and later invert.
import numpy as np
from skimage import data, img_as_float
from scipy.signal import convolve2d
rng = np.random.default_rng(0)
clean = img_as_float(data.camera()) # 512x512 float image in [0, 1]
def gaussian_psf(size=9, sigma=2.0):
ax = np.arange(size) - size // 2
xx, yy = np.meshgrid(ax, ax)
# TODO: build a 2D Gaussian from xx and yy, then normalize it to sum to 1
# so the blur preserves overall brightness. Hint: np.exp(-(xx**2+yy**2)/(2*sigma**2)).
...
psf = gaussian_psf()
print("psf sums to", psf.sum()) # should print 1.0Hint
k = np.exp(-(xx**2 + yy**2) / (2 * sigma**2)); return k / k.sum(). Normalizing to sum one keeps the blurred image at the same average brightness as the original.
Step 2: Run the image forward through the degradation model
Apply the forward model of Section 7.1: convolve with the PSF, then corrupt with photon (Poisson) noise and a small additive Gaussian read noise. This is the dirty image your pipeline must repair.
def degrade(image, psf, photons=120.0, read_sigma=0.01):
blurred = convolve2d(image, psf, mode="same", boundary="symm")
# TODO: add Poisson shot noise, then additive Gaussian read noise.
# Scale to a photon count first: rng.poisson(blurred * photons) / photons
# gives signal-dependent noise; then add rng.normal(0, read_sigma, ...).
...
return np.clip(noisy, 0, 1)
observed = degrade(clean, psf)Hint
shot = rng.poisson(np.clip(blurred, 0, 1) * photons) / photons turns brightness into a photon count and back, so bright regions get more noise. Then noisy = shot + rng.normal(0, read_sigma, shot.shape).
Step 3: Write the scorecard function
You cannot improve what you cannot measure. Wrap PSNR and SSIM (Chapter 1) into one helper so every later stage reports the same two numbers against the clean reference.
from skimage.metrics import peak_signal_noise_ratio as psnr
from skimage.metrics import structural_similarity as ssim
def score(reference, candidate, label):
# TODO: compute psnr and ssim of candidate against reference (data_range=1.0),
# print them on one line prefixed by label, and return the (psnr, ssim) tuple.
...
score(clean, observed, "degraded ")Hint
p = psnr(reference, candidate, data_range=1.0); s = ssim(reference, candidate, data_range=1.0); print(f"{label}: PSNR {p:5.2f} dB SSIM {s:.3f}"); return p, s.
Step 4: Denoise with non-local means
Apply the self-similarity prior of Section 7.2. Estimate the noise level, then run non-local means, which averages patches that look alike anywhere in the image rather than just nearby pixels.
from skimage.restoration import denoise_nl_means, estimate_sigma
sigma_est = estimate_sigma(observed)
# TODO: call denoise_nl_means on observed with h tied to sigma_est
# (h = 0.8 * sigma_est is a good start), patch_size=5, patch_distance=6,
# fast_mode=True. Store the result in `denoised`, then score it.
...
score(clean, denoised, "denoised ")Hint
denoised = denoise_nl_means(observed, h=0.8 * sigma_est, patch_size=5, patch_distance=6, fast_mode=True). The h parameter sets how aggressively similar patches are averaged.
Step 5: Deblur with Richardson-Lucy
Now invert the blur with the iterative method of Section 7.3. Richardson-Lucy multiplies the estimate toward agreement with the observation; the iteration count is itself a regularizer, since too many iterations re-amplify noise.
from skimage.restoration import richardson_lucy
# TODO: run richardson_lucy on the denoised image using the same psf,
# with num_iter=20. Clip the result to [0, 1] and score it as "restored".
...
score(clean, restored, "restored ")Hint
restored = np.clip(richardson_lucy(denoised, psf, num_iter=20), 0, 1). Try sweeping num_iter over 5, 20, and 100 and watch SSIM rise then fall, the classic semi-convergence of Section 7.3.
Step 6: Assemble a comparison panel and save the artifact
Stack the four stages side by side, clean, degraded, denoised, restored, into one image you can keep and show. This panel plus the printed scorecard is the deliverable.
from skimage.io import imsave
from skimage import img_as_ubyte
divider = np.ones((clean.shape[0], 4)) # thin white separator column
panel = np.hstack([clean, divider, observed, divider, denoised, divider, restored])
# TODO: convert panel to 8-bit and save it as "restoration_panel.png".
# Hint: img_as_ubyte after clipping to [0, 1], then imsave.
...Hint
imsave("restoration_panel.png", img_as_ubyte(np.clip(panel, 0, 1))). Clipping before img_as_ubyte avoids the warning it raises on out-of-range floats.
Expected Output
One PNG file, restoration_panel.png, showing the clean, degraded, denoised, and restored images in a row, and a four-line scorecard on the console. The numbers will look roughly like degraded : PSNR 20.81 dB SSIM 0.402, denoised : PSNR 24.93 dB SSIM 0.642, restored : PSNR 26.10 dB SSIM 0.704 (exact values depend on the noise draw and your library version). The key lesson is the direction of travel: denoising should raise both metrics, and deconvolution should raise them again, confirming that each stage of model, prior, invert adds recoverable quality. If Richardson-Lucy lowers SSIM, you are over-iterating and re-amplifying noise.
Stretch Goals
- Library shortcut: replace your non-local means call with
skimage.restoration.denoise_tv_chambolle(total-variation denoising, a sparse-gradient prior) and compare the scorecard. This mirrors the book's Right Tool principle, swapping one hand-named prior for another in a single line. - Sweep the deconvolution: run Richardson-Lucy for 5, 20, 50, and 100 iterations, plot SSIM against iteration count, and locate the peak that demonstrates the semi-convergence of Section 7.3.
- Reach forward to Chapter 33: note that your Poisson-plus-Gaussian degradation is a single fixed noise level, whereas a diffusion model is trained to denoise across an entire schedule of noise levels. Re-run Step 4 at three different
read_sigmavalues to feel why a learned denoiser needs the whole range.
Complete Solution
import numpy as np
from skimage import data, img_as_float, img_as_ubyte
from skimage.io import imsave
from skimage.restoration import (denoise_nl_means, estimate_sigma,
richardson_lucy)
from skimage.metrics import peak_signal_noise_ratio as psnr
from skimage.metrics import structural_similarity as ssim
from scipy.signal import convolve2d
rng = np.random.default_rng(0)
clean = img_as_float(data.camera())
def gaussian_psf(size=9, sigma=2.0):
ax = np.arange(size) - size // 2
xx, yy = np.meshgrid(ax, ax)
k = np.exp(-(xx**2 + yy**2) / (2 * sigma**2))
return k / k.sum()
def degrade(image, psf, photons=120.0, read_sigma=0.01):
blurred = convolve2d(image, psf, mode="same", boundary="symm")
shot = rng.poisson(np.clip(blurred, 0, 1) * photons) / photons
noisy = shot + rng.normal(0, read_sigma, shot.shape)
return np.clip(noisy, 0, 1)
def score(reference, candidate, label):
p = psnr(reference, candidate, data_range=1.0)
s = ssim(reference, candidate, data_range=1.0)
print(f"{label}: PSNR {p:5.2f} dB SSIM {s:.3f}")
return p, s
psf = gaussian_psf()
observed = degrade(clean, psf)
score(clean, observed, "degraded ")
sigma_est = estimate_sigma(observed)
denoised = denoise_nl_means(observed, h=0.8 * sigma_est,
patch_size=5, patch_distance=6, fast_mode=True)
score(clean, denoised, "denoised ")
restored = np.clip(richardson_lucy(denoised, psf, num_iter=20), 0, 1)
score(clean, restored, "restored ")
divider = np.ones((clean.shape[0], 4))
panel = np.hstack([clean, divider, observed, divider, denoised, divider, restored])
imsave("restoration_panel.png", img_as_ubyte(np.clip(panel, 0, 1)))
print("wrote restoration_panel.png")What's Next?
You can put the whole chapter into practice in the Hands-On Lab above, which carries one image forward through degradation and backward through the full model, prior, invert recipe with a PSNR and SSIM scorecard. This chapter completes the conceptual core of Part I: you can now describe, filter, transform, measure, and repair an image. Chapter 8: Tools of the Trade: The Image Processing Stack steps back from algorithms to engineering: how OpenCV, scikit-image, Pillow, and friends divide the territory, how to choose among them, and how to assemble the techniques of Chapters 0 through 7 into pipelines that survive production. Every restoration method you just learned has a battle-tested library implementation; Chapter 8 shows you where each one lives and what it costs.
Bibliography & Further Reading
Foundational Papers
Buades, A., Coll, B., and Morel, J.-M. "A Non-Local Algorithm for Image Denoising." CVPR (2005). IEEE Xplore
The paper that redefined a pixel's "neighborhood" as every similar patch in the image. Section 7.2's centerpiece, and the conceptual ancestor of attention mechanisms.
Dabov, K., Foi, A., Katkovnik, V., and Egiazarian, K. "Image Denoising by Sparse 3-D Transform-Domain Collaborative Filtering." IEEE TIP (2007). IEEE Xplore
BM3D: group similar patches, denoise them jointly in a transform domain. The classical state of the art for a full decade and still the baseline deep denoisers must beat.
Richardson, W. H. "Bayesian-Based Iterative Method of Image Restoration." JOSA 62(1) (1972), Optica Publishing; Lucy, L. B. "An Iterative Technique for the Rectification of Observed Distributions." The Astronomical Journal 79, 745 (1974). NASA ADS
The two independent derivations of the Richardson-Lucy deconvolution of Section 7.3: a multiplicative update that respects the Poisson statistics of photon counting. Astronomers still run it nightly, half a century later.
Bertalmio, M., Sapiro, G., Caselles, V., and Ballester, C. "Image Inpainting." SIGGRAPH (2000). ACM Digital Library
Brought the word "inpainting" from art conservation into vision, framing hole filling as a PDE that continues isophotes into the missing region (Section 7.4).
Telea, A. "An Image Inpainting Technique Based on the Fast Marching Method." Journal of Graphics Tools 9(1) (2004). Taylor & Francis
The fast, simple inpainting algorithm behind OpenCV's cv2.INPAINT_TELEA flag: fill the hole inward from the boundary, in the order fast marching dictates.
Criminisi, A., Pérez, P., and Toyama, K. "Region Filling and Object Removal by Exemplar-Based Image Inpainting." IEEE TIP (2004). IEEE Xplore
Inpainting by copying patches, with a priority rule that lets structure propagate before texture fills in. The classical method that handles large holes credibly.
Debevec, P. and Malik, J. "Recovering High Dynamic Range Radiance Maps from Photographs." SIGGRAPH (1997). ACM Digital Library
How to recover a camera's response curve and true scene radiance from a handful of bracketed exposures. The founding paper of Section 7.6's HDR pipeline.
Mertens, T., Kautz, J., and Van Reeth, F. "Exposure Fusion." Pacific Graphics (2007). IEEE Xplore
Skip radiance recovery entirely: blend the best-exposed pixels from each bracket with a Laplacian pyramid. The pragmatic algorithm inside most phone "HDR" modes.
Dong, C., Loy, C. C., He, K., and Tang, X. "Image Super-Resolution Using Deep Convolutional Networks." (2015). arXiv:1501.00092
SRCNN: the three-layer network that ended the classical era of Section 7.5 and opened the deep one. Reads as a direct translation of sparse-coding super-resolution into convolutions.
Books
Gonzalez, R. and Woods, R. "Digital Image Processing," 4th ed. Pearson (2018). Book website
Chapter 5 ("Image Restoration and Reconstruction") is the canonical textbook treatment of noise models, inverse filtering, and Wiener theory, with worked frequency-domain examples.
Szeliski, R. "Computer Vision: Algorithms and Applications," 2nd ed. Springer (2022). Free online edition
Sections on image processing, computational photography, and HDR cover this chapter's territory with modern references; free PDF from the author.
Tools & Libraries
OpenCV: Computational Photography module (photo). Official API reference
Home of fastNlMeansDenoising, inpaint, the HDR calibration and merge classes, and every tone-mapping operator used in this chapter.
scikit-image: skimage.restoration API. Official documentation
Clean NumPy implementations of Wiener and Richardson-Lucy deconvolution, non-local means, biharmonic inpainting, and noise estimation, ideal for studying the algorithms.
Recent Research (2024-2026)
Yu, F. et al. "Scaling Up to Excellence: Practicing Model Scaling for Photo-Realistic Image Restoration In the Wild" (SUPIR). CVPR (2024). arXiv:2401.13627
Restoration in the diffusion era: a generative prior so strong it reinvents missing content. The direct descendant of every prior in this chapter, with the same hallucination risks scaled up.
Wang, X. et al. "Real-ESRGAN: Training Real-World Blind Super-Resolution with Pure Synthetic Data." ICCV Workshops (2021). GitHub repository
The degradation-pipeline insight of Section 7.1 weaponized: train a restorer entirely on synthetically damaged images, and it generalizes to real damage. Actively maintained and widely deployed.