Sub-Pixel Image Processing

Contents

• Introduction and overview
• 0. Review: Fourier analysis and differential calculus
• 1. The image formation process
• 2. Shannon sampling and Fourier interpolation
• 3. Splines and lossless geometric transforms
• 4. Image reduction: aliasing, ringing and blur
• 5. Differential operators and iterative filtering
• 6. Sub-pixel image geometry
• 7. The phase spectrum: textures and sharpness
• 8. Research directions and open problems
• Exercises
• Projects and assessment
• References

Introduction and overview

Audience: a student with a good M1 in applied mathematics, computer science or physics who knows the Fourier transform and wants image algorithms that stay exact at the scale of a fraction of a pixel.

Empirical first: rotate a test image 36 times by 10° with bilinear interpolation and the full turn leaves it visibly blurred (an RMS error of 5.1 grey levels for an image of RMS 12, Code 3.11); the same full turn computed with FFT-based shears returns it to within 4 \times 10^{-8}. A striped shirt photographed with a phone shows moiré for a related reason: frequencies above the sensor’s Nyquist limit fold back as false patterns.

A digital image is an array of numbers, but the scene, the optics and most image-processing models are continuous. This course builds the bridge in both directions. It explains which continuous image the pixels are samples of, how to turn a continuous operation (a translation, a rotation, a derivative, a diffusion) into an exact or controlled discrete algorithm, and how to read geometry from pixels at a precision far below the grid step.

The tools are Fourier analysis, sampling theory, splines and a little probability. They reappear in Mathematical Models for Digital Photography, Image Denoising, Inverse Problems and Imaging, Remote Sensing Data and Generative Models for Images.

What you will be able to do

• Model the acquisition chain of a digital image — projection, diffraction, pixel integration, sampling and noise — and predict where it aliases.
• State and prove the Poisson summation formula and the Shannon–Whittaker theorem, and bound the error made when an image is not band-limited.
• Implement exact sub-pixel translations and zooms with the FFT, including the delicate Nyquist coefficient of even sizes.
• Remove the border artefacts of the DFT with the periodic plus smooth decomposition.
• Build spline interpolators by prefiltering and compare interpolators through their frequency responses.
• Rotate images losslessly with three FFT-based shears, and explain why bilinear interpolation cannot do this.
• Choose a prefilter for image reduction by quantifying the aliasing, ringing and blur it produces.
• Design consistent, isotropic discrete derivatives and stable iterative filters, from the heat equation to curvature flows.
• Estimate translations, edges and shapes with sub-pixel accuracy, and compare the result with the Cramér–Rao bound.
• Use the phase spectrum to synthesise micro-textures and to compute phase-based sharpness indices.

Official objectives and syllabus (2026–2027)

The course description on master-mva.com sets the objective of exploring ways, and applications, of linking continuous and discrete image models. It asks two questions: how to translate a continuous image-processing model into a discrete numerical algorithm, and, conversely, how to extract geometric information from a discrete array of pixels.

The same page announces sessions split half into lectures and half into exercise or computer practice, and an exam. Its two references are L. Moisan’s lecture notes Modeling and Image Processing and his 2011 paper on the periodic plus smooth image decomposition, both available on the course web page. The table matches each official item with the section of these notes that teaches it.

Beyond the syllabus, these notes add the Cramér–Rao analysis of sub-pixel registration (ch. 6.2) and a self-study chapter on research directions: Shannon-consistent variational models, aliasing in neural networks, continuous neural image representations and multi-frame super-resolution (ch. 8).

Where this course fits in the programme

The official prerequisites are the basics of Fourier analysis and differential calculus.

• Preparation — the refresher Linear Algebra, Analysis and Optimization covers the analysis (Taylor formulas, integrals, Fourier series); Scientific Python, NumPy and PyTorch covers the FFT routines used in every lab; Probability and Statistics covers the Gaussian vectors, the maximum likelihood and the Cramér–Rao bound used in ch. 6–7.
• Related courses of the same semester — Mathematical Models for Digital Photography (camera pipeline, noise, demosaicking), Image Denoising, 3D Computer Vision (pinhole geometry, homographies, stereo), The Machine Intelligence of Images and Medical Image Analysis Based on Generative, Geometric and Biophysical Models.
• Courses that build on it (semester 2) — Inverse Problems and Imaging (deconvolution, sampling models), Generative Models for Images (texture synthesis), Deep Learning for Image Restoration and Synthesis, Remote Sensing Data (satellite sampling, sub-pixel registration), Deformable Models and Geodesic Methods and Audio Signal Processing: Time–Frequency Analysis (the same Fourier tools in one dimension).

Course outline

• Review: Fourier analysis and differential calculus — the Fourier transform, the Poisson summation formula, the DFT and the curvature of level lines.
• The image formation process — pinhole geometry, diffraction, pixel integration, sampling, aliasing and noise, combined into one acquisition model.
• Shannon sampling and Fourier interpolation — the Shannon–Whittaker theorem, the interpolation of finite images by the DFT, exact translations and zooms, and the periodic plus smooth decomposition.
• Splines and lossless geometric transforms — interpolation kernels, B-splines and prefiltering, frequency responses, and rotations by three FFT shears.
• Image reduction: aliasing, ringing and blur — three measures of the defects of a prefilter, the Gibbs overshoot, and a theorem showing that no positive filter escapes the trade-off.
• Differential operators and iterative filtering — consistency and isotropy of discrete derivatives, stability of the heat equation, curvature motion.
• Sub-pixel image geometry — registration at the Cramér–Rao bound, phase correlation, peak and edge localisation, level lines.
• The phase spectrum: textures and sharpness — what the phase carries, random phase textures, and the phase-coherence sharpness index.
• Research directions and open problems — self-study: Shannon-consistent variational models, aliasing in neural networks, continuous representations, super-resolution.

Session plan

36 h in class = 12 sessions of 3 h. The volume is the official one of the Université Paris-Saclay programme table (lectures 20 h, labs 16 h); master-mva.com organises the course as half lecture sessions and half exercise or computer-practice sessions. The split below is a proposal: each session pairs about 1 h 40 of lecture with 1 h 20 of exercises or lab, which gives the official 20 h and 16 h over the twelve sessions.

Level 1 = foundations, 2 = core results and algorithms, 3 = advanced and research-level material; the level never decreases.

Self-study, outside the 36 h (about 4 h in total): ch. 8 Research directions and open problems (≈ 4 h, with Ex. 16).

How to study this module

• Read ch. 0 before the first session even if Fourier analysis feels familiar: the conventions fixed there (2\pi in the exponent, angular frequencies, unnormalised DFT) are used everywhere afterwards.
• Run every code block and change one parameter at a time; the numbers quoted in the text were produced by these blocks. They need Python 3 with NumPy, SciPy and Matplotlib (the level-line extraction of ch. 6 uses contourpy, which Matplotlib installs).
• Experiment with odd image sizes first, where the FFT tools are exactly invertible; then study the even case (ch. 2.3–2.4, Exercise 7).
• Do the exercises listed for each session before the next one; the ★★ proofs of ch. 2, 4 and 5 are typical exam questions.
• Keep two test images throughout the course: a smooth, well-sampled one and a sharp, aliased one. Most phenomena of the course appear as a difference between the two.

0. Review: Fourier analysis and differential calculus

Every chapter of the course moves between a continuous image and its samples. This review fixes the conventions of the Fourier transform, proves the one identity that links the two worlds (the Poisson summation formula), and recalls the differential geometry of level lines used in ch. 5 and ch. 6.

0.1 The Fourier transform

We use angular frequencies, so that the sampling grid \mathbb{Z}^d corresponds to the frequency band [-\pi, \pi]^d.

The rules below are used without comment in the rest of the page. Each one follows from a change of variables or an integration by parts in the definition.

Two facts matter for imaging. First, a convolution — blur, pixel integration, any linear shift-invariant filter — acts frequency by frequency as a multiplication. Second, a translation changes only the phase of the transform, never its modulus; ch. 6 and ch. 7 exploit this.

0.2 Fourier series and the Poisson summation formula

A 1-periodic function F that is regular enough is the sum of its Fourier series. Its coefficients are c_m(F) = \int_0^1 F(x)\, e^{-2 i \pi m x}\, dx, \qquad F(x) = \sum_{m \in \mathbb{Z}} c_m(F)\, e^{2 i \pi m x} .

The second form is the heart of sampling theory. The left side is computed from the samples only; the right side is the spectrum of f periodised with period 2\pi. Sampling a function therefore superposes copies of its spectrum. The identity extends to d dimensions with \mathbb{Z}^d in place of \mathbb{Z}, by applying it variable by variable.

0.3 The discrete Fourier transform

The fast Fourier transform computes the DFT in \mathcal{O}(M \log M) operations instead of \mathcal{O}(M^2). Every “exact” algorithm of ch. 2 and ch. 3 costs a few FFTs.

The index p should be read as a signed frequency. The coefficient \hat u_p with p > M/2 describes the frequency p - M, and numpy.fft.fftfreq returns these signed values.

For even M, the coefficient p = M/2 is the Nyquist frequency, shared between +M/2 and -M/2; ch. 2.3 shows that it needs special care.

0.4 Differential calculus of level lines

An image u : \mathbb{R}^2 \to \mathbb{R} of class C^2 has gradient \nabla u = (u_x, u_y) and Hessian D^2 u. Its level lines are the curves \{u = \lambda\}.

Where \nabla u \neq 0, the implicit function theorem says that the level line through a point is locally a C^2 curve, orthogonal to \nabla u.

The numerator is |\nabla u|^2 times the second derivative of u in the direction of the level line. Ch. 5.5 discretises it, and ch. 6.5 measures shapes from level lines.

Code 0.6 — Checking the conventions. The script verifies the Poisson summation formula on a Gaussian and the DFT conventions of NumPy.

import numpy as np

# Poisson summation for a Gaussian: sum_k g(k) e^{-i k xi} = sum_m ghat(xi + 2 pi m)
s = 0.7                                   # standard deviation (pixels)
g = lambda x: np.exp(-x**2 / (2 * s**2)) / np.sqrt(2 * np.pi * s**2)
ghat = lambda xi: np.exp(-s**2 * xi**2 / 2)
k, m = np.arange(-60, 61), np.arange(-20, 21)
for xi in [0.0, 1.0, np.pi]:
    lhs = np.sum(g(k) * np.exp(-1j * k * xi)).real
    rhs = np.sum(ghat(xi + 2 * np.pi * m))
    print(f"xi = {xi:.3f}:  samples side {lhs:.10f}   periodized side {rhs:.10f}")

# DFT conventions of numpy: forward unnormalized, inverse with 1/M
u = np.random.default_rng(0).normal(size=8)
U = np.fft.fft(u)
print("Parseval:", np.allclose(np.sum(u**2), np.sum(abs(U)**2) / 8))
print("Hermitian symmetry:", np.allclose(U[1:], np.conj(U[1:][::-1])))

Both sides agree to ten digits at \xi = 0, 1 and \pi.

Note the value at \xi = 0: the sum of the samples of a Gaussian of unit mass is 1.000126, not 1. The excess is the contribution of the spectral copies \hat g(\pm 2\pi), a first glimpse of aliasing.

1. The image formation process

Before processing an image we need to know what its pixels measure. This chapter follows light from the scene to the sensor: the geometry of the projection, the blur imposed by diffraction, the integration over each photosite, the sampling on a grid, and the noise. The result is an acquisition model that the rest of the course inverts or exploits.

1.1 Geometry: the pinhole camera

Homographies, and their affine approximations for small fields of view, are the geometric transforms that an image-processing chain must apply to pixels. Lens distortion adds a radial term, often modelled as x_d = x (1 + k_1 r^2 + k_2 r^4), which is estimated by calibration and then undone by resampling — again an interpolation problem (ch. 2–3).

1.2 Diffraction and the optical cut-off frequency

Light is a wave. A lens with a finite aperture cannot focus a point of the scene into a point of the sensor; it produces a small diffraction pattern instead. The remarkable consequence is that optical images are band-limited.

The first zero of J_1 is at v \approx 3.832, so the dark ring of the Airy pattern has radius 1.22\, \lambda N. A sensor samples the irradiance without aliasing only if its Nyquist frequency 1/(2p), for a pixel pitch p, exceeds \nu_c. That is, only if p \le \frac{\lambda N}{2}.

1.3 Pixel integration and the total point-spread function

With a fill factor of 100 % (c = 1), the pixel alone attenuates the Nyquist frequency \xi = \pi by \sin(\pi/2)/(\pi/2) = 2/\pi \approx 0.637. Some cameras add a birefringent optical low-pass filter in front of the sensor to attenuate aliasing further.

Because K is a convolution of several kernels, the central limit theorem suggests a Gaussian approximation K \approx g_\sigma, with \sigma measured in pixels. Its gain at the Nyquist frequency is e^{-\sigma^2 \pi^2 / 2}:

An image whose PSF has \sigma \approx 0.8 pixel is nearly free of aliasing but slightly blurred; with \sigma = 0.4, it is sharp but aliased. This tension reappears as the central trade-off of ch. 4.

1.4 Sampling and aliasing

Code 1.9 — The acquisition budget. The script evaluates the diffraction transfer function, the pixel attenuation and the folding of a pure frequency.

import numpy as np

# 1. Diffraction-limited OTF of a circular pupil (incoherent light)
lam, N = 0.55e-3, 2.8                     # wavelength (mm), f-number
nu_c = 1 / (lam * N)                      # cut-off frequency (cycles/mm)
def otf(nu):
    s = np.clip(nu / nu_c, 0, 1)
    return 2 / np.pi * (np.arccos(s) - s * np.sqrt(1 - s**2))
print(f"cut-off {nu_c:.0f} cycles/mm, first Airy zero {1.22*lam*N*1e3:.2f} um")
for p_um in [1.0, 2.0, 4.0]:              # pixel pitch (micrometres)
    nyq = 1 / (2 * p_um * 1e-3)           # Nyquist frequency (cycles/mm)
    pix = np.sinc(nyq * p_um * 1e-3)      # 100 % fill-factor pixel: sinc(1/2)
    print(f"pitch {p_um} um: Nyquist {nyq:.0f} c/mm, OTF {otf(nyq):.3f}, "
          f"pixel {pix:.3f}, total {otf(nyq)*pix:.3f}")

# 2. Aliasing of a pure frequency: 0.7 cycles/pixel is seen as 0.3
k = np.arange(64)
x = np.cos(2 * np.pi * 0.7 * k)
print("max |cos(2pi 0.7 k) - cos(2pi 0.3 k)| =", np.max(abs(x - np.cos(2*np.pi*0.3*k))))

# 3. Gaussian PSF: gain at Nyquist (xi = pi) as a function of sigma (pixels)
for sigma in [0.4, 0.6, 0.8, 1.0]:
    print(f"sigma = {sigma}: ghat(pi) = {np.exp(-sigma**2*np.pi**2/2):.3f}")

For pitches of 1, 2 and 4 µm, the product of the optical and pixel attenuations at the Nyquist frequency is 0.081, 0.332 and 0.482. Large pixels give a sharper-looking but more aliased image; small pixels approach Shannon sampling at the cost of contrast and of light per pixel.

1.5 Noise and quantisation

Quantisation with step \Delta adds an error that is close to uniform on [-\Delta/2, \Delta/2] when the signal varies over many steps. Its variance is \frac{1}{\Delta} \int_{-\Delta/2}^{\Delta/2} t^2\, dt = \frac{\Delta^2}{12} .

Noise sets the ultimate precision of every sub-pixel measurement: ch. 6.2 turns this into a lower bound, and ch. 7.3 shows how noise lowers the phase coherence of an image.

1.6 The acquisition model

The rest of the course follows this model. Interpolation (ch. 2–3) reconstructs K * u from well-sampled data. Reduction (ch. 4) produces well-sampled images from finer ones. Derivatives (ch. 5), sub-pixel geometry (ch. 6) and the phase analysis (ch. 7) are all computed on the reconstructed continuous image.

2. Shannon sampling and Fourier interpolation

When can the samples of an image determine the continuous image, and how do we compute it between the pixels? This chapter answers the question for infinite signals (the Shannon–Whittaker theorem), then for the finite images we actually store, where the DFT gives exact sub-pixel translations and zooms. The last section removes the artefacts created by the borders of a finite image.

2.1 The Shannon–Whittaker theorem

The theorem extends to two dimensions with the tensor kernel \mathrm{sinc}(x)\, \mathrm{sinc}(y) and the band [-\pi, \pi]^2. The band is sharp: the function \sin(\pi x) has frequency exactly \pi and vanishes at every integer, so the samples cannot distinguish it from zero.

The theorem has two practical weaknesses. The kernel decays only like 1/|x|, so a truncated series converges slowly. And real images are finite and not exactly band-limited. Sections 2.2 and 2.3 address both.

2.2 What happens without band-limitation

The error has two halves of equal size: the high frequencies that are lost, and the same high frequencies that come back, folded, as false low frequencies.

import numpy as np
from scipy.special import erfc

k = np.arange(-60, 61)                          # enough samples: g decays fast
x = np.linspace(-8, 8, 16001)                   # fine evaluation grid
S = np.sinc(x[:, None] - k[None, :])            # Shannon kernel sinc(x - k)
for s in [0.6, 1.0]:
    bound = 2 * erfc(s * np.pi / np.sqrt(2))    # Proposition 2.2, relative to g(0)
    for c in [0.0, 0.5]:                        # spot on a sample, or halfway
        g = lambda t: np.exp(-(t - c)**2 / (2 * s**2)) / np.sqrt(2 * np.pi * s**2)
        err = np.max(np.abs(g(x) - S @ g(k))) * np.sqrt(2 * np.pi * s**2)
        print(f"sigma = {s}, centre {c}: error {100 * err:.2f} %, bound {100 * bound:.2f} %")

2.3 Discrete Shannon interpolation of a finite image

A stored image has finitely many pixels. The natural continuous model is a trigonometric polynomial whose DFT coefficients are those of the image.

The even case is not a detail. Images from cameras have even sizes, and the choice of the cosine for the Nyquist term is the only one that makes U real, interpolating and symmetric. It has a price, shown in Proposition 2.6.

2.4 Exact translation and zoom with the FFT

The translation of item 1 is the reference sub-pixel translation of the course. It is exact for the band-limited periodic model, it can be undone, and translations compose without accumulating any blur. Item 2 is the reason why the code below works with odd sizes, or treats the Nyquist coefficient explicitly.

Code 2.7 — Shannon interpolation, translation and zoom. The script evaluates U at arbitrary points, translates by FFT, checks invertibility for odd and even sizes, and zooms by zero-padding.

import numpy as np

def shannon_eval(u, x):
    """Evaluate the real trigonometric (Shannon) interpolate of u at the points x."""
    M = len(u)
    U = np.fft.fft(u)
    p = np.fft.fftfreq(M, d=1 / M)                 # signed integer frequencies
    w = np.ones(M)
    if M % 2 == 0:                                 # split the Nyquist coefficient
        p, U, w = np.append(p, M // 2), np.append(U, U[M // 2]), np.append(w, 0.5)
        w[M // 2] = 0.5
    E = np.exp(2j * np.pi * np.outer(x, p) / M)
    return (E @ (w * U)).real / M

def shannon_shift(u, tau):
    """Exact sub-pixel translation k -> U(k - tau) (Proposition 2.6)."""
    M = len(u)
    p = np.fft.fftfreq(M, d=1 / M)
    phase = np.exp(-2j * np.pi * p * tau / M)
    if M % 2 == 0:
        phase[M // 2] = np.cos(np.pi * tau)        # Nyquist coefficient
    return np.fft.ifft(np.fft.fft(u) * phase).real

def fourier_zoom(u, z):
    """Zoom by an integer factor z (M odd): zero-padding of the centred spectrum."""
    M = len(u)
    pad = np.zeros(z * M, complex)
    c = z * M // 2 - M // 2
    pad[c:c + M] = np.fft.fftshift(np.fft.fft(u))
    return np.fft.ifft(np.fft.ifftshift(pad)).real * z

rng = np.random.default_rng(1)
for M in (9, 8):
    u = rng.normal(size=M)
    print(f"M = {M}: interpolates {np.allclose(shannon_eval(u, np.arange(M)), u)}, "
          f"shift = evaluation {np.allclose(shannon_shift(u, 0.3), shannon_eval(u, np.arange(M) - 0.3))}")
    back = shannon_shift(shannon_shift(shannon_shift(u, 0.5), 0.5), -1.0)
    print(f"        error after +1/2, +1/2, -1: {np.max(abs(back - u)):.2e}")
u = rng.normal(size=9)
print("zero-padding zoom = Shannon evaluation:",
      np.allclose(fourier_zoom(u, 4), shannon_eval(u, np.arange(36) / 4)))

For M = 9, two half-pixel translations followed by a translation of -1 return the signal to 4 \times 10^{-16}. For M = 8, the error is 0.30: the first half-pixel translation destroyed the Nyquist coefficient. Exercise 7 asks for a proof and a fix.

2.5 Borders: the periodic plus smooth decomposition

The DFT treats an image as periodic. Its left and right borders, top and bottom borders, are glued together, and the jump between them acts like a sharp edge. In the spectrum this jump creates a bright cross along the frequency axes. In interpolation it creates ringing along the borders, and in translation it makes content wrap around.

The decomposition of Moisan (2011), one of the official references of the course, splits the image into a periodic component, which carries the content and has a clean spectrum, and a smooth component, which absorbs the border jumps.

Code 2.10 — Periodic plus smooth decomposition. The image below is deliberately non-periodic: a ramp, an oblique texture and a disc.

import numpy as np

def per_smooth(u):
    """Periodic plus smooth decomposition u = p + s (Moisan, 2011)."""
    M, N = u.shape
    v = np.zeros_like(u)
    v[0, :] += u[-1, :] - u[0, :]            # jumps across the periodic boundary
    v[-1, :] += u[0, :] - u[-1, :]
    v[:, 0] += u[:, -1] - u[:, 0]
    v[:, -1] += u[:, 0] - u[:, -1]
    q = np.arange(M)[:, None]
    r = np.arange(N)[None, :]
    den = 2 * np.cos(2 * np.pi * q / M) + 2 * np.cos(2 * np.pi * r / N) - 4
    den[0, 0] = 1.0                          # avoid 0/0; the mean of s is set to 0
    s_hat = np.fft.fft2(v) / den
    s_hat[0, 0] = 0.0
    s = np.fft.ifft2(s_hat).real
    return u - s, s

M = N = 128
y, x = np.mgrid[0:M, 0:N].astype(float)
rng = np.random.default_rng(0)
u = 0.8 * x + 0.5 * y + 20 * np.sin(0.9 * x + 0.4 * y) + 60 * ((x - 70)**2 + (y - 50)**2 < 300)
u += rng.normal(0, 2, u.shape)
p, s = per_smooth(u)

def axis_energy(img):
    """Fraction of non-DC spectral energy lying on the two frequency axes."""
    E = np.abs(np.fft.fft2(img - img.mean()))**2
    return (E[0, 1:].sum() + E[1:, 0].sum()) / E.sum()

print(f"energy on the axes: u {axis_energy(u):.3f}   p {axis_energy(p):.3f}")
lap = np.roll(s, 1, 0) + np.roll(s, -1, 0) + np.roll(s, 1, 1) + np.roll(s, -1, 1) - 4 * s
print(f"max |Laplacian(s)| inside the image: {np.abs(lap[1:-1, 1:-1]).max():.1e}")
print(f"boundary jump of u: {np.abs(u[0] - u[-1]).mean():.1f}, of p: {np.abs(p[0] - p[-1]).mean():.1f}")

The share of spectral energy on the two axes falls from 0.797 for u to 0.180 for p: the cross has disappeared. The mean jump across the top and bottom borders falls from 63.4 to 3.0 grey levels. The Laplacian of s vanishes inside the image to 10^{-13}, as Theorem 2.9 predicts.

In practice, FFT-based operations are applied to p, while s — smooth and harmonic — is handled by any simple method, or added back unchanged when the operation is a filtering with small support. Two alternatives are common: mirror symmetrisation, which doubles each dimension of the image and keeps continuity but not the derivative, and apodisation by a window, which alters the content near the borders.

3. Splines and lossless geometric transforms

Shannon interpolation is exact for band-limited periodic images but global: every output value depends on every pixel. Local kernels are cheaper and handle non-periodic images, but they blur. This chapter compares interpolation kernels through their frequency responses, builds the B-spline family, and then asks which geometric transforms can be applied to an image without any loss.

3.1 Interpolation kernels

The classical kernels are nearest neighbour (\beta^0, the box), linear interpolation (\beta^1, the hat), the cubic convolution kernel of Keys (1981), cubic splines, and the sinc of ch. 2. With the parameter a = -1/2, the Keys kernel is \varphi_{\mathrm{K}}(x) = \begin{cases} \tfrac32 |x|^3 - \tfrac52 |x|^2 + 1 & |x| \le 1, \\ -\tfrac12 |x|^3 + \tfrac52 |x|^2 - 4|x| + 2 & 1 < |x| < 2, \\ 0 & |x| \ge 2 . \end{cases}

Reproducing constants is the minimum one asks of an interpolator: a flat image must stay flat. More generally, when \hat \varphi(0) \neq 0 and \hat \varphi has zeros of order L at every 2\pi m, m \neq 0, the integer shifts of the kernel span the polynomials of degree less than L.

A smooth image sampled with step h is then approximated with an error of order h^L. These are the Strang–Fix conditions, reviewed by Thévenaz, Blu and Unser (2000).

3.2 B-splines

3.3 Spline interpolation by prefiltering

The cubic B-spline does not interpolate, since \beta^3(\pm 1) \neq 0. Spline interpolation looks for coefficients c such that the spline U = \sum_k c_k \beta^3(\cdot - k) passes through the samples.

The recursions cost a few operations per sample, independently of the size, and are applied along each axis of an image. They need an initial value at each end: Unser, Aldroubi and Eden (1993) derive it for mirror boundary conditions, and with periodic conditions the filter is a simple division by B in the DFT domain.

The equivalent interpolating kernel is the cardinal cubic spline \eta^3(x) = \sum_k (b^{-1})_k\, \beta^3(x - k), \qquad \hat \eta^3(\xi) = \frac{\hat \beta^3(\xi)}{B(\xi)} .

It has infinite support but decays like 0.268^{|x|}. As the degree tends to infinity, the cardinal splines converge to the sinc kernel (Aldroubi, Unser and Eden, 1992): splines interpolate between local kernels and Shannon interpolation.

Code 3.7 — Spline interpolation, checked against SciPy. The prefilter is applied in Fourier with periodic boundary conditions; the result coincides with scipy.ndimage.shift in mode grid-wrap.

import numpy as np
from scipy import ndimage

def bspline3(x):
    x = np.abs(x)
    return np.where(x <= 1, 2/3 - x**2 + x**3 / 2, np.where(x < 2, (2 - x)**3 / 6, 0.0))

def spline_coefficients(u):
    """Periodic cubic-spline prefilter: solve b * c = u by division in Fourier."""
    M = len(u)
    B = 2/3 + (1/3) * np.cos(2 * np.pi * np.arange(M) / M)   # DFT of (1/6, 2/3, 1/6)
    return np.fft.ifft(np.fft.fft(u) / B).real

def spline_shift(u, tau):
    """Evaluate the periodic cubic-spline interpolate at the points k - tau, |tau| < 1."""
    c, M = spline_coefficients(u), len(u)
    k = np.arange(M)
    return sum(c[(k - j) % M] * bspline3(j - tau) for j in range(-2, 3))

rng = np.random.default_rng(2)
u = rng.normal(size=32)
print("interpolates:", np.allclose(spline_shift(u, 0.0), u))
print("matches scipy.ndimage.shift(order=3, mode='grid-wrap'):",
      np.allclose(spline_shift(u, 0.3), ndimage.shift(u, 0.3, order=3, mode="grid-wrap")))
delta = np.zeros(32); delta[0] = 1.0
z1 = np.sqrt(3) - 2
print("impulse response = sqrt(3) z1^|k|:",
      np.allclose(spline_coefficients(delta)[:5], np.sqrt(3) * z1**np.arange(5)))

3.4 Comparing interpolators in the frequency domain

Translating samples by \tau with a kernel \varphi is a discrete filter. Its frequency response tells which frequencies survive. For a direct kernel, the response is H_\tau(\xi) = \sum_k \varphi(k - \tau)\, e^{-i k \xi}, and for a spline it is divided by the prefilter symbol B(\xi). The Shannon translation has |H_\tau| = 1 on the whole band. The half-pixel translation \tau = 1/2 is the worst case for local kernels.

Code 3.8 — Gain of the half-pixel translation.

import numpy as np

def keys(x, a=-0.5):
    x = np.abs(x)
    return np.where(x <= 1, (a + 2) * x**3 - (a + 3) * x**2 + 1,
           np.where(x < 2, a * x**3 - 5 * a * x**2 + 8 * a * x - 4 * a, 0.0))

def bspline3(x):
    x = np.abs(x)
    return np.where(x <= 1, 2/3 - x**2 + x**3 / 2, np.where(x < 2, (2 - x)**3 / 6, 0.0))

def gain(weights_at, xi, tau=0.5, prefilter=None):
    """|frequency response| of the filter that shifts samples by tau."""
    k = np.arange(-4, 5)
    H = np.sum(weights_at(k - tau)[:, None] * np.exp(-1j * np.outer(k, xi)), axis=0)
    if prefilter is not None:
        H = H / np.sum(prefilter(k)[:, None] * np.exp(-1j * np.outer(k, xi)), axis=0)
    return np.abs(H)

xi = np.array([np.pi / 4, np.pi / 2, 3 * np.pi / 4, 0.95 * np.pi])
print("xi/pi:               ", np.round(xi / np.pi, 2))
print("linear (hat)         ", np.round(gain(lambda x: np.maximum(1 - abs(x), 0), xi), 3))
print("Keys cubic a=-1/2    ", np.round(gain(keys, xi), 3))
print("cubic B-spline interp", np.round(gain(bspline3, xi, prefilter=bspline3), 3))
print("Shannon              ", np.ones(4))

All local kernels attenuate the upper half of the band. One application costs little, but the attenuation compounds when transforms are chained: after ten half-pixel translations with linear interpolation, the gain at \xi = \pi/2 is 0.707^{10} \approx 0.031.

3.5 Lossless geometric transforms

A discrete geometric transform is lossless when it is invertible and its inverse is a transform of the same family; ideally it is unitary, so that it preserves the energy and the noise level. A chain of lossless transforms can always be undone, and iterating them accumulates no blur.

FFT translations of odd-sized images (Proposition 2.6), integer zooms by zero-padding (which have an exact left inverse, ch. 4.6), flips and rotations by multiples of 90° are lossless. Interpolation with any kernel whose gain is below 1 somewhere in the band is not. Rotations by arbitrary angles become lossless thanks to a factorisation into shears.

Three remarks keep the method honest.

• The shears are periodic, so content pushed across a border wraps around: the image must be padded so that every intermediate image fits. For a square image of half-size L and |\theta| \le \pi/4, following the corners shows that all three intermediate images fit in the bounding box of the rotated square, of half-size L(|\cos\theta| + |\sin\theta|) (Exercise 5).
• The discrete rotation is not exactly the rotation of the trigonometric polynomial of the image, because a sheared trigonometric polynomial is no longer one of the same degree. The difference is tiny for well-sampled images, as the experiment shows, and the transform is exactly invertible in any case.
• For larger angles, one first rotates by a multiple of 90°, which is exact, so that the remaining angle satisfies |\theta| \le \pi/4 and |t| \le \tan(\pi/8) \approx 0.41.

Code 3.11 — A full turn in 36 steps. The test image is a sum of smooth Gaussian spots inside a zero margin, of odd size 129.

import numpy as np
from scipy import ndimage

def shear_rows(u, a):
    """Row y is translated by a*(y - c) with the exact Fourier (Shannon) translation."""
    M, N = u.shape                                  # N odd: no Nyquist coefficient
    c = (M - 1) / 2
    p = np.fft.fftfreq(N, d=1 / N)
    tau = a * (np.arange(M) - c)
    phase = np.exp(-2j * np.pi * np.outer(tau, p) / N)
    return np.fft.ifft(np.fft.fft(u, axis=1) * phase, axis=1).real

def rotate_3shears(u, theta):
    t, s = np.tan(theta / 2), np.sin(theta)
    u = shear_rows(u, -t)
    u = shear_rows(u.T, s).T                        # vertical shear
    return shear_rows(u, -t)

n = 129
y, x = np.mgrid[0:n, 0:n] - (n - 1) / 2
rng = np.random.default_rng(3)
u = np.zeros((n, n))
for _ in range(25):
    cx, cy = rng.uniform(-30, 30, 2)
    u += rng.uniform(50, 100) * np.exp(-((x - cx)**2 + (y - cy)**2) / (2 * 2.0**2))

steps, theta = 36, 2 * np.pi / 36
v_lin, v_fft = u.copy(), u.copy()
for _ in range(steps):
    v_lin = ndimage.rotate(v_lin, np.degrees(theta), reshape=False, order=1)
    v_fft = rotate_3shears(v_fft, theta)
rms = lambda e: np.sqrt(np.mean(e**2))
print(f"after a full turn in {steps} steps: RMSE bilinear {rms(v_lin - u):.2f}, "
      f"three shears {rms(v_fft - u):.2e}  (image RMS {rms(u):.2f})")
w = rotate_3shears(rotate_3shears(u, 0.7), -0.7)
print(f"rotate by 0.7 rad then back: max error {np.abs(w - u).max():.1e}")
print("energy preserved by one rotation:", np.isclose(np.sum(rotate_3shears(u, 0.7)**2), np.sum(u**2)))

After the full turn, bilinear interpolation leaves an RMS error of 5.12 grey levels, for an image of RMS 12.04: the spots are visibly blurred. The three-shear rotation returns the image to within 3.9 \times 10^{-8}. A rotation by 0.7 radian followed by its inverse is exact to 10^{-13}, and one rotation preserves the energy exactly. This is the experiment announced in the introduction.

4. Image reduction: aliasing, ringing and blur

Reducing an image by a factor z — for a thumbnail, a pyramid, or a coarser analysis scale — means sampling it on a grid z times coarser. The new grid can carry only the frequencies below \pi / z. Whatever we do with the frequencies above, we pay in aliasing, in ringing or in blur. This chapter measures the three defects and proves that no filter without ringing can avoid trading blur against aliasing.

4.1 The reduction model

On the band [-\pi, \pi]^d of the coarse grid, the term m = 0 carries the frequencies of U below \pi / z, attenuated by \hat h.

The terms m \neq 0 carry frequencies above \pi / z, folded. An ideal prefilter would have \hat h = 1 below \pi / z and \hat h = 0 above.

4.2 Three measures of the defects

The three numbers are dimensionless and lie in [0, 1] for reasonable filters. A filter with h \ge 0 has a non-decreasing step response with values in [0, 1], so it has no ringing at all.

4.3 The ideal low-pass filter and the Gibbs overshoot

This is the Gibbs phenomenon. Reducing the factor z narrows the oscillations but never lowers them: every sharp edge of an image reduced with the ideal filter is surrounded by ripples of 9 % of its contrast. Exercise 10 adds a second obstruction: a filter with compact support cannot have \hat h = 0 on a whole interval, so zero aliasing needs an infinitely extended kernel.

4.4 No escape without ringing

The constant 4 cannot be improved: for a Gaussian, the ratio of the two sides tends to 1 as \xi \to 0. The theorem is the mathematical content of the trade-off of the syllabus. A filter without ringing that keeps the band within \varepsilon passes the next octave within 4\varepsilon, which then folds. To remove aliasing one must either blur — a wide Gaussian — or accept negative lobes and therefore ringing.

4.5 The trade-off in numbers

Code 4.6 — Measuring prefilters for z = 2. The integrals of Definition 4.3 are computed on a fine frequency grid, and the ringing of the ideal filter by the formula of Proposition 4.4.

import numpy as np
from scipy.special import sici

z = 2.0                                          # reduction factor
xi = np.linspace(-80, 80, 800001)                # angular frequencies
x = np.linspace(-40, 40, 160001)                 # positions, in fine pixels

def aliasing_blur(H):
    E = np.abs(H)**2
    band = np.abs(xi) <= np.pi / z
    alias = E[~band].sum() / E.sum()             # energy that will fold back
    blur = (z / (2 * np.pi)) * np.sum(1 - E[band]) * (xi[1] - xi[0])
    return alias, blur

rows = []
for c in [0.3, 0.5, 0.8, 1.0]:                   # Gaussian, sigma = c z
    rows.append((f"Gaussian sigma = {c} z", np.exp(-(c * z)**2 * xi**2 / 2), 0.0))
rows.append(("box of width z", np.sinc(z * xi / (2 * np.pi)), 0.0))
step = 0.5 + sici(np.pi * x / z)[0] / np.pi      # step response of the ideal filter
rows.append(("ideal low-pass", (np.abs(xi) <= np.pi / z) * 1.0, step.max() - 1))
print(f"{'prefilter':24s} aliasing   blur  ringing")
for name, H, ring in rows:
    a, b = aliasing_blur(H)
    print(f"{name:24s} {a:8.4f} {b:6.3f} {ring:8.4f}")

Each row sits on one side of the triangle. The Gaussian family slides from aliasing to blur as \sigma grows; the common choice \sigma \approx 0.8 z makes aliasing negligible but blurs noticeably. Averaging z \times z blocks — what binning sensors and many resizing routines do — is the worst compromise: it blurs as much as a Gaussian of width 0.3z and aliases more. The ideal filter removes both but rings. Windowed sinc kernels, such as the Lanczos kernels, sit in between, with small negative lobes and small ringing.

4.6 Reduction algorithms in practice

• Fourier cropping keeps the central frequencies of the DFT: the ideal low-pass filter on the periodic model. It has no aliasing but rings, and it should be applied to the periodic component of ch. 2.5 so that the border jump does not ring as well.
• Gaussian filtering then subsampling with \sigma \approx 0.8 z has no ringing, almost no aliasing and visible blur. It is the standard choice for scale-space pyramids.
• Block averaging is cheap but aliases, as the table shows.
• Least-squares spline reduction projects the spline model of the fine image on the coarser spline space (Unser, Aldroubi and Eden, 1995); it sits between the Gaussian and the ideal filter.

Code 4.8 — Reducing a zone plate. A zone plate \cos(\alpha r^2) has a local frequency 2\alpha r that grows with the radius. Reduced by z = 2, its outer ring, with local frequency above the new Nyquist limit, should disappear, while its centre, well inside the band, should keep its contrast (RMS 0.707).

import numpy as np
from scipy import ndimage

n, z = 256, 2
y, x = np.mgrid[0:n, 0:n] - n / 2
r = np.hypot(x, y)
alpha = np.pi / n                             # local frequency 2*alpha*r = pi at r = n/2
u = np.cos(alpha * r**2) * (r < 0.48 * n)     # zone plate, cut before the fine Nyquist

def fourier_reduce(u, z):
    """Keep the central (n/z)^2 DFT coefficients (ideal low-pass + subsampling)."""
    m = u.shape[0] // z
    U = np.fft.fftshift(np.fft.fft2(u))
    c = u.shape[0] // 2
    V = U[c - m // 2:c + m // 2, c - m // 2:c + m // 2].copy()
    V[0, :] = 0; V[:, 0] = 0                  # drop the (ambiguous) Nyquist row/column
    return np.fft.ifft2(np.fft.ifftshift(V)).real / z**2

out = {
    "naive subsampling": u[::z, ::z],
    "Gaussian sigma=0.8z": ndimage.gaussian_filter(u, 0.8 * z, mode="wrap")[::z, ::z],
    "Fourier cropping": fourier_reduce(u, z),
}
rc = r[::z, ::z]                              # radius of the coarse samples (fine units)
outer = (rc > 80) & (rc < 115)                # local frequency above the coarse Nyquist
inner = rc < 40                               # local frequency well inside the band
for name, v in out.items():
    print(f"{name:20s} RMS outer ring {np.sqrt(np.mean(v[outer]**2)):.3f}   "
          f"RMS inner disc {np.sqrt(np.mean(v[inner]**2)):.3f}")

The ringing of Fourier cropping is especially visible here. For a chirp, frequency and position are tied, so cutting the spectrum at \pi / z cuts the pattern sharply at the radius where its local frequency crosses that value, r = 64. The sinc tail of that cut oscillates across the outer ring.

5. Differential operators and iterative filtering

Edge detectors, curvature estimates and diffusion filters all need derivatives of an image that only exists on a grid. This chapter measures the quality of a discrete derivative by its symbol, compares finite differences with spectral derivatives, designs isotropic stencils, and studies the stability of filters applied many times in a row — linear (the heat equation) and nonlinear (curvature motion, Perona–Malik).

5.1 Symbols and consistency

On a well-sampled image, the relative error of L at frequency \xi is |\hat L(\xi) - P(i\xi)| / |P(i\xi)|. It is small at low frequencies and large near the Nyquist frequency.

At \xi_1 = \pi/2, the centred difference returns \sin(\pi/2) = 1 instead of \pi/2 \approx 1.571, an error of 36 %. At \xi_1 = \pi it returns 0: the checkerboard pattern is invisible to it.

5.2 Spectral derivatives

Spectral derivatives are exact for the continuous model but inherit its weaknesses. On a non-periodic image, the border jump makes the derivative ring along the borders; one differentiates the periodic component p of ch. 2.5 and the smooth component s separately. And differentiation multiplies the noise at frequency \xi by |\xi|. In practice one differentiates a slightly smoothed image, for example a Gaussian of width \sigma \approx 1, which is the derivative-of-Gaussian filter.

5.3 Isotropic gradients and Laplacians

The norm of the gradient should not depend on the orientation of an edge. On the plane wave of frequency \xi, the exact gradient has norm |\xi|. A stencil is isotropic to some order if its symbol depends only on |\xi| up to that order.

The Sobel operator uses w = 2; the weights (3, 10, 3) used by the Scharr operator in OpenCV correspond to w = 10/3.

Code 5.4 — Isotropy of gradient stencils. For each stencil, the gain |\hat \nabla(\xi)| / |\xi| is computed over all directions, at two radii.

import numpy as np

def grad_gain(w, r, theta):
    """|symbol of the smoothed centred-difference gradient| / r on the plane wave
    of angular frequency r in direction theta; smoothing (1, w, 1)/(w + 2)."""
    x1, x2 = r * np.cos(theta), r * np.sin(theta)
    sm = lambda t: (w + 2 * np.cos(t)) / (w + 2)
    gx, gy = np.sin(x1) * sm(x2), np.sin(x2) * sm(x1)
    return np.hypot(gx, gy) / r

theta = np.linspace(0, np.pi / 4, 91)
print(f"{'smoothing weights':22s}  r = pi/4: min/max   r = pi/2: min/max")
for name, w in [("none (1D centred)", np.inf), ("Sobel (1,2,1)", 2.0),
                ("Scharr (3,10,3)", 10 / 3), ("(1,4,1)", 4.0)]:
    out = []
    for r in [np.pi / 4, np.pi / 2]:
        if np.isinf(w):
            g = np.hypot(np.sin(r * np.cos(theta)), np.sin(r * np.sin(theta))) / r
        else:
            g = grad_gain(w, r, theta)
        out.append(f"{g.min():.3f}/{g.max():.3f}")
    print(f"{name:22s}  {out[0]:>17s}   {out[1]:>17s}")

The (1, 4, 1) weights are the most isotropic at low frequency, as Proposition 5.3 predicts. The Scharr weights, tuned on a wider band, are the most isotropic at \pi/2. All of them underestimate the gradient at mid frequencies, a consequence of the centred difference itself; a spectral or spline-based derivative removes this bias.

The same analysis applies to the Laplacian.

5.4 Iterative linear filtering: the heat equation

The simplest iterative filter repeats a small averaging. Its continuous model is the heat equation \partial_t u = \Delta u, whose solution at time t is the convolution with the Gaussian of standard deviation \sigma = \sqrt{2t}. The explicit scheme with time step \delta t reads u^{n+1} = u^n + \delta t\, \Delta_5 u^n .

Two exact alternatives avoid the time step altogether. The semi-discrete equation \partial_t u = \Delta_5 u is solved exactly by multiplying the DFT by e^{t \hat L_5(\xi)}, a “discrete Gaussian”. The continuous heat equation on the Shannon interpolate is solved by the multiplier e^{-t|\xi|^2}.

Code 5.7 — Stability in practice. The explicit scheme is run up to t = 2 (\sigma = 2 pixels) on white noise, and compared with the exact semi-discrete solution.

import numpy as np

def lap5(u):
    return (np.roll(u, 1, 0) + np.roll(u, -1, 0) + np.roll(u, 1, 1)
            + np.roll(u, -1, 1) - 4 * u)

def heat_explicit(u, t, dt):
    for _ in range(int(round(t / dt))):
        u = u + dt * lap5(u)
    return u

def heat_exact_discrete(u, t):
    """Exact semigroup of the 5-point Laplacian: multiply the DFT by exp(t * symbol)."""
    M, N = u.shape
    q = 2 * np.pi * np.fft.fftfreq(M)[:, None]
    r = 2 * np.pi * np.fft.fftfreq(N)[None, :]
    sym = 2 * np.cos(q) + 2 * np.cos(r) - 4
    return np.fft.ifft2(np.fft.fft2(u) * np.exp(t * sym)).real

def heat_gaussian(u, t):
    """Continuous heat equation on the Shannon interpolate: exp(-t |xi|^2)."""
    M, N = u.shape
    q = 2 * np.pi * np.fft.fftfreq(M)[:, None]
    r = 2 * np.pi * np.fft.fftfreq(N)[None, :]
    return np.fft.ifft2(np.fft.fft2(u) * np.exp(-t * (q**2 + r**2))).real

rng = np.random.default_rng(4)
u = rng.normal(size=(65, 65))
t = 2.0                                          # sigma = sqrt(2t) = 2 pixels
ref = heat_exact_discrete(u, t)
for dt in [0.05, 0.125, 0.25, 0.30]:
    v = heat_explicit(u, t, dt)
    print(f"dt = {dt:5.3f}: max |explicit - exact discrete| = {np.abs(v - ref).max():.2e}")
print(f"discrete vs continuous Gaussian semigroup: {np.abs(ref - heat_gaussian(u, t)).max():.2e}")

With \delta t = 0.05 and 0.125, the explicit scheme is within 0.009 and 0.023 of the exact solution: the error of the time discretisation grows with \delta t.

At \delta t = 1/4 it is stable but wrong by 0.46: the checkerboard has g = -1 and is never damped. At \delta t = 0.3 it explodes (error 5.1 after only seven steps). Stability is necessary, not sufficient.

One explicit heat step has kernel weights 1 - 4\,\delta t at the centre and \delta t at the four neighbours, hence s^2 = 2\,\delta t per axis. After n steps, s^2 n = 2 t, as predicted by the heat equation. The Gaussian is the universal limit of iterated linear smoothing, which is why it dominates scale-space theory (Koenderink, 1984; Lindeberg, 1994).

5.5 Nonlinear iterative filtering: curvature motion and Perona–Malik

Linear smoothing blurs edges. Nonlinear filters can smooth along edges without blurring across them.

Contrast invariance means that the filter acts on the level lines, not on the grey levels: changing the exposure or the gamma of the image commutes with the filtering. Alvarez, Guichard, Lions and Morel (1993) show that, together with a few other axioms, it singles out the curvature motions and their affine-invariant variant, \partial_t u = |\nabla u|\, \mathrm{curv}(u)^{1/3}.

The discretisation uses the formula of Definition 5.9 with finite differences — preferably the isotropic stencils of ch. 5.3 — and a small \varepsilon^2 added to the denominator to avoid division by zero. Explicit steps with \delta t \approx 0.1 are common. Near flat regions the equation degenerates, and Moisan (1998) proposes a scheme for the affine variant that is consistent in all cases.

Backward diffusion sharpens edges whose slope exceeds K, which is the purpose of the filter (Perona and Malik, 1990). It also makes the continuous problem ill-posed; discrete schemes behave well thanks to the implicit regularisation of the grid, at the price of “staircasing” of smooth ramps.

6. Sub-pixel image geometry

A pixel is not the limit of geometric precision. When an image is well sampled, its continuous model is known everywhere, and positions, displacements and shapes can be measured to a small fraction of the grid step. The limit is set by noise and by model errors such as aliasing. This chapter derives that limit and the estimators that reach it.

6.1 Why sub-pixel accuracy matters

The same arithmetic drives satellite imaging, where small-baseline stereo needs disparities accurate to a few hundredths of a pixel (Delon and Rougé, 2007), and also image stabilisation, camera calibration, motion estimation and the measurement of displacement fields in mechanics.

6.2 Registration and the Cramér–Rao bound

The Gauss–Newton method linearises the residual around the current estimate. Since \partial_\tau U(k - \tau) = -U'(k - \tau), one step reads \tau \leftarrow \tau - \frac{\sum_k \big( v(k) - U(k - \tau) \big)\, U'(k - \tau)}{\sum_k U'(k - \tau)^2} .

In two dimensions, with the gradient in place of U' and a 2 \times 2 system, this is the Lucas–Kanade algorithm (Lucas and Kanade, 1981). The translations and derivatives are computed exactly by the FFT tools of ch. 2 and ch. 5.

The bound says that precision grows with contrast and with the number of edges, and decreases with noise. When the reference image is noisy too, the variance roughly doubles. In two dimensions, a single straight edge gives a singular structure matrix F: the displacement along the edge cannot be measured. This is the aperture problem, studied in Exercise 12. Robinson and Milanfar (2004) analyse these limits for image registration in detail.

6.3 Phase correlation

The phase correlation method uses the fact that a translation acts only on the phase of the DFT (ch. 0.1).

This two-sample formula is due to Foroosh, Zerubia and Berthod (2002). Guizar-Sicairos, Thurman and Fienup (2008) locate the peak instead by evaluating the cross-correlation on a finely upsampled grid around it, with small matrix-multiplication DFTs.

Phase correlation is insensitive to changes of contrast and finds large displacements without any initial guess. But the normalisation gives the same weight to every frequency, including those where the signal is weaker than the noise. It is therefore a good initialisation and a poor final estimator.

Code 6.6 — Phase correlation, Gauss–Newton and the Cramér–Rao bound. A smooth signal of length 129 is translated by 0.37 pixel, and noise of standard deviation 2 is added; 400 trials are run.

import numpy as np

M = 129                                           # odd length: no Nyquist issue
p = np.fft.fftfreq(M, d=1 / M)

def shift(u, tau):
    return np.fft.ifft(np.fft.fft(u) * np.exp(-2j * np.pi * p * tau / M)).real

def derivative(u):
    return np.fft.ifft(np.fft.fft(u) * 2j * np.pi * p / M).real

def phase_correlation(u, v):
    """Integer peak of the phase correlation, then Foroosh's two-sample formula."""
    R = np.fft.fft(v) * np.conj(np.fft.fft(u))
    r = np.fft.ifft(R / np.maximum(np.abs(R), 1e-12)).real
    k0 = int(np.argmax(r))
    k1 = (k0 + 1) % M if r[(k0 + 1) % M] > r[k0 - 1] else (k0 - 1) % M
    d = r[k1] / (r[k0] + r[k1])                   # fraction of the way from k0 to k1
    t = k0 + d if k1 == (k0 + 1) % M else k0 - d
    return (t + M / 2) % M - M / 2

def gauss_newton(u, v, tau0=0.0, iters=5):
    """Minimise sum (v - u(. - tau))^2 over tau with Shannon interpolation."""
    tau = tau0
    for _ in range(iters):
        w = shift(u, tau)
        dw = derivative(w)                        # d/dtau u(k - tau) = -u'(k - tau)
        tau += np.sum((v - w) * (-dw)) / np.sum(dw**2)
    return tau

# a well-sampled smooth signal: Gaussian-blurred random steps
rng = np.random.default_rng(5)
u0 = np.repeat(rng.uniform(0, 100, 16), 9)[:M]
u = np.fft.ifft(np.fft.fft(u0) * np.exp(-(2 * np.pi * p / M)**2 * 1.5**2 / 2)).real
tau, sigma, trials = 0.37, 2.0, 400
est_pc, est_gn = [], []
for _ in range(trials):
    v = shift(u, tau) + rng.normal(0, sigma, M)
    t0 = phase_correlation(u, v)
    est_pc.append(t0)
    est_gn.append(gauss_newton(u, v, t0))
crb = sigma / np.sqrt(np.sum(derivative(u)**2))
print(f"true shift {tau}")
print(f"phase correlation: mean {np.mean(est_pc):.4f}, std {np.std(est_pc):.4f}")
print(f"Gauss-Newton:      mean {np.mean(est_gn):.4f}, std {np.std(est_gn):.4f}")
print(f"Cramer-Rao bound on the std: {crb:.4f}")
v = shift(u, tau)
print(f"noise-free: phase correlation {phase_correlation(u, v):.4f}, "
      f"Gauss-Newton {gauss_newton(u, v, 0.0):.6f}")

Without noise, phase correlation returns 0.3714 (the small bias comes from the finite length, s_M \neq \mathrm{sinc}) and Gauss–Newton returns 0.370000. With noise, Gauss–Newton reaches the Cramér–Rao bound, while phase correlation is ten times less precise. The good practice is the combination: phase correlation to find the integer part, least squares to reach the bound.

6.4 Peaks and edges

Many sub-pixel methods locate the maximum of a function known only at integer points: a correlation peak, a gradient norm across an edge, a spot in a fluorescence image.

The parabola is exact only for parabolic peaks. For a Gaussian peak of width one pixel, its bias reaches 0.048 pixel, and for a sinc peak 0.119 pixel (Exercise 13). The estimates are pulled towards integer positions, a phenomenon called pixel locking in particle image velocimetry. Fitting the parabola to the logarithm of the values is exact for Gaussian peaks.

Devernay (1995) makes Canny’s detector sub-pixel with Proposition 6.7. At each pixel where the gradient norm is maximal along the gradient direction, he fits a parabola through the three norms along that direction and takes the vertex as the edge position. On a well-sampled image, one can instead search for the zero of U_{\eta\eta} along the gradient line on the Shannon interpolate, with a few Newton steps.

6.5 Level lines and shape measurement

A level line of the continuous image is an exact sub-pixel object: its points are defined at any precision. The simplest continuous model that can be computed cell by cell is the bilinear interpolation.

The marching-squares algorithm joins these edge crossings, cell by cell, into polygons. When a cell has four crossings, the saddle value decides how they connect: the two branches of the hyperbola are separated by the lines x = x_s and y = y_s, and the sign of \lambda - u_s tells which pairs of crossings belong to the same branch. Monasse and Guichard (2000) compute the tree of the level sets of an image, and Caselles and Monasse (2010) describe its version for the bilinear interpolate.

This is the curvature motion of ch. 5.5 in disguise: by Example 5.11, a disc evolved for a time t = \sigma^2/2 has radius \sqrt{R^2 - \sigma^2} \approx R - \sigma^2/(2R).

Code 6.11 — Measuring a disc to a hundredth of a pixel. A disc of radius 10.3 is integrated over each pixel, blurred by a Gaussian of width 0.6 pixel, and corrupted by a small noise.

import numpy as np
import contourpy

# A disc of radius R, integrated over each pixel (16 x 16 supersampling), then blurred
n, R, c0 = 64, 10.3, (31.6, 32.2)
ss = 16
f = (np.arange(n * ss) + 0.5) / ss - 0.5
yy, xx = np.meshgrid(f, f, indexing="ij")
fine = ((xx - c0[0])**2 + (yy - c0[1])**2 < R**2).astype(float)
u = fine.reshape(n, ss, n, ss).mean(axis=(1, 3))          # pixel integration
q = 2 * np.pi * np.fft.fftfreq(n)
u = np.fft.ifft2(np.fft.fft2(u) * np.exp(-0.6**2 * (q[:, None]**2 + q[None, :]**2) / 2)).real
rng = np.random.default_rng(6)
u = u + rng.normal(0, 0.01, u.shape)

# 1. pixel counting: area of {u > 1/2}
R_count = np.sqrt(np.sum(u > 0.5) / np.pi)
# 2. level line of the bilinear interpolate (marching squares), polygon area
line = contourpy.contour_generator(z=u).lines(0.5)[0]       # (x, y) vertices
xs, ys = line[:, 0], line[:, 1]
area = 0.5 * abs(np.dot(xs, np.roll(ys, 1)) - np.dot(ys, np.roll(xs, 1)))
R_line = np.sqrt(area / np.pi)
# 3. total mass (the integral of a disc of height 1 is pi R^2)
R_mass = np.sqrt(u.sum() / np.pi)
print(f"true radius {R}")
print(f"pixel counting     {R_count:.3f}")
print(f"bilinear level line {R_line:.3f}")
print(f"integral of u       {R_mass:.3f}")
cross = xs * np.roll(ys, -1) - np.roll(xs, -1) * ys        # shoelace terms
A = cross.sum() / 2
cx = np.sum((xs + np.roll(xs, -1)) * cross) / (6 * A)
cy = np.sum((ys + np.roll(ys, -1)) * cross) / (6 * A)
print(f"polygon centroid ({cx:.3f}, {cy:.3f}), true {c0}")
s2 = 0.6**2 + 1 / 12                                        # blur + pixel variance
print(f"curvature shift sigma^2/(2R) = {s2 / (2 * R):.3f} pixel")

The level line is about twice as accurate as pixel counting: its error is 0.031 pixel, against 0.067.

Most of its remaining error is the curvature shift of Proposition 6.10: the total variance is 0.6^2 + 1/12 (blur plus pixel integration), which predicts a shift of 0.022 pixel, hence a radius of 10.278. The last 0.009 pixel comes from the polygon being inscribed in the curve, and from the noise. The centroid of the polygon, (31.595, 32.205), is within 0.005 pixel of the true centre.

7. The phase spectrum: textures and sharpness

The DFT of an image has a modulus and a phase. The modulus is easy to interpret — how much energy at each frequency — and the phase looks like noise. Yet the phase carries the geometry: positions, edges, the alignment of structures. This chapter makes the statement precise and turns it into two tools of the syllabus: texture synthesis by random phases, and sharpness metrics based on phase coherence.

7.1 Phase versus modulus

Code 7.1 — Swapping modulus and phase. Two images made of different shapes exchange their DFT moduli; the result is compared with both sources.

import numpy as np

n = 128
y, x = np.mgrid[0:n, 0:n]
A = 100.0 * ((x - 40)**2 + (y - 50)**2 < 400) + 60.0 * ((abs(x - 90) < 15) & (abs(y - 80) < 25))
B = 80.0 * ((x - 85)**2 + (y - 35)**2 < 250) + 70.0 * ((abs(x - 35) < 20) & (abs(y - 95) < 12))
FA, FB = np.fft.fft2(A), np.fft.fft2(B)
mag_A_phase_B = np.fft.ifft2(np.abs(FA) * np.exp(1j * np.angle(FB))).real
corr = lambda s, t: np.corrcoef(s.ravel(), t.ravel())[0, 1]
print(f"|A| with phase of B: correlation with A {corr(mag_A_phase_B, A):.2f}, "
      f"with B {corr(mag_A_phase_B, B):.2f}")

The image built from the modulus of A and the phase of B looks like B: its correlation is -0.02 with the first image and 0.78 with the second.

This is the classical experiment of Oppenheim and Lim (1981). The next proposition explains why the modulus alone cannot carry positions.

An edge, in the Fourier view, is a set of frequencies whose phases align at the position of the edge: all the sinusoids reach an extremum or a zero crossing at the same place. Kovesi (1999) builds an edge detector, phase congruency, on this observation. Natural images have moduli that decay roughly like a power of |\xi|, much alike from one image to the next; what distinguishes them is the phase.

7.2 Random phase textures

If the phase carries the geometry, an image whose phase is replaced by noise should keep the “texture” of the original and lose its structures. This gives a remarkably simple texture synthesis method for textures without structure, called micro-textures.

The phase of the DFT of a real Gaussian white noise has exactly the distribution of Definition 7.3, which gives the simplest way of sampling \theta.

RPN and ADSN differ only by the modulus: RPN keeps |\hat t| exactly, while ADSN multiplies it by a random Rayleigh factor. Galerne, Gousseau and Morel (2011) study both, show that they synthesise convincing micro-textures (sand, granite, fabric seen from afar), and apply them to the periodic component of ch. 2.5 to avoid the cross artefacts of the border jump. Both fail on textures made of distinct objects, whose structure lives in the phase.

Code 7.6 — Random phase synthesis. The script checks that RPN preserves the DFT modulus exactly, and uses the kurtosis of the grey-level histogram (3 for a Gaussian) as a crude witness of structure.

import numpy as np

rng = np.random.default_rng(7)
n = 128
y, x = np.mgrid[0:n, 0:n]

def random_phase_noise(u, rng):
    """RPN: keep |DFT|, draw a uniform random phase with Hermitian symmetry."""
    m = u.mean()
    theta = np.angle(np.fft.fft2(rng.normal(size=u.shape)))  # odd and uniform phase field
    return m + np.fft.ifft2(np.abs(np.fft.fft2(u - m)) * np.exp(1j * theta)).real

def kurtosis(v):
    v = v - v.mean()
    return np.mean(v**4) / np.mean(v**2)**2

# micro-texture: oriented filtered noise ; structured texture: random discs
W = rng.normal(size=(n, n))
q = np.fft.fftfreq(n)[:, None]; r = np.fft.fftfreq(n)[None, :]
micro = np.fft.ifft2(np.fft.fft2(W) * np.exp(-((q * 0.8 - r * 0.6) * 60)**2
                                             - ((q * 0.6 + r * 0.8) * 12)**2)).real
discs = np.zeros((n, n))
for cx, cy in rng.uniform(0, n, (40, 2)):                     # 40 random discs of radius 4
    discs[((x - cx + n / 2) % n - n / 2)**2 + ((y - cy + n / 2) % n - n / 2)**2 < 16] = 1.0
for name, u in [("micro-texture", micro), ("random discs", discs)]:
    v = random_phase_noise(u, rng)
    same = np.allclose(np.abs(np.fft.fft2(v - v.mean())), np.abs(np.fft.fft2(u - u.mean())))
    print(f"{name:14s}: same |DFT| {same}, kurtosis {kurtosis(u):.2f} -> {kurtosis(v):.2f}")

Both syntheses keep the DFT modulus exactly. The micro-texture, already Gaussian-like (kurtosis 2.92), stays so (3.15): the synthesis is statistically faithful. The binary discs (kurtosis 6.50) become a Gaussian-looking field (2.91): their shapes, which were in the phase, are gone.

7.3 Phase coherence and sharpness metrics

A sharp image has strongly aligned phases. Blur and noise both weaken this alignment. Blanchet, Moisan and Rougé (2008) turned this into a measure by asking how unlikely the image would be if its phases were random.

Random phases spread the energy of edges over the whole image and increase the total variation. A sharp image therefore has a much smaller TV than its random-phase versions, and a large GPC. The probability must be estimated by Monte Carlo, which cannot reach values like 10^{-300}. Blanchet and Moisan (2012) replace RPN by ADSN and approximate the distribution of the TV by a Gaussian, whose two moments are explicit.

Code 7.10 — The sharpness index. The script computes \mu and \nu by Proposition 7.9, checks them by Monte Carlo on 2000 ADSN samples, then evaluates SI on a sharp synthetic image, two blurred versions and a noisy one.

import numpy as np
from scipy.special import log_ndtr

def dx(u): return np.roll(u, -1, axis=1) - u          # periodic differences
def dy(u): return np.roll(u, -1, axis=0) - u
def tv(u): return np.abs(dx(u)).sum() + np.abs(dy(u)).sum()

def tv_moments_adsn(u):
    """Mean and variance of TV(ADSN(u)) (Proposition 7.9)."""
    MN = u.size
    g = {"x": dx(u - u.mean()), "y": dy(u - u.mean())}
    F = {a: np.fft.fft2(g[a]) for a in g}
    sd = {a: np.sqrt(np.sum(g[a]**2) / MN) for a in g}
    mu = MN * np.sqrt(2 / np.pi) * (sd["x"] + sd["y"])
    omega = lambda t: t * np.arcsin(t) + np.sqrt(1 - t**2) - 1
    var = 0.0
    for a in g:
        for b in g:
            gamma = np.fft.ifft2(np.conj(F[a]) * F[b]).real      # cross-correlation
            rho = np.clip(gamma / (MN * sd[a] * sd[b]), -1, 1)
            var += MN * (2 / np.pi) * sd[a] * sd[b] * omega(rho).sum()
    return mu, var

def sharpness_index(u):
    mu, var = tv_moments_adsn(u)
    return -log_ndtr(-(mu - tv(u)) / np.sqrt(var)) / np.log(10)   # -log10 P(N >= t)

rng = np.random.default_rng(8)
n = 64
y, x = np.mgrid[0:n, 0:n]
u = 100.0 * ((x - 30)**2 + (y - 34)**2 < 150) + 50.0 * ((abs(x - 45) < 8) & (abs(y - 15) < 6))
u = u + 30.0 * (x > 52)

mu, var = tv_moments_adsn(u)
t = u - u.mean()
samples = [tv(np.fft.ifft2(np.fft.fft2(t) * np.fft.fft2(rng.normal(size=t.shape))).real / n)
           for _ in range(2000)]
print(f"mean TV of ADSN: formula {mu:.0f}, Monte Carlo {np.mean(samples):.0f}")
print(f"std  TV of ADSN: formula {np.sqrt(var):.0f}, Monte Carlo {np.std(samples):.0f}")

q = 2 * np.pi * np.fft.fftfreq(n)
def blur(v, s):
    return np.fft.ifft2(np.fft.fft2(v) * np.exp(-s**2 * (q[:, None]**2 + q[None, :]**2) / 2)).real
for name, v in [("sharp", u), ("blur sigma=1", blur(u, 1.0)), ("blur sigma=2", blur(u, 2.0)),
                ("noise std 20", u + rng.normal(0, 20, u.shape))]:
    print(f"{name:13s} SI = {sharpness_index(v):7.1f}")

The formulas agree with the simulation: mean 79 484 against 79 459, standard deviation 1707 against 1678 (the Monte Carlo standard error on the latter is about 27).

An index of 300 means that, in the Gaussian approximation of Definition 7.8, the ADSN of the image — same autocorrelation, random phases — would have a total variation as small as the observed one with probability 10^{-300}. A blur of one pixel divides the index by 15, and a strong noise almost cancels it. The index needs no reference image, which makes it useful for autofocus, quality assessment and the tuning of restoration methods.

Leclaire and Moisan (2015) analyse these phase-based indices, propose a simplified variant, and use them for no-reference quality assessment and blind deblurring: among candidate deblurring kernels, the one that maximises the index is selected. Exercise 15 reproduces a small version of this experiment. Related ideas measure phase coherence locally, in a wavelet decomposition (Wang and Simoncelli, 2003).

8. Research directions and open problems

This chapter is for self-study. It points to current research that extends the questions of the course — how to discretise continuous models consistently, and how to measure geometry below the pixel — with the main references. The descriptions are short; the papers are the place to go further.

8.1 Shannon-consistent variational models

Many restoration methods minimise an energy that involves the total variation. Its usual discretisation by finite differences is not invariant under sub-pixel translations or rotations, and its minimisers show pixel-scale artefacts that make them hard to interpolate afterwards. Abergel and Moisan (2017) define a Shannon total variation, the total variation of the Shannon interpolate estimated by oversampling the gradient, and show that its minimisers are well interpolable and free of these artefacts. Condat (2017) proposes another discrete total variation with better isotropy.

The open question is general: for a non-smooth functional of a continuous image, which discretisation commutes with the geometric transforms of ch. 3 and keeps the result well sampled?

8.2 Aliasing inside neural networks

Convolutional networks subsample their feature maps with strides and pooling, without the prefilter that ch. 4 requires. Azulay and Weiss (2019) show that their outputs can change markedly when the input moves by one pixel. Zhang (2019) inserts a small blur before each subsampling (“anti-aliased” or “blur-pool” networks) and improves both shift-invariance and accuracy: Theorem 4.5 in a modern setting.

Karras and coauthors (2021) go further for generative networks. They treat feature maps as samples of continuous band-limited signals, apply carefully designed low-pass filters around every nonlinearity, and obtain generated images whose details move coherently with the content instead of sticking to the pixel grid. Pointwise nonlinearities create frequencies beyond the band, which is the mechanism behind Exercise 16.

Open problems include the best compromise between blur and ringing for learned filters, and the analysis of attention-based architectures, whose patch embeddings are strided operations without prefilter.

8.3 Continuous image representations

Several recent methods represent an image by a function of continuous coordinates. Chen, Liu and Wang (2021) learn a local implicit function that decodes features at any position, giving super-resolution at arbitrary, non-integer scales. Networks with Fourier features (Tancik and coauthors, 2020) or with periodic activations (Sitzmann and coauthors, 2020) fit images and scenes as coordinate-based functions; their ability to represent high frequencies is governed by the spectrum of the input encoding.

These models raise the questions of ch. 2 in a new form. What is the effective band of a learned representation? Does it interpolate the samples, and is it consistent with the acquisition model of ch. 1?

8.4 Multi-frame super-resolution

A single aliased image has lost information, but several images shifted by sub-pixel amounts sample the scene on a finer combined grid. Farsiu and coauthors (2004) formulate robust multi-frame super-resolution as a variational problem. Wronski and coauthors (2019) use the natural tremor of a hand holding a phone: the burst of frames, registered with sub-pixel accuracy, replaces demosaicking and increases resolution in commercial smartphones.

Here aliasing is the signal, and sub-pixel registration (ch. 6) is the bottleneck. Open questions concern the attainable resolution gain as a function of noise, of the number of frames and of the registration error, in the spirit of the Cramér–Rao analysis of ch. 6.2.

8.5 Texture synthesis beyond random phases

Random phase textures (ch. 7.2) are exact and fast but limited to micro-textures. Portilla and Simoncelli (2000) match a richer set of statistics of a wavelet decomposition. Gatys, Ecker and Bethge (2015) match the correlations of the features of a pretrained network, at a much higher cost. Galerne, Leclaire and Moisan (2017) design a texton, a small image whose spot noise reproduces a Gaussian texture, for fast on-the-fly synthesis in graphics. An open problem is to find models that keep the explicit, analysable character of RPN while capturing some phase structure.

8.6 Open problems

• Optimal prefilters for real sensors. Theorem 4.5 states a trade-off; which prefilter is optimal for a given PSF, noise level and downstream task is an open design question.
• Sub-pixel accuracy under aliasing. The Cramér–Rao bound of ch. 6.2 assumes a well-sampled model; quantifying the bias created by aliasing, such as pixel locking, remains largely empirical.
• Interpolation on non-periodic domains. The periodic plus smooth decomposition handles borders for the DFT. A theory of exact, lossless transforms for images with borders, masks or irregular sampling is incomplete.
• Reference-free quality. Phase-coherence indices measure sharpness; extending them to other degradations, such as compression, ringing or aliasing itself, is an active topic.

Exercises

1. ★ Poisson summation for a Gaussian. Using Theorem 0.2, show that the samples of the Gaussian density g_\sigma satisfy \sum_k g_\sigma(k) = \sum_m e^{-2\pi^2 \sigma^2 m^2} . Evaluate the deviation from 1 for \sigma = 1 and \sigma = 0.7.

   Solution: apply the first form of Theorem 0.2 at x = 0, with \hat g_\sigma(2\pi m) = e^{-\sigma^2 (2\pi m)^2 / 2}.

   The dominant correction is 2 e^{-2\pi^2\sigma^2}: about 5.4 \times 10^{-9} for \sigma = 1 and 1.26 \times 10^{-4} for \sigma = 0.7. The second value is the excess 1.000126 printed by Code 0.6.

2. ★ Folding frequencies. A signal sampled at the integers contains cosines of frequencies 0.7, 1.2 and -0.45 cycle per pixel. At which frequencies are they seen? In 2D, stripes have frequency (0.55, -0.3) cycle per pixel: give the apparent frequency, period and orientation.

   Solution: subtracting the nearest integer gives -0.3 (seen as 0.3 for a cosine), 0.2, and -0.45, which is already in the band. In 2D, the stripes are seen at (-0.45, -0.3), equivalent to (0.45, 0.3) for a cosine.

   The period changes from 1/0.627 = 1.60 to 1/0.541 = 1.85 pixels, and the direction of the frequency vector changes from -28.6° to 33.7°: the stripes turn by about 62°.

3. ★ A diffraction and pixel budget. A lens at N = 4 images light of wavelength 0.5 µm on a sensor of pitch 3 µm with a 100 % fill factor. Compute the optical cut-off, the Nyquist frequency, the optical and pixel attenuations at the Nyquist frequency, and the largest pitch that would avoid aliasing.

   Solution: \nu_c = 1/(0.5 \times 10^{-3} \times 4) = 500 cycles/mm, and the Nyquist frequency is 1/(2 \times 0.003) = 166.7 cycles/mm. With s = 1/3, Theorem 1.4 gives H = 0.584; the pixel gives 2/\pi = 0.637; the product is 0.372.

   Frequencies between 166.7 and 500 cycles/mm pass the optics and alias. Aliasing-free sampling needs p \le \lambda N / 2 = 1 µm.

4. ★ Cardinal cubic spline coefficients. Check that c_k = \sqrt 3\, z_1^{|k|}, with z_1 = \sqrt 3 - 2, solves b * c = \delta for b = (\tfrac16, \tfrac23, \tfrac16). Compute c_0, \dots, c_3 and the value of the cardinal spline \eta^3(1/2).

   Solution: at k = 0, \frac23 c_0 + \frac13 c_1 = \sqrt 3\, (2 + z_1)/3 = 1, since 2 + z_1 = \sqrt 3.

   At k \ge 1, the value is \sqrt 3\, z_1^{k-1} (1 + 4 z_1 + z_1^2) / 6 = 0, because z_1 is a root of z^2 + 4z + 1.

   Numerically, c = (1.7321, -0.4641, 0.1244, -0.0333). Then \eta^3(1/2) = \beta^3(1/2)(c_0 + c_1) + \beta^3(3/2)(c_1 + c_2) \approx 0.6005, to compare with \mathrm{sinc}(1/2) = 0.6366, Keys 0.5625 and linear 0.5.

5. ★ Shear amounts and padding. For a rotation by 30° with Theorem 3.10, compute the shear parameters. For an image of size 513 \times 513, find the largest displacement of the first shear and a padding that avoids any wrap-around.

   Solution: t = \tan 15° = 0.268 and s = \sin 30° = 0.5. Rows at distance 256 from the centre move by 0.268 \times 256 = 68.6 pixels in the first shear.

   Following the four corners through the three shears, every intermediate image stays within 256 (\cos 30° + \sin 30°) \approx 349.7 pixels of the centre in both coordinates. Zero-padding to 701 \times 701 (an odd size) therefore avoids wrap-around.

6. ★ How many heat iterations? One wants a Gaussian blur of \sigma = 3 pixels with the explicit heat scheme. How many iterations are needed with \delta t = 0.2? With \delta t = 1/4? Which time step damps all frequencies without sign change, and at what cost?

   Solution: t = \sigma^2 / 2 = 4.5, so 22.5 steps are needed; 23 steps give \sigma = 3.03. With \delta t = 1/4, 18 steps suffice, but the checkerboard is not damped (Theorem 5.6). The condition g \ge 0 requires \delta t \le 1/8, hence 36 steps.

7. ★★ Even sizes and the Nyquist coefficient. Prove item 2 of Proposition 2.6. Show that T_{1/2} T_{1/2} \neq T_1 for even M, and propose two remedies.

   Solution: the Nyquist term of U(k - \tau) is \frac1M \hat u_{M/2} \cos(\pi(k - \tau)) = \frac1M \hat u_{M/2} (-1)^k \cos(\pi\tau) .

   The integer translation T_1 multiplies the Nyquist coefficient by \cos\pi = -1, as a circular shift should, whereas T_{1/2} T_{1/2} multiplies it by \cos^2(\pi/2) = 0.

   Remedies: work with odd sizes, by cropping or padding one row and one column; or set the Nyquist coefficient to zero once, a negligible low-pass, after which all translations are unitary. This explains the error 0.30 of Code 2.7.

8. ★★ Periodic plus smooth in one dimension. For a signal u \in \mathbb{R}^M, write the one-dimensional version of Theorem 2.9. Show that the smooth component is the affine ramp s_k = b \Big( k - \frac{M-1}{2} \Big), \qquad b = \frac{u_{M-1} - u_0}{M} , and that the periodic jump of p is that of u divided by M.

   Solution: the boundary signal is v_0 = u_{M-1} - u_0, v_{M-1} = -v_0 and zero elsewhere. Inside, \Delta s = 0 forces s to be affine, s_k = a + bk.

   At k = 0, the periodic Laplacian equals s_{M-1} + s_1 - 2 s_0 = bM, so b = v_0 / M; the equation at k = M-1 then holds automatically, and the zero mean fixes a. Finally p_0 - p_{M-1} = (u_0 - u_{M-1}) + b(M-1) = (u_0 - u_{M-1}) / M.

9. ★★ Half-pixel gain of the Keys kernel. For the Keys kernel with parameter a, show that \varphi(1/2) = (4 - a)/8 and \varphi(3/2) = a/8. Deduce the half-pixel response, its value at \xi = \pi/2, and explain why it vanishes at \xi = \pi for every symmetric kernel.

   Solution: substituting into the two cubic pieces gives the two values, and pairing symmetric taps, |H_{1/2}(\xi)| = \Big| \frac{4 - a}{4} \cos\frac{\xi}{2} + \frac{a}{4} \cos\frac{3\xi}{2} \Big| . At \xi = \pi/2 this is \frac{\sqrt 2}{2}(1 - a/2): 0.707 for a = 0, 0.884 for a = -1/2, and 1.061 for a = -1, which amplifies and therefore sharpens.

   At \xi = \pi, the taps k and 1 - k have equal weights \varphi(k - 1/2) = \varphi(1/2 - k) and opposite signs (-1)^k, so they cancel.

10. ★★ No compact support on both sides. Let h \in L^1(\mathbb{R}) have compact support. Show that if \hat h vanishes on an open interval, then h = 0. Deduce that no prefilter with finite support has zero aliasing, \mathcal{A}(h) = 0.

    Solution: with compact support, \hat h(\zeta) = \int h(x) e^{-ix\zeta}\, dx is defined for every complex \zeta and is holomorphic, by differentiation under the integral. An entire function that vanishes on an interval vanishes everywhere (identity theorem), so \hat h = 0 and h = 0.

    Zero aliasing means \hat h = 0 on the open set |\xi| > \pi/z, which would force h = 0, contradicting \int h = 1.

11. ★★ The isotropic nine-point Laplacian. Prove Proposition 5.5. Expand the symbols of L_+ and L_\times to fourth order and find the unique combination that is isotropic at that order.

    Solution: for the five-point Laplacian, \hat L_+ = 2\cos\xi_1 + 2\cos\xi_2 - 4 = -|\xi|^2 + \frac{\xi_1^4 + \xi_2^4}{12} + O(|\xi|^6) .

    For the diagonal Laplacian, \hat L_\times = 2\cos\xi_1 \cos\xi_2 - 2, and expanding the product of cosines, \hat L_\times = -|\xi|^2 + \frac{\xi_1^4 + \xi_2^4}{12} + \frac{\xi_1^2 \xi_2^2}{2} + O(|\xi|^6) . The combination \alpha L_+ + (1 - \alpha) L_\times has the quartic part \frac{\xi_1^4 + \xi_2^4}{12} + \frac{(1-\alpha)\, \xi_1^2 \xi_2^2}{2} .

    It is proportional to |\xi|^4 exactly when (1 - \alpha)/2 = 1/6, that is \alpha = 2/3, which gives the stencil and the symbol -|\xi|^2 + |\xi|^4/12.

12. ★★ The aperture problem. Compute the structure matrix F of Theorem 6.3 for an image that depends only on \langle x, n \rangle for a unit vector n (a straight edge), then for a corner made of two straight edges with normals n_1 and n_2. What can be measured in each case?

    Solution: for a straight edge, \nabla U = U'(\langle x, n \rangle)\, n, so F = \big(\sum_k U'^2\big)\, n n^\top has rank one. The Fisher information along n^\perp is zero: only the displacement across the edge can be measured.

    For a corner, F = \alpha_1 n_1 n_1^\top + \alpha_2 n_2 n_2^\top is invertible when n_1 and n_2 are not parallel. Its smallest eigenvalue sets the worst-case precision, which is the quantity examined by the corner detector of Harris and Stephens (1988).

13. ★★ Bias of the parabolic peak fit (Python). For a peak f(x) = e^{-(x - \delta)^2 / 2} and for f(x) = \mathrm{sinc}(x - \delta), with \delta \in [-1/2, 1/2], compute the estimate of Proposition 6.7 from the samples at -1, 0, 1 and plot the bias. Repeat with the parabola fitted to \log f for the Gaussian.

    Hints and expected result: the bias is an odd, S-shaped function of \delta, zero at 0 and \pm 1/2. Its maximum is 0.048 pixel for the Gaussian and 0.119 pixel for the sinc; the logarithmic fit is exact for the Gaussian, up to round-off.

    import numpy as np
    
    def parabola_vertex(a, b, c):
        return (a - c) / (2 * (a - 2 * b + c))
    
    deltas = np.linspace(-0.5, 0.5, 101)
    for name, f in [("Gaussian peak, s = 1", lambda t: np.exp(-t**2 / 2)),
                    ("sinc peak", lambda t: np.sinc(t))]:
        est = np.array([parabola_vertex(f(-1 - d), f(-d), f(1 - d)) for d in deltas])
        print(f"{name}: max |bias| of the parabola fit = {np.max(np.abs(est - deltas)):.3f} pixel")
    g = lambda t: np.exp(-t**2 / 2)
    est_log = np.array([parabola_vertex(*np.log([g(-1 - d), g(-d), g(1 - d)])) for d in deltas])
    print(f"log-parabola on the Gaussian: max |bias| = {np.max(np.abs(est_log - deltas)):.1e} pixel")

14. ★★ The mean total variation of ADSN. Prove the formula for \mu in Proposition 7.9. Then take the M \times M image u = A\, \mathbf{1}[x < M/2], which has two vertical edges on the periodic domain, and compare \mathrm{TV}(u) with \mu.

    Solution: the formula follows from \mathbb{E}|X| = s \sqrt{2/\pi} for X \sim \mathcal{N}(0, s^2), summed over pixels and directions.

    For the edge image, \partial_y t = 0 and \partial_x t equals \pm A on two columns, so \|\partial_x t\|^2 = 2 M A^2 and \sigma_x = A\sqrt{2/M}.

    Hence \mu = 2 A M^{3/2} / \sqrt\pi, while \mathrm{TV}(u) = 2 M A.

    The ratio \mathrm{TV}(u) / \mu = \sqrt{\pi / M} equals 0.22 for M = 64: random phases spread the two edges over the whole image and multiply the total variation by more than four.

15. ★★★ Blind deblurring with the sharpness index (Python). Build a sharp image of random rectangles (96 \times 96), blur it with a Gaussian of \sigma_0 = 1.5 and add noise of standard deviation 1. For \sigma from 0 to 2.5, deblur with the regularised inverse filter \hat g_\sigma / (\hat g_\sigma^2 + 10^{-2}) and compute the sharpness index of the result (Code 7.10). Which \sigma maximises it? Discuss the shape of the curve.

    Hints and expected result: the script below reuses the formulas of Code 7.10. The index starts at about 12 for \sigma = 0, dips to about 8 near \sigma = 0.5, peaks at about 33 for \sigma = 1.5 — the true value — and falls to about 10 at \sigma = 2.5. Over-deblurring creates ringing, which lowers the phase coherence as surely as blur does. Extensions: motion-blur kernels, a finer grid, and the comparison with the approach of Leclaire and Moisan (2015).

    import numpy as np
    from scipy.special import log_ndtr
    
    def dx(u): return np.roll(u, -1, axis=1) - u
    def dy(u): return np.roll(u, -1, axis=0) - u
    def tv(u): return np.abs(dx(u)).sum() + np.abs(dy(u)).sum()
    
    def sharpness_index(u):
        MN = u.size
        g = [dx(u - u.mean()), dy(u - u.mean())]
        F = [np.fft.fft2(a) for a in g]
        sd = [np.sqrt(np.sum(a**2) / MN) for a in g]
        mu = MN * np.sqrt(2 / np.pi) * sum(sd)
        om = lambda t: t * np.arcsin(t) + np.sqrt(1 - t**2) - 1
        var = sum(MN * (2 / np.pi) * sd[i] * sd[j]
                  * om(np.clip(np.fft.ifft2(np.conj(F[i]) * F[j]).real / (MN * sd[i] * sd[j]), -1, 1)).sum()
                  for i in range(2) for j in range(2))
        return -log_ndtr(-(mu - tv(u)) / np.sqrt(var)) / np.log(10)
    
    rng = np.random.default_rng(10)
    n = 96
    u = np.zeros((n, n))
    for _ in range(30):                                  # random rectangles: a sharp scene
        x0, y0 = rng.integers(0, n - 10, 2); w, h = rng.integers(4, 25, 2)
        u[y0:y0 + h, x0:x0 + w] = rng.uniform(0, 255)
    q = 2 * np.pi * np.fft.fftfreq(n)
    G = lambda s: np.exp(-s**2 * (q[:, None]**2 + q[None, :]**2) / 2)
    sigma0 = 1.5
    v = np.fft.ifft2(np.fft.fft2(u) * G(sigma0)).real + rng.normal(0, 1.0, u.shape)
    
    def wiener(v, s, eps=1e-2):
        H = G(s)
        return np.fft.ifft2(np.fft.fft2(v) * H / (H**2 + eps)).real
    
    grid = np.arange(0.0, 2.6, 0.25)
    scores = [sharpness_index(wiener(v, s)) for s in grid]
    for s, sc in zip(grid, scores):
        print(f"sigma = {s:4.2f}: SI = {sc:8.1f}")
    print("argmax:", grid[int(np.argmax(scores))], " true blur:", sigma0)

16. ★★★ Aliasing inside a strided layer (Python). A one-dimensional “layer” applies a ReLU, an optional blur, and keeps one sample out of two. For band-limited random inputs, measure how far the output of the input shifted by one sample is from the exact half-sample shift of the original output. Compare no blur, the blur (1, 2, 1)/4 and the blur (1, 4, 6, 4, 1)/16, and relate the result to ch. 4 and ch. 8.2.

    Hints and expected result: with the code below, the relative errors are about 0.14, 0.05 and 0.02. The ReLU creates frequencies beyond the band, which the stride folds back; a prefilter reduces the folding (Theorem 4.5 says it cannot remove it without blur or ringing). This is the mechanism studied by Zhang (2019). Open-ended continuation: place the blur before the ReLU (the error is then about 0.13, barely better than no blur), and explain why.

    import numpy as np
    
    rng = np.random.default_rng(9)
    M = 254                                                # output length 127 (odd)
    
    def shift(u, tau):
        """Exact Shannon translation of a periodic signal (Nyquist coefficient dropped)."""
        p = np.fft.fftfreq(len(u), d=1 / len(u))
        H = np.exp(-2j * np.pi * p * tau / len(u)) * (np.abs(p) < len(u) / 2)
        return np.fft.ifft(np.fft.fft(u) * H).real
    
    def layer(x, blur):
        y = np.maximum(x, 0)                               # ReLU creates new high frequencies
        if blur is not None:
            k = np.array(blur, float) / np.sum(blur)
            y = sum(w * np.roll(y, j - len(k) // 2) for j, w in enumerate(k))
        return y[::2]                                      # stride 2
    
    p = np.fft.fftfreq(M, d=1 / M)
    for name, blur in [("stride only", None), ("blur [1,2,1]", [1, 2, 1]),
                       ("blur [1,4,6,4,1]", [1, 4, 6, 4, 1])]:
        e = []
        for _ in range(200):
            x = np.fft.ifft(np.fft.fft(rng.normal(size=M)) * (np.abs(p) < M / 8)).real
            y0 = layer(x, blur)
            y1 = layer(shift(x, 1.0), blur)                # input moved by 1 sample
            e.append(np.linalg.norm(y1 - shift(y0, 0.5)) / np.linalg.norm(y0 - y0.mean()))
        print(f"{name:18s} relative error to an exact half-sample shift: {np.mean(e):.3f}")

Projects and assessment

The official assessment is an exam (master-mva.com, “Mode de validation”); the Université Paris-Saclay course detail describes it as “Evaluation Continue non Intégrale”, a continuous assessment that is not the only component of the grade, and no weighting is displayed. The projects below are proposals for the continuous part or for personal training, in the MVA tradition of studying and reimplementing a research paper; the written exam proposal shows what a final exam on these notes could look like.

Solo projects

1. [Solo · ≈ 20 h] Paper study: the periodic plus smooth decomposition (Moisan, 2011). How much do border jumps degrade DFT-based processing, and how much does the decomposition remove?
   • Implement: Theorem 2.9 with two FFTs, its one-dimensional case (Exercise 8), and two alternatives — mirror symmetrisation and apodisation by a Hann window.
   • Data: synthetic non-periodic images and the sample photographs bundled with scikit-learn (sklearn.datasets.load_sample_image).
   • Measure: the energy on the frequency axes (Code 2.10), the ringing of a Fourier zoom by 4 near the borders, and the error of a translation by 0.5 pixel followed by its inverse.
   • Success criterion: the axis energy divided by at least 4, and a Laplacian of s below 10^{-10} inside the image.
   • Deliverables and evaluation: code, a 4–6 page report, a 10-minute talk, graded with the rubric below. Chapters used: ch. 2, 3.
2. [Solo · ≈ 20 h] Sub-pixel registration at the Cramér–Rao bound. How close to the bound of Theorem 6.3 do standard estimators get, and how does aliasing bias them?
   • Implement: in 2D, phase correlation with the two-sample formula (Foroosh, Zerubia and Berthod, 2002), upsampled cross-correlation (Guizar-Sicairos, Thurman and Fienup, 2008), and Gauss–Newton with exact Shannon translations.
   • Data: synthetic images rendered with pixel integration and a Gaussian PSF of width \sigma \in \{0.3, 0.5, 0.8, 1.0\}, plus noise; the scikit-learn sample images.
   • Success criterion: Gauss–Newton within 10 % of the bound for \sigma \ge 0.8; a plot of the bias against the true shift for \sigma \le 0.5 showing pixel locking.
   • Deliverables and evaluation: code, report, talk; rubric below. Chapters used: ch. 1, 2, 6.
3. [Solo · ≈ 15 h] Paper study: high-quality rotation by shears (Unser, Thévenaz and Yaroslavsky, 1995). Which rotation algorithm survives repeated application?
   • Implement: three-shear rotations with FFT shears and with cubic-spline shears, and direct rotation with bilinear and Keys kernels.
   • Experiment: the iterated-rotation test of Code 3.11 with n \in \{6, 12, 36, 72\} steps per full turn; RMSE, PSNR and spectra after the turn.
   • Success criterion: an FFT three-shear error below 10^{-6} on band-limited inputs, and a clear ranking of the other methods.
   • Deliverables and evaluation: code, report, talk; rubric below. Chapters used: ch. 2, 3.

Group projects

1. [Group of 2–3 · ≈ 30 h each] Paper study: random phase textures (Galerne, Gousseau and Morel, 2011). Which textures does a random phase reproduce, and why do the others fail?
   • Implement: RPN and ADSN for grey-level and colour images, applied to the periodic component of ch. 2.5. For colour, add the same random phase to the three channels, which keeps the phase differences between channels.
   • Data: textures cropped from the scikit-learn sample images and synthetic textures (filtered noise, random discs, regular patterns).
   • Evaluate: exact preservation of the power spectrum, histogram and kurtosis, the sharpness index of ch. 7.3, and a small blind visual comparison.
   • Split: member 1 implements grey RPN and ADSN with the decomposition; member 2 handles colour and the convergence of spot noise to ADSN as the number of spots grows; member 3 runs the evaluation and documents the failure cases.
   • Deliverables and evaluation: code, a joint 6-page report, a 15-minute talk; rubric below, with an individual question per member. Chapters used: ch. 2, 7.
2. [Group of 2–3 · ≈ 30 h each] Paper study: the sharpness index and blind deblurring (Blanchet and Moisan, 2012; Leclaire and Moisan, 2015). Can phase coherence select a deblurring kernel without any reference image?
   • Implement: SI with the closed-form moments of Proposition 7.9, checked by Monte Carlo; GPC by Monte Carlo for moderate values.
   • Experiments: SI as a function of blur, noise and over-deblurring; blind selection of the width of a Gaussian blur, then of the length of a horizontal motion blur, by maximising SI over regularised inverse filters (Exercise 15).
   • Success criterion: the selected Gaussian width within \pm 0.25 pixel of the truth on at least 80 % of a set of 20 test images.
   • Split: theory and verification of the moments; quality experiments; deblurring pipeline.
   • Deliverables and evaluation: code, joint report, talk; rubric below. Chapters used: ch. 4, 5, 7.
3. [Group of 2–3 · ≈ 30 h each] Paper study: anti-aliased downsampling in a small CNN (Zhang, 2019). Does a prefilter before each stride make a network more shift-consistent, and at what cost in accuracy?
   • Implement: a small PyTorch CNN, trained on CPU, with three downsampling variants: strided convolution, max-pooling, and blur-pooling with binomial filters of sizes 3 and 5.
   • Data: a synthetic data set of 32 \times 32 shapes (discs, squares, triangles) rendered at random sub-pixel positions with pixel integration (ch. 1.3), or the scikit-learn digits zoomed to 32 \times 32 by Fourier zero-padding.
   • Measure: accuracy, the probability that the predicted class is unchanged under one-pixel shifts, and the feature-level error of Exercise 16.
   • Split: data generation; model training; equivariance analysis in the light of Theorem 4.5.
   • Deliverables and evaluation: code, joint report, talk; rubric below. Chapters used: ch. 1, 2, 4, 8.

Written exam (proposal)

Two hours, one handwritten A4 sheet allowed, no electronic device; four exercises on 20 points.

1. Question (5 points). State the Shannon–Whittaker theorem with its hypotheses. Using the Poisson summation formula, explain what the samples of a non-band-limited image represent. Explain why diffraction makes optical images band-limited, why digital cameras nevertheless alias, and give the pitch condition for aliasing-free sampling at f-number N and wavelength \lambda. Expected: the band [-\pi, \pi] and the sinc series with its convergence; the periodisation formula of Proposition 1.7 and the folding of frequencies; the cut-off 1/(\lambda N), the condition p \le \lambda N / 2, and the remark that real pitches are larger, with attenuation by the optics and the pixel.

2. Question (5 points). For the explicit heat scheme with the nine-point Laplacian of Proposition 5.5, compute the symbol, its minimum on [-\pi, \pi]^2, the largest time step for \ell^2 stability, and the largest one for the maximum principle. How many steps give \sigma = 2 at the latter time step?

   Expected: the symbol is \frac16 \big( -20 + 8\cos\xi_1 + 8\cos\xi_2 + 4\cos\xi_1\cos\xi_2 \big), bilinear in the two cosines, with minimum -16/3 at (\pi, \pi).

   Stability holds for \delta t \le 3/8, and the maximum principle for \delta t \le 3/10, since the centre weight 1 - \frac{10}{3}\delta t must be non-negative.

   One step has variance 2\,\delta t per axis, as for the 5-point scheme, so \sigma = 2 needs t = 2.

   Seven steps of 0.3 give t = 2.1 and \sigma \approx 2.05; six steps of 0.3 followed by one of 0.2 give exactly t = 2.

3. Question (5 points). Prove the three-shear factorisation of a rotation. Show that a shear computed with Fourier translations is unitary on images of odd width and height, and deduce the inverse of the three-shear rotation. Explain why the half-pixel Fourier translation is not invertible for even sizes. Expected: the matrix computation with ts = 1 - \cos\theta (Theorem 3.10); unitarity by Parseval, since each row is multiplied by unimodular factors; the inverse is the rotation by -\theta; for even sizes, the Nyquist coefficient is multiplied by \cos(\pi/2) = 0 (Proposition 2.6).

4. Question (5 points). An edge U(x) = A\, \Phi(x / s) is observed with Gaussian noise of standard deviation \sigma after an unknown translation \tau. Derive the Cramér–Rao bound for \tau. Evaluate it for a contrast of 50, a blur s = 0.8 and a noise of standard deviation 1, first for one line and then for 64 lines crossing the edge. Explain why phase correlation alone does not reach the bound and how to reach it.

   Expected: the Fisher information is \frac{1}{\sigma^2} \sum_k U'(k - \tau)^2 \approx \frac{A^2}{2 s \sqrt\pi\, \sigma^2} , which gives a standard deviation of about 0.034 pixel for one line and 0.0042 for 64 lines.

   The normalisation of phase correlation weights noisy frequencies equally. Gauss–Newton least squares with exact translations, started from the phase-correlation estimate, reaches the bound.

Grading rubric (proposal)

References

Official references of the course (master-mva.com)

• L. Moisan. Modeling and Image Processing. Lecture notes of the course, available on the course web page (no date displayed).
• L. Moisan. Periodic plus smooth image decomposition. Journal of Mathematical Imaging and Vision, 39(2):161–179, 2011.

Textbooks

• J. W. Goodman. Introduction to Fourier Optics, 3rd edition. Roberts & Company, 2005.
• R. N. Bracewell. The Fourier Transform and Its Applications, 3rd edition. McGraw-Hill, 2000.
• S. Mallat. A Wavelet Tour of Signal Processing: The Sparse Way, 3rd edition. Academic Press, 2009.
• R. Hartley and A. Zisserman. Multiple View Geometry in Computer Vision, 2nd edition. Cambridge University Press, 2004.
• T. Lindeberg. Scale-Space Theory in Computer Vision. Kluwer Academic Publishers, 1994.
• V. Caselles and P. Monasse. Geometric Description of Images as Topographic Maps. Lecture Notes in Mathematics 1984, Springer, 2010.

Sampling, interpolation and geometric transforms

• C. E. Shannon. Communication in the presence of noise. Proceedings of the IRE, 37(1):10–21, 1949.
• M. Unser. Sampling — 50 years after Shannon. Proceedings of the IEEE, 88(4):569–587, 2000.
• R. G. Keys. Cubic convolution interpolation for digital image processing. IEEE Transactions on Acoustics, Speech, and Signal Processing, 29(6):1153–1160, 1981.
• A. Aldroubi, M. Unser and M. Eden. Cardinal spline filters: stability and convergence to the ideal sinc interpolator. Signal Processing, 28(2):127–138, 1992.
• M. Unser, A. Aldroubi and M. Eden. B-spline signal processing: Part I — Theory; Part II — Efficient design and applications. IEEE Transactions on Signal Processing, 41(2):821–848, 1993.
• M. Unser, A. Aldroubi and M. Eden. Enlargement or reduction of digital images with minimum loss of information. IEEE Transactions on Image Processing, 4(3):247–258, 1995.
• M. Unser. Splines: a perfect fit for signal and image processing. IEEE Signal Processing Magazine, 16(6):22–38, 1999.
• P. Thévenaz, T. Blu and M. Unser. Interpolation revisited. IEEE Transactions on Medical Imaging, 19(7):739–758, 2000.
• A. W. Paeth. A fast algorithm for general raster rotation. Proceedings of Graphics Interface ’86, 77–81, 1986.
• M. Unser, P. Thévenaz and L. Yaroslavsky. Convolution-based interpolation for fast, high-quality rotation of images. IEEE Transactions on Image Processing, 4(10):1371–1381, 1995.

Differential operators and iterative filtering

• J. J. Koenderink. The structure of images. Biological Cybernetics, 50(5):363–370, 1984.
• P. Perona and J. Malik. Scale-space and edge detection using anisotropic diffusion. IEEE Transactions on Pattern Analysis and Machine Intelligence, 12(7):629–639, 1990.
• L. Alvarez, F. Guichard, P.-L. Lions and J.-M. Morel. Axioms and fundamental equations of image processing. Archive for Rational Mechanics and Analysis, 123(3):199–257, 1993.
• L. Moisan. Affine plane curve evolution: a fully consistent scheme. IEEE Transactions on Image Processing, 7(3):411–420, 1998.

Sub-pixel geometry

• B. D. Lucas and T. Kanade. An iterative image registration technique with an application to stereo vision. Proceedings of the International Joint Conference on Artificial Intelligence (IJCAI), 674–679, 1981.
• J. Canny. A computational approach to edge detection. IEEE Transactions on Pattern Analysis and Machine Intelligence, 8(6):679–698, 1986.
• C. Harris and M. Stephens. A combined corner and edge detector. Proceedings of the 4th Alvey Vision Conference, 147–151, 1988.
• F. Devernay. A non-maxima suppression method for edge detection with sub-pixel accuracy. INRIA Research Report RR-2724, 1995.
• P. Monasse and F. Guichard. Fast computation of a contrast-invariant image representation. IEEE Transactions on Image Processing, 9(5):860–872, 2000.
• H. Foroosh, J. B. Zerubia and M. Berthod. Extension of phase correlation to subpixel registration. IEEE Transactions on Image Processing, 11(3):188–200, 2002.
• D. Robinson and P. Milanfar. Fundamental performance limits in image registration. IEEE Transactions on Image Processing, 13(9):1185–1199, 2004.
• J. Delon and B. Rougé. Small baseline stereovision. Journal of Mathematical Imaging and Vision, 28(3):209–223, 2007.
• M. Guizar-Sicairos, S. T. Thurman and J. R. Fienup. Efficient subpixel image registration algorithms. Optics Letters, 33(2):156–158, 2008.

Phase, textures and sharpness

• A. V. Oppenheim and J. S. Lim. The importance of phase in signals. Proceedings of the IEEE, 69(5):529–541, 1981.
• J. J. van Wijk. Spot noise texture synthesis for data visualization. Computer Graphics (Proceedings of SIGGRAPH ’91), 25(4):309–318, 1991.
• P. Kovesi. Image features from phase congruency. Videre: Journal of Computer Vision Research, 1(3), 1999.
• Z. Wang and E. P. Simoncelli. Local phase coherence and the perception of blur. Advances in Neural Information Processing Systems 16 (NIPS 2003), 2004.
• G. Blanchet, L. Moisan and B. Rougé. Measuring the global phase coherence of an image. Proceedings of the IEEE International Conference on Image Processing (ICIP), 2008.
• B. Galerne, Y. Gousseau and J.-M. Morel. Random phase textures: theory and synthesis. IEEE Transactions on Image Processing, 20(1):257–267, 2011.
• G. Blanchet and L. Moisan. An explicit sharpness index related to global phase coherence. Proceedings of the IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2012.
• A. Leclaire and L. Moisan. No-reference image quality assessment and blind deblurring with sharpness metrics exploiting Fourier phase information. Journal of Mathematical Imaging and Vision, 52(1):145–172, 2015.

Research directions (ch. 8)

• J. Portilla and E. P. Simoncelli. A parametric texture model based on joint statistics of complex wavelet coefficients. International Journal of Computer Vision, 40(1):49–70, 2000.
• S. Farsiu, M. D. Robinson, M. Elad and P. Milanfar. Fast and robust multiframe super resolution. IEEE Transactions on Image Processing, 13(10):1327–1344, 2004.
• L. A. Gatys, A. S. Ecker and M. Bethge. Texture synthesis using convolutional neural networks. Advances in Neural Information Processing Systems (NIPS), 2015.
• R. Abergel and L. Moisan. The Shannon total variation. Journal of Mathematical Imaging and Vision, 59(2):341–370, 2017.
• L. Condat. Discrete total variation: new definition and minimization. SIAM Journal on Imaging Sciences, 10(3):1258–1290, 2017.
• B. Galerne, A. Leclaire and L. Moisan. Texton noise. Computer Graphics Forum, 36(8), 2017.
• A. Azulay and Y. Weiss. Why do deep convolutional networks generalize so poorly to small image transformations? Journal of Machine Learning Research, 20, 2019.
• R. Zhang. Making convolutional networks shift-invariant again. Proceedings of the International Conference on Machine Learning (ICML), 2019.
• B. Wronski, I. Garcia-Dorado, M. Ernst, D. Kelly, M. Krainin, C.-K. Liang, M. Levoy and P. Milanfar. Handheld multi-frame super-resolution. ACM Transactions on Graphics (Proceedings of SIGGRAPH), 38(4), 2019.
• M. Tancik, P. P. Srinivasan, B. Mildenhall, S. Fridovich-Keil, N. Raghavan, U. Singhal, R. Ramamoorthi, J. T. Barron and R. Ng. Fourier features let networks learn high frequency functions in low dimensional domains. Advances in Neural Information Processing Systems (NeurIPS), 2020.
• V. Sitzmann, J. N. P. Martel, A. W. Bergman, D. B. Lindell and G. Wetzstein. Implicit neural representations with periodic activation functions. Advances in Neural Information Processing Systems (NeurIPS), 2020.
• Y. Chen, S. Liu and X. Wang. Learning continuous image representation with local implicit image function. Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), 2021.
• T. Karras, M. Aittala, S. Laine, E. Härkönen, J. Hellsten, J. Lehtinen and T. Aila. Alias-free generative adversarial networks. Advances in Neural Information Processing Systems (NeurIPS), 2021.