Detection Theory and Industrial Applications

Contents

• Introduction and overview
• 0. Review: tails, union bounds and p-values
• 1. Perception, Gestalt theory and the discrete image
• 2. Statistical principles of detection
• 3. The a-contrario framework: background model and NFA
• 4. Detection of geometric structures
• 5. Anomaly detection in images and hyperspectral data
• 6. Detection in time series and recurrent satellite imaging
• 7. Detection of image forgeries
• 8. Research directions and open problems
• Exercises
• Projects and assessment
• References

Introduction and overview

Audience: a student with a good M1 in mathematics, computer science or physics who knows elementary probability and wants detection algorithms whose thresholds are justified rather than tuned.

Empirical first: threshold a white-noise image of 512 \times 512 pixels at three standard deviations above its mean and about 350 pixels light up, although there is nothing to see. The eye, on the same image, sees nothing at all. A detector must behave like the eye: it must count how many times the observed event could happen by chance among all the events it looked at, and stay silent on noise.

Detection is the most frequent task of image, video and time-series analysis in industry: find the defect on a wafer, the change in a satellite time series, the pasted region in a photograph. It differs from classification because the objects to detect are rare and seldom come with enough examples to train a classifier.

The course builds one answer, the a-contrario theory. It models the absence of structure (the background), not the structure, and declares a detection when the expected number of such observations under the background model, the number of false alarms (NFA), is small.

The tools come from Probability and Statistics and from the image-processing courses of semester 1, notably Sub-Pixel Image Processing and Mathematical Models for Digital Photography. The applications meet Remote Sensing Data and Machine Learning for Time Series.

What you will be able to do

• Explain the Gestalt grouping laws and the Helmholtz principle, and turn a perceptual grouping into a statistical test.
• Bound binomial and Gaussian tails, and choose between exact tails, Chernoff and Hoeffding bounds.
• Prove the Neyman–Pearson lemma and derive the matched filter for a known shape in Gaussian noise.
• Control the mean number of false discoveries of a family of tests with the Bonferroni rule, and compare it with false discovery rate control.
• Define a background model and a number of false alarms, and prove that the expected number of false detections is at most \varepsilon.
• Count the tests of a detection problem and compute the resulting detection thresholds.
• Detect alignments and line segments with a false-detection guarantee, as the LSD algorithm does.
• Detect anomalies in multi-band images with RX, mixture and matched-filter detectors, with thresholds set by the NFA.
• Detect, and classify by type, changes in satellite image time series after a cloud screening.
• Detect forged regions of an image from inconsistencies of its JPEG grid.

Official objectives and syllabus (2026–2027)

The course description on master-mva.com, written by the instructors’ team at the Centre Borelli, starts from industrial problems whose goal is the detection of objects or events in images or image series. It stresses that detection is not classification, because detection usually lacks a sufficient number of examples to be reduced to a classification problem.

It introduces the a-contrario theory, a statistical framework that selects detection thresholds automatically, so as to control the expected number of false detections under given background models. It covers the theory and practice of image and video processing needed for automatic detection, through several examples that include real-world applications.

The Université Paris-Saclay course detail sets the learning objective of training students in detection, which it calls one of the most frequent applications of image, video and time-series analysis. Its programme has eight items; master-mva.com lists five themes. The table matches every official item with the section of these notes that teaches it.

Beyond the syllabus, these notes add the exact Student test for changes in a seasonal series (ch. 6.2), the Hotelling correction of local RX (ch. 5.2), and a self-study chapter on learned background models, conformal p-values and the forensics of generated images (ch. 8).

Where this course fits in the programme

The official pages display no prerequisites. The course uses probability at the level of a good M1: discrete and Gaussian distributions, conditional independence, hypothesis tests, and a little linear algebra.

• Preparation — the refresher Probability and Statistics covers the tests, Gaussian vectors and the chi-square distribution of ch. 0 and ch. 2; Linear Algebra, Analysis and Optimization covers orthogonal projections and the Mahalanobis distance; Scientific Python, NumPy and PyTorch covers the vectorised code of every chapter.
• Semester-1 courses that help — Sub-Pixel Image Processing (Shannon sampling, used in ch. 1 and 4), Mathematical Models for Digital Photography (camera pipeline, demosaicking and noise, used in ch. 7), Image Denoising (noise models and residuals), Machine Learning for Time Series (ch. 6) and Statistical Learning with Extreme Values (tails of rare events).
• Related courses of the same semester — Remote Sensing Data (satellite sensors and time series), Inverse Problems and Imaging (statistical image models), Generative Models for Images (texture models, and the generated images that forensics must detect) and Interactions: Complexity, Stochasticity, Extrema and Rare Events.
• Afterwards — the Research Internship: detection problems with few examples are frequent internship subjects in industry and in imaging laboratories.

Course outline

• Review: tails, union bounds and p-values — binomial and Gaussian tails, Chernoff and Hoeffding bounds, the union bound without independence, super-uniform p-values.
• Perception, Gestalt theory and the discrete image — what vision groups, the Gestalt laws, the Helmholtz principle, Shannon sampling and Attneave’s view of visual information.
• Statistical principles of detection — detection versus classification, the Neyman–Pearson lemma, the matched filter, multiple testing and the mean number of false discoveries.
• The a-contrario framework — background models, the number of false alarms, the fundamental theorem, thresholds, Gaussian models and maximality.
• Detection of geometric structures — meaningful alignments, the LSD line segment detector and other geometric structures.
• Anomaly detection in images and hyperspectral data — families of anomaly detectors, the RX detector and its variants, the a-contrario reduction of anomalies to detection in noise.
• Detection in time series and recurrent satellite imaging — exact change tests in seasonal series, cloud screening, detection and categorisation of changes.
• Detection of image forgeries — JPEG grid, demosaicking, noise and copy-move traces, and the control of false alarms on authentic images.
• Research directions and open problems — self-study: dependence, learned background models, conformal p-values, generated images.

Session plan

24 h in class = 8 sessions of 3 h. The volume is the official one of the Université Paris-Saclay programme table (lectures 24 h). master-mva.com announces weekly online lectures of 3 hours combining theory and the discussion of particular cases; the Université Paris-Saclay detail announces 7 to 9 three-hour sessions, each split into theory and a case analysis, and lists distance teaching as “non”.

The split below is a proposal. Each session gives about two hours of theory and one hour of case analysis with the code of the chapter. Level 1 = foundations, 2 = core results and algorithms, 3 = advanced and research-level material; the level never decreases.

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

How to study this module

• Read ch. 0 before the first session: every NFA of the course is a tail probability multiplied by a number of tests, and the bounds of ch. 0 tell you how such numbers behave.
• After each session, run the code of the chapter and change one parameter at a time (number of tests, precision, contrast). The numbers quoted in the text were produced by these scripts; they need Python 3 with NumPy, SciPy and scikit-learn.
• For each detector of the course, write down three things before reading the algorithm: the background model, the statistic, and the number of tests.
• Do the exercises listed for each session before the next one; the ★★ proofs are typical homework questions.
• For the paper reports, pick an article of the IPOL journal (ipol.im) with an online demo, as the Université Paris-Saclay detail suggests, and test it on images of your own.

0. Review: tails, union bounds and p-values

Every detection threshold of this course is a tail probability multiplied by a number of tests. This review gathers the M1-level tools that compute such tails and that count false alarms without any independence assumption between the tests.

0.1 Binomial tails

Many statistics of the course count how many of n independent observations are “aligned”, “voting” or “changed”, each with probability p under the background model.

The exact tail is cheap to compute (SciPy: binom.sf(k - 1, n, p)), but bounds explain how thresholds grow with the number of tests.

0.2 Gaussian tails and quadratic forms

Let \varphi be the standard normal density and \Phi^c(z) = \mathbb{P}(Z \ge z) the Gaussian tail.

For instance, \Phi^c(4) = 3.17 \times 10^{-5} lies between the two bounds 3.15 \times 10^{-5} and 3.35 \times 10^{-5}.

In several bands, the Gaussian tail becomes the tail of a quadratic form. If X \sim \mathcal{N}(\mu, \Sigma) in \mathbb{R}^d with \Sigma invertible, then Y = \Sigma^{-1/2}(X - \mu) is standard normal. Hence the squared Mahalanobis distance (X - \mu)^\top \Sigma^{-1} (X - \mu) = \|Y\|^2 follows the chi-square distribution with d degrees of freedom. This is the basis of the RX detector of ch. 5.

0.3 Counting events without independence

The central computation of the course is the expected number of events that occur among many.

The second inequality is the union bound. Neither statement needs the events to be independent: this is why the a-contrario theory can test millions of overlapping candidates.

0.4 Super-uniform p-values

The proposition matters for discrete statistics such as binomial counts. The test “reject when the p-value is at most \alpha” never has a level larger than \alpha, even though the p-value takes only finitely many values.

Code 0.8 — Binomial tails and their bounds. The script reproduces the table of Example 0.3.

import numpy as np
from scipy.stats import binom

def tail(n, k, p):
    """B(n, k, p) = P(S >= k) for S ~ Binomial(n, p)."""
    return binom.sf(k - 1, n, p)

def kl(a, p):
    """Kullback-Leibler divergence between Bernoulli(a) and Bernoulli(p)."""
    return a * np.log(a / p) + (1 - a) * np.log((1 - a) / (1 - p))

n, p = 100, 1 / 8
for k in (20, 25, 30, 40):
    a = k / n
    print(f"k={k}: exact {tail(n, k, p):.2e}   "
          f"Chernoff {np.exp(-n * kl(a, p)):.2e}   "
          f"Hoeffding {np.exp(-2 * n * (a - p) ** 2):.2e}")

1. Perception, Gestalt theory and the discrete image

This chapter asks what a detector should detect. Human vision gives the qualitative answer, the Gestalt laws; the Helmholtz principle makes it quantitative; the discrete nature of images tells which measurements can be treated as independent.

1.1 What vision does

An image is an array of a few million numbers. Seeing consists in interpreting this array in terms of objects, their boundaries and their arrangement. The first step of this interpretation is a grouping: local elements (dots, edge pieces, regions of similar colour) are gathered into wholes that are perceived at once.

Human vision has striking properties that a detection theory must explain.

• It sees alignments, symmetries and closed curves immediately, with no conscious effort.
• It sees nothing in white noise, although noise contains every local pattern somewhere.
• It can be forced to see what is not there: in Kanizsa’s triangle, three black discs with a missing sector create a white triangle whose edges are drawn nowhere.
• It can fail to see what is there, when a shape is embedded in a stronger structure (masking).

1.2 The Gestalt laws

Gestalt psychology, from Wertheimer’s 1923 article to the work of Kanizsa (1979), described the grouping laws of perception. Each law, called a partial gestalt, groups elements that share one quality.

• Proximity: elements close to each other are grouped.
• Similarity: elements of similar colour, size or orientation are grouped.
• Continuity of direction (good continuation): elements that continue the same direction form a curve; aligned elements form a line.
• Parallelism, symmetry, constant width: pairs of curves with these relations are grouped.
• Closure and convexity: closed and convex curves are perceived as the boundaries of objects.
• Amodal completion: at a T-junction, the occluded edge is perceived as continuing behind the occluder.
• Common fate: elements moving together are grouped.

Grouping is recursive. A group, for instance a line formed by aligned dots, becomes an element of a higher-level group, for instance a pair of parallel lines.

Several laws can act together and reinforce each other, or conflict. Kanizsa showed that a figure disappears when it is masked, for example when it is embedded in a texture or when added lines create a stronger gestalt that absorbs it.

The Gestalt laws are qualitative. They do not say how close is “close” or how aligned is “aligned”. A computational theory must provide these thresholds, and it must do so without tuning them for each image.

1.3 The Helmholtz principle

The quantitative answer of Desolneux, Moisan and Morel (2007) is a principle that they attribute to Helmholtz: we perceive at once any structure that is too regular to be there by chance. The same idea appears in computer vision as the non-accidentalness of configurations (Lowe, 1985).

The example contains the whole method of the course: a background model (uniform points), a list of tests (the strips), a statistic (the count) and a threshold derived from a bound on the expected number of false detections.

1.4 The discrete nature of images

A digital image is a set of samples u(k), k \in \mathbb{Z}^2, of a continuous image blurred by the optics. If the continuous image is band-limited, its spectrum being supported in [-\pi, \pi]^2, the Shannon–Whittaker theorem says that it is determined by its samples: u(x) = \sum_{k \in \mathbb{Z}^2} u(k)\, \mathrm{sinc}(x_1 - k_1)\, \mathrm{sinc}(x_2 - k_2) . Here \mathrm{sinc}(t) = \sin(\pi t) / (\pi t). The theorem and its proof are in Sub-Pixel Image Processing.

Three consequences matter for detection.

• An image of N \times M pixels has NM degrees of freedom: the number of pixels is the natural unit when tests are counted.
• Local measurements, such as a gradient computed on a 2 \times 2 window, use several samples. Measurements on overlapping windows are dependent, so a statistic that assumes independence must use measurements on disjoint supports.
• A measured orientation has a limited precision, and a detector must use a tolerance rather than exact equality.

The gradient used by the LSD detector of ch. 4 makes the second point exact. With a = u(i, j), b = u(i, j+1), c = u(i+1, j) and d = u(i+1, j+1), it is g_x = \frac{b + d - a - c}{2}, \qquad g_y = \frac{c + d - a - b}{2} .

This is why the alignment detector of ch. 4.1 samples orientations at points two pixels apart. The script below checks the three statements: the empirical covariance matrix is close to the identity, the orientation histogram passes a chi-square test of uniformity, and the correlation of \cos\theta is -0.41 at distance 1 and 0.00 at distance 2.

Code 1.5 — Orientations of a noise image.

import numpy as np
from scipy.stats import chisquare

rng = np.random.default_rng(0)
u = rng.normal(size=(512, 512))                      # white Gaussian noise image

# 2x2 gradient (the scheme of LSD): differences averaged over a 2x2 window
gx = (u[:-1, 1:] + u[1:, 1:] - u[:-1, :-1] - u[1:, :-1]) / 2
gy = (u[1:, :-1] + u[1:, 1:] - u[:-1, :-1] - u[:-1, 1:]) / 2
theta = np.arctan2(gy, gx)

print("var gx, var gy, cov:", np.round([gx.var(), gy.var(), np.mean(gx * gy)], 3))
counts, _ = np.histogram(theta, bins=16, range=(-np.pi, np.pi))
print("chi-square p-value of uniformity: %.2f" % chisquare(counts).pvalue)

# dependence: correlation of cos(theta) between neighbours at distance 1 and 2
c = np.cos(theta)
for d in (1, 2):
    r = np.corrcoef(c[:, :-d].ravel(), c[:, d:].ravel())[0, 1]
    print(f"corr(cos theta) at horizontal distance {d}: {r:+.3f}")

1.5 Attneave’s informational view of perception

Attneave (1954) applied Shannon’s information theory to perception. His observation is that images are highly redundant. A uniform region is predictable from a few of its pixels; information concentrates where the image changes, on the contours. Along a contour, it concentrates where the direction changes most, at the points of maximal curvature.

He illustrated the claim with a drawing of a sleeping cat made of 38 points of maximal curvature joined by straight segments, which remains immediately recognisable.

Two lessons follow for detection.

• Detectors should work on geometric features such as level lines, gradient orientations and corners, rather than on raw grey levels. Orientations have an extra advantage: they do not change under a monotone change of contrast.
• The surprise of an event, in Shannon’s sense, is -\log \mathbb{P}(\text{event}). The a-contrario theory measures the surprise of a group by -\log_{10} of its number of false alarms (ch. 3), which is this self-information minus the logarithm of the number of groups examined.

2. Statistical principles of detection

Detection is a statistical decision. This chapter recalls the classical theory of tests, shows what it gives when the object to detect can be modelled, and explains why testing many candidates at once changes the question from “is this test significant?” to “how many false detections do I expect?”.

2.1 Detection is not classification

A classifier learns from labelled examples of every class and minimises a misclassification risk. Detection problems in industry rarely fit this frame.

• The background is abundant, but the objects to detect (defects, changes, forgeries) are rare and variable; there are seldom enough examples to learn them.
• The two errors have very different costs. A miss loses one object; false alarms, repeated over millions of pixels and thousands of images, make the detector useless to the operator who must check them.
• The decision is taken for many candidates at once: every position, scale and orientation of the image is a candidate.

The theory of this course therefore models the background only, and controls the false alarms first.

2.2 Tests and the Neyman–Pearson lemma

2.3 Detecting a known shape: the matched filter

When the object can be modelled, the lemma gives the optimal detector.

Two extensions are used later.

• Coloured noise. If w \sim \mathcal{N}(0, C), whitening by C^{-1/2} reduces to the white case. The statistic becomes T(x) = \frac{s^\top C^{-1} x}{\sqrt{s^\top C^{-1} s}} \sim \mathcal{N}(0, 1) \text{ under } H_0 .
• Unknown position or scale. The shape s_i is known up to an index i (position, scale, orientation). The generalised likelihood ratio test takes the maximum of T_i over i. This is a family of tests, one per candidate, and its false alarms must be counted globally (§ 2.4).

When nothing is known about the object, only H_0 remains. Fisher’s significance test then rejects H_0 when the p-value of a statistic is small. Proposition 0.7 guarantees its level, whatever the alternative.

2.4 Multiple testing

An image of 512 \times 512 pixels tested at four scales gives about a million tests. At level 10^{-3} each, a thousand false alarms are expected on pure noise. This is the phenomenon of the empirical-first example.

The FDR adapts to the number of discoveries: it rejects more when many signals are present. The Bonferroni rule is often called conservative, but this depends on which error rate it is asked to control.

2.5 The mean number of false discoveries (Gordon et al., 2007)

Gordon, Glazko, Qiu and Yakovlev (2007) argued that the reputation of Bonferroni comes from asking it to control the FWER at a small level. They advocate reading the Bonferroni rule as a control of the PFER, \mathbb{E}[V] \le \varepsilon. A level that is not small, such as \varepsilon = 1, then makes sense.

They also reported the first study of the stability of the two procedures. The Bonferroni rule showed a smaller variance of the number of true discoveries and of the total number of discoveries than Benjamini–Hochberg, which matters most when the p-values are correlated.

Code 2.8 reproduces the phenomenon on 2000 Gaussian tests, 100 of them with a signal of 3.5 standard deviations, with the Bonferroni rule at \varepsilon = 1 and Benjamini–Hochberg at q = 0.05. The correlation is introduced by a common factor shared by all the tests.

The Bonferroni rule keeps \mathbb{E}[V] = N_0 / N = 0.95 whatever the correlation, as Proposition 2.5 predicts. Its number of discoveries is three times more stable than that of Benjamini–Hochberg under correlation, at the price of fewer discoveries.

The mean number of false discoveries is exactly the quantity that the a-contrario theory controls in images, where the tests are always strongly correlated.

Code 2.8 — Bonferroni and Benjamini–Hochberg under correlation.

import numpy as np
from scipy.stats import norm

rng = np.random.default_rng(1)
m, m1, shift, runs = 2000, 100, 3.5, 4000
eps, q = 1.0, 0.05                                   # PFER target, FDR target

def bh(pv, q):
    """Benjamini-Hochberg step-up: boolean mask of rejections."""
    order = np.argsort(pv)
    ok = pv[order] <= q * np.arange(1, len(pv) + 1) / len(pv)
    k = ok.nonzero()[0].max() + 1 if ok.any() else 0
    rej = np.zeros(len(pv), bool)
    rej[order[:k]] = True
    return rej

truth = np.zeros(m, bool); truth[:m1] = True
for rho in (0.0, 0.5):
    V = {"Bonferroni": [], "BH": []}; R = {"Bonferroni": [], "BH": []}
    for _ in range(runs):
        common = rng.normal()                         # one-factor equicorrelation
        z = np.sqrt(rho) * common + np.sqrt(1 - rho) * rng.normal(size=m) + shift * truth
        pv = norm.sf(z)                               # one-sided p-values
        for name, rej in (("Bonferroni", pv <= eps / m), ("BH", bh(pv, q))):
            V[name].append(np.sum(rej & ~truth)); R[name].append(rej.sum())
    for name in V:
        print(f"rho={rho}  {name:10s}  E[V]={np.mean(V[name]):5.2f}  sd(V)={np.std(V[name]):5.2f}"
              f"  E[R]={np.mean(R[name]):6.1f}  sd(R)={np.std(R[name]):5.1f}")

3. The a-contrario framework: background model and NFA

This chapter turns the Helmholtz principle of ch. 1 and the Bonferroni rule of ch. 2 into one definition, the number of false alarms. It proves the fundamental property of the definition and computes the thresholds it produces.

3.1 The background model

A good background model follows three principles.

• It is as weak as possible, typically independence and invariance, so that rejecting it says something about the data.
• It makes the probability of the observed statistic computable.
• It does not contain the structure sought.

Typical choices are independent uniform orientations (justified for noise by Proposition 1.4), independent uniform points, Gaussian white or coloured noise, and exchangeable residuals.

3.2 The number of false alarms

Let the candidates be indexed by i \in \{1, \dots, N_T\}, where N_T is the number of tests. To each candidate is attached a statistic X_i, larger when the candidate is more structured.

More generally, any function of the candidate and the data is called an NFA if it satisfies the property of the next theorem.

The NFA is the Bonferroni rule of Proposition 2.5 written as a single number per candidate. The threshold \varepsilon is the mean number of false detections that the user accepts per image; the default is \varepsilon = 1.

Different families of candidates often have different sizes, for instance segments of different lengths. The tests can then be weighted.

3.3 What the NFA is and is not

• The NFA is not a probability, and it is not the probability that the object exists. It is an expected count under the background model; it can be larger than 1.
• The guarantee holds under H_0 only. When the detector fires on a real image, either a structure is present or the background model is wrong there. The detection is a local rejection of H_0, hence the name a contrario.
• Independence is needed inside a test, to compute the tail of X_i. It is not needed between tests.
• -\log_{10} \mathrm{NFA} is a measure of meaningfulness: a segment with \mathrm{NFA} = 10^{-12} is twelve orders of magnitude beyond chance.

3.4 Thresholds and their stability

Detection thresholds depend on \varepsilon and on N_T only through \log(N_T / \varepsilon).

For N_T = 65\,536 tests, the Gaussian thresholds are:

Multiplying \varepsilon by 10^4 moves the threshold by about two standard deviations. This weak dependence is why \varepsilon = 1 is used everywhere, and why a-contrario detectors are often said to have no detection parameter: the user fixes the accepted number of false alarms, not a threshold on the statistic.

3.5 Gaussian models of the background and of the shape

The Université Paris-Saclay programme asks for Gaussian modelling of the background and of the shape to be detected. Both appear naturally.

• Background. The image u is a Gaussian white noise of mean \mu and variance \sigma^2, estimated robustly on the image; or a coloured Gaussian noise of covariance C, a model of textured backgrounds used by Grosjean and Moisan (2009).
• Shape. The object is a known profile s_i (a disc of radius r at position x, a bar, a blob) added with an amplitude a: u = \mu + a s_i + w.

The background alone defines the NFA. For a template s_i in white noise, T_i = \frac{\langle s_i, u - \mu \rangle}{\sigma \|s_i\|}, \qquad \mathrm{NFA}(i) = N_T\, \Phi^c(T_i) . For coloured noise, T_i is the whitened statistic of ch. 2.3.

The shape model gives the detectability: under the shape model, T_i \sim \mathcal{N}(a\|s_i\|/\sigma, 1). The probability of detection at level \varepsilon is \mathbb{P}(\text{detect}) = \Phi^c\big(z(\varepsilon) - a \|s_i\| / \sigma\big), and an object is detected with probability one half when its amplitude reaches a_{\min} = \sigma z(\varepsilon) / \|s_i\|.

Code 3.9 checks these numbers on simulated images; its output is discussed in § 3.6.

3.6 Redundancy and maximality

Overlapping candidates see the same object several times. A bright spot is detected at several neighbouring centres and radii; a long alignment contains many meaningful sub-segments.

Reporting maximal candidates does not increase the number of false detections, since it only removes candidates from the list.

In Code 3.9, the mean number of 1-meaningful tests on pure noise is 1.05, close to \varepsilon, as Theorem 3.3 predicts. After grouping connected detections, the mean number of detected spots is 0.80, because some false detections come in clusters.

With a spot of contrast 0.5, below a_{\min}, the true spot is missed (\log_{10} \mathrm{NFA} = 0.77) and one false spot appears elsewhere. With contrast 1, the spot is found with \log_{10} \mathrm{NFA} = -7.2 and no false spot appears.

Code 3.9 — A-contrario detection of spots in Gaussian noise.

import numpy as np
from scipy.stats import norm
from scipy.ndimage import label

rng = np.random.default_rng(5)
n, radii = 128, (1, 2, 3, 4)
yy, xx = np.mgrid[:n, :n]
dist2 = np.minimum(yy, n - yy) ** 2 + np.minimum(xx, n - xx) ** 2   # periodic distance to 0
discs = [(dist2 <= r * r).astype(float) for r in radii]
fft_discs = [np.fft.fft2(d) for d in discs]
N_tests = n * n * len(radii)

def log10_nfa(u):
    """log10 NFA of a bright disc centred at each pixel, for each radius (sigma = 1 known)."""
    U = np.fft.fft2(u)
    out = []
    for d, D in zip(discs, fft_discs):
        t = np.real(np.fft.ifft2(U * np.conj(D))) / np.sqrt(d.sum())   # ~ N(0,1) under H0
        out.append(np.log10(N_tests) + norm.logsf(t) / np.log(10))
    return np.array(out)

print("N_tests =", N_tests, " threshold z(eps=1) = %.2f" % norm.isf(1 / N_tests))
for eps in (0.1, 1.0, 10.0):
    tests, groups = [], []
    for _ in range(200):
        det = log10_nfa(rng.normal(size=(n, n))) <= np.log10(eps)
        tests.append(det.sum()); groups.append(label(det.any(axis=0))[1])
    print(f"eps={eps:5}: mean detected tests {np.mean(tests):6.3f}, "
          f"mean detected spots {np.mean(groups):6.3f}")

# a faint spot of radius 3 and contrast a (in units of the noise sigma)
for a in (0.5, 1.0):
    u = rng.normal(size=(n, n)); u[(xx - 64) ** 2 + (yy - 64) ** 2 <= 9] += a
    L = log10_nfa(u)
    near = (xx - 64) ** 2 + (yy - 64) ** 2 <= 25            # within 5 pixels of the spot
    lab, k = label((L <= 0).any(axis=0))
    false = len(set(np.unique(lab[~near])) - set(np.unique(lab[near])) - {0})
    print(f"contrast {a}: best log10 NFA near the spot = {L[:, near].min():5.2f}, "
          f"false spots elsewhere = {false}")

4. Detection of geometric structures

Alignments are the simplest Gestalt and the most frequent structure of man-made scenes. This chapter detects them a contrario, first as alignments of orientations (Desolneux, Moisan and Morel, 2000), then as line segments with the LSD algorithm.

4.1 Meaningful alignments

Consider an image of N \times N pixels and its level-line orientation \theta(x), the gradient direction rotated by \pi/2. Orientations are sampled at points two pixels apart, so that by Proposition 1.4 they are independent on noise.

• A segment S joins two sample points; it contains l sample points.
• A point is aligned with S, at precision p, if its orientation differs from the direction of S by at most p\pi.
• The statistic is the number k of aligned points among the l.
• The background model: the orientations are independent and uniform on the circle, so each point is aligned with probability p.

By Theorem 3.3, fewer than \varepsilon segments are \varepsilon-meaningful on average in a noise image.

The threshold grows only logarithmically with the size of the image: at p = 1/16, doubling the image side adds one point.

4.2 Maximal meaningful segments

When a long alignment is present, many of its sub-segments and super-segments are meaningful too. Following Definition 3.8, a meaningful segment is maximal if it contains no segment with a strictly smaller NFA, and no segment containing it has an NFA smaller than or equal to its own.

Maximal segments locate the endpoints of the alignment: extending a segment into non-aligned points increases its NFA, and cutting aligned points also increases it. The script of § 4.4 reports, for a planted alignment, the segment of smallest NFA, which is the maximal one.

4.3 The LSD line segment detector

LSD (Grompone von Gioi, Jakubowicz, Morel and Randall, 2010) detects line segments in linear time, with the a-contrario validation of § 4.1 adapted to rectangles of pixels.

For a 512 \times 512 image, N_T = 3.9 \times 10^{14}. A perfectly aligned rectangle one pixel wide must contain at least 17 pixels, since \log_{10} N_T / \log_{10} 8 = 16.2. A rectangle of 60 pixels needs at least 34 aligned ones.

LSD counts every pixel of the rectangle, although adjacent 2 \times 2 gradients are dependent (Proposition 1.4, point 3). The experiments of its authors on noise images show essentially no detection, and the method has no parameter to tune for each image.

4.4 Other geometric structures

The same three ingredients, background model, statistic and number of tests, give detectors for many structures.

In each case, most of the work goes into the number of tests: it must count every candidate the algorithm could have examined, including the parameters tuned on the data.

Code 4.6 — Alignments on an orientation grid. The script tests every segment of every row, column and diagonal of a 64 \times 64 grid of orientations, in both directions, at precision p = 1/8. On pure noise, the mean number of 1-meaningful segments over 50 grids is 0.20, below \varepsilon = 1. A planted diagonal segment of 30 sites, with orientations spread by 0.25 rad around the right direction, is detected with \log_{10} \mathrm{NFA} = -13.3 on exactly its 30 sites.

import numpy as np
from scipy.stats import binom

rng = np.random.default_rng(3)
n, p = 64, 1 / 8                     # grid of n x n orientations, precision p
logB = np.full((n + 1, n + 1), 0.0)  # log10 B(l, k, p), lookup table
for l in range(1, n + 1):
    logB[l, :l + 1] = binom.logsf(np.arange(l + 1) - 1, l, p) / np.log(10)

def grid_lines(n):
    """Rows, columns, diagonals and anti-diagonals: (rows, cols) of the sites, direction."""
    i = np.arange(n)
    out = []
    for r in range(n):
        out += [((np.full(n, r), i), 0.0), ((i, np.full(n, r)), np.pi / 2)]
    for d in range(-n + 1, n):
        j = i[(i + d >= 0) & (i + d < n)]
        out.append(((j, j + d), np.pi / 4))            # y - x = -d
    for s in range(2 * n - 1):
        j = i[(s - i >= 0) & (s - i < n)]
        out.append(((j, s - j), 3 * np.pi / 4))        # y + x = s
    return out

LINES = grid_lines(n)
N_tests = sum(2 * len(y) * (len(y) + 1) // 2 for (y, x), _ in LINES)

def meaningful(theta, eps=1.0):
    """List of (log10 NFA, line index, start, end) of eps-meaningful segments."""
    out = []
    for idx, ((y, x), phi) in enumerate(LINES):
        for direction in (phi, phi + np.pi):
            diff = np.angle(np.exp(1j * (theta[y, x] - direction)))
            c = np.concatenate([[0], np.cumsum(np.abs(diff) <= p * np.pi)])
            s, e = np.triu_indices(len(y))            # segment = sites s..e
            l, k = e - s + 1, c[e + 1] - c[s]
            lnfa = np.log10(N_tests) + logB[l, k]
            for t in np.nonzero(lnfa <= np.log10(eps))[0]:
                out.append((lnfa[t], idx, s[t], e[t]))
    return out

print("number of tests N_T =", N_tests)
counts = [len(meaningful(rng.uniform(0, 2 * np.pi, (n, n)))) for _ in range(50)]
print("pure noise: mean number of 1-meaningful segments = %.2f" % np.mean(counts))

theta = rng.uniform(0, 2 * np.pi, (n, n))
j = np.arange(10, 40)                                  # a diagonal segment of 30 sites
theta[j, j] = np.pi / 4 + rng.normal(0, 0.25, j.size)   # ~88% within p*pi = 22.5 deg
best = min(meaningful(theta))
print("planted segment: best log10 NFA = %.1f, sites %d..%d" % (best[0], best[2], best[3]))

5. Anomaly detection in images and hyperspectral data

An anomaly is whatever the model of normal data does not explain: a defect on a textile, a lesion, a rare material in a hyperspectral scene. This chapter reviews the families of anomaly detectors, compares the classical Gaussian detectors of hyperspectral imaging, and shows what the a-contrario approach adds: a threshold that means the same thing on every data set.

5.1 A generic problem and its families of methods

Ehret, Davy, Morel and Delbracio (2019) review image anomaly detectors and group them by the model of normality they use.

• Probabilistic background models: Gaussian, Gaussian mixtures, Markov random fields; the anomaly has a low likelihood.
• Homogeneity and stationarity: the normal image is a stationary texture, described by its spectrum or its local statistics.
• Sparsity: normal patches are sparse in a dictionary learned on the image; anomalies are not.
• Self-similarity: a normal patch has many similar patches elsewhere in the image; an anomaly has few.
• Reconstruction and one-class learning: an auto-encoder or a one-class classifier trained on normal data.

Every family produces an anomaly score per pixel. The decision threshold is the weak point: it is usually tuned on each data set, often with labelled anomalies that detection problems lack.

5.2 The Gaussian background in several bands: RX

In a hyperspectral image, each pixel is a vector x \in \mathbb{R}^d of reflectances in d bands, with d from tens to hundreds.

The a-contrario threshold follows: with N pixels, \mathrm{NFA}(x) = N \cdot \mathbb{P}(\chi^2_d \ge D(x)). For N = 10^4 and d = 10, the threshold is D \ge 35.6, for a mean of 10 under H_0.

The global RX detector uses one mean and one covariance for the whole image. Local RX estimates them in a window around each pixel, which adapts to a changing background but adds estimation noise.

With a window of n = 50 pixels and d = 10 bands, the threshold for a false alarm probability of 10^{-4} is \hat D \ge 62.9, not 35.6. Using the chi-square threshold would give a false alarm rate of about 1 %, a hundred times the nominal one; a simulation of 2 \times 10^4 draws confirms the Fisher law.

5.3 Comparing classical detectors on a hyperspectral example

Two more detectors complete the classical toolbox.

• Matched filter. If the spectral signature t of the target is known, the whitened matched filter of ch. 2.3 applies: T(x) = \frac{(t - \mu)^\top \Sigma^{-1} (x - \mu)}{\sqrt{(t - \mu)^\top \Sigma^{-1} (t - \mu)}} \sim \mathcal{N}(0, 1) \text{ under } H_0 .
• Mixture RX. The background is a mixture of K Gaussian components, one per material, fitted by EM. The statistic is the smallest Mahalanobis distance D_{\min} to a component. If a background pixel comes from component c^*, then D_{\min} \le D_{c^*}, so \mathbb{P}_{H_0}(D_{\min} \ge t) \le \mathbb{P}(\chi^2_d \ge t), and the chi-square NFA remains valid.

The lesson of the comparison is that anomaly is not the same as target. An anomaly detector is only as good as its background model, and the NFA tells when that model is wrong: on a target-free scene, many more than \varepsilon detections mean that the background model must be enriched.

Code 5.4 — RX, mixture RX and matched filter with a-contrario thresholds.

import numpy as np
from scipy.stats import chi2, norm
from sklearn.mixture import GaussianMixture
from sklearn.metrics import roc_auc_score

rng = np.random.default_rng(4)
H = W = 100; d = 10; N = H * W
bands = np.linspace(0, 1, d)
spec = lambda c, w: 0.4 + 0.3 * np.exp(-(bands - c) ** 2 / w)  # smooth reflectance curves
mA, mB, mC, tgt = spec(0.2, 0.05), spec(0.7, 0.08), spec(0.5, 0.3), spec(0.9, 0.02)

cls = np.where(np.arange(W) < W // 2, 0, 1)[None, :].repeat(H, 0)  # two fields
cls[:, 30] = 2                                                     # a road: 1 % of pixels
M = np.stack([mA, mB, mC])[cls]                                    # H x W x d
X = M * (1 + 0.05 * rng.normal(size=(H, W, 1))) + 0.01 * rng.normal(size=(H, W, d))
truth = np.zeros((H, W), bool)
pos = rng.choice(np.flatnonzero(cls.ravel() != 2), 12, replace=False)
truth.flat[pos] = True
X[truth] = 0.85 * X[truth] + 0.15 * tgt                            # sub-pixel targets
X, truth = X.reshape(N, d), truth.ravel()

def rx(X, mu, S):
    Y = X - mu
    return np.einsum("ij,jk,ik->i", Y, np.linalg.inv(S), Y)

def report(name, log10nfa, score):
    det = log10nfa <= 0                                            # eps = 1
    print(f"{name:12s} targets found {np.sum(det & truth):2d}/12, false alarms "
          f"{np.sum(det & ~truth):4d}, AUC {roc_auc_score(truth, score):.3f}")

D = rx(X, X.mean(0), np.cov(X.T))                       # global RX, chi2(d) under H0
report("global RX", np.log10(N) + chi2.logsf(D, d) / np.log(10), D)

K = 3                                                   # Gaussian mixture background
gm = GaussianMixture(K, covariance_type="full", n_init=3, random_state=0).fit(X)
Dk = np.min([rx(X, gm.means_[c], gm.covariances_[c]) for c in range(K)], axis=0)
report("mixture RX", np.log10(N) + chi2.logsf(Dk, d) / np.log(10), Dk)

mu, Si = X.mean(0), np.linalg.inv(np.cov(X.T))          # matched filter, known target
s = tgt - mu
T = (X - mu) @ Si @ s / np.sqrt(s @ Si @ s)
report("matched filt", np.log10(N) + norm.logsf(T) / np.log(10), T)

5.4 The a-contrario contribution: reducing anomalies to detection in noise

Davy, Ehret, Morel and Delbracio (2018) proposed a generic way to give any background model a calibrated threshold.

1. Model the background with a method that reconstructs normal content well and anomalies badly, for instance by replacing each patch with an average of its most similar patches elsewhere in the image (self-similarity), possibly in the feature space of a neural network.
2. Compute the residual, the image minus its reconstruction. Under the background, the residual is close to a structureless noise; normalise it to zero mean and unit variance.
3. Detect a contrario in the residual, as in ch. 3.5: average the residual over discs of several radii; a deviation in either direction is tested with \mathrm{NFA} = N_{\text{pix}}\, N_{\text{scales}} \cdot 2\, \Phi^c(|T|) .

The reconstruction step carries the knowledge of the data; the detection step is always the same and is calibrated by \varepsilon. When the residual is a coloured rather than white noise, it is whitened first, as Grosjean and Moisan (2009) do for spots in textured backgrounds.

6. Detection in time series and recurrent satellite imaging

Satellite constellations now image every place on Earth every few days; the two Sentinel-2 satellites, for instance, revisit the equator every five days. Each pixel becomes a time series, and the questions become temporal: was there a change, when, and of which kind? This chapter builds exact a-contrario tests for changes in seasonal series, screens clouds, and categorises the changes.

6.1 Background models of a time series

Let y = (y_1, \dots, y_T) be the series of a vegetation index at one pixel, for instance the NDVI, the normalised difference of the near-infrared and red reflectances.

• Background. Without change, the series is a seasonal cycle plus noise: y = A\beta + w, where the columns of A \in \mathbb{R}^{T \times q} are a constant and a few harmonics of the yearly period, \beta is unknown, and w \sim \mathcal{N}(0, \sigma^2 I_T) with \sigma unknown.
• Clouds. A cloud perturbs one date strongly; in the NDVI it looks like a large negative outlier.
• Changes. A change adds a profile \delta c to the series, where c is known up to a few parameters (a date, a duration) and \delta is unknown.

The parameters \beta and \sigma are nuisances: the test must be exact whatever their values.

6.2 An exact test for a change

The statistic is the t-statistic of the coefficient of c in the regression on [A, c]. Over N_{\text{pix}} pixels, K families of change profiles and n_f profiles in family f, Proposition 3.4 gives the two-sided NFA \mathrm{NFA} = N_{\text{pix}}\, K\, n_f \cdot 2\, \mathbb{P}(t_{T-q-1} \ge |z|) .

Classical alternatives exist. CUSUM (Page, 1954) detects changes sequentially, PELT (Killick, Fearnhead and Eckley, 2012) segments a series into many pieces, and BFAST (Verbesselt et al., 2010) and CCDC (Zhu and Woodcock, 2014) are standard tools for satellite series. The a-contrario version adds a calibrated threshold across millions of pixel series.

6.3 Families of changes and their categorisation

Three families of profiles describe the usual changes, each indexed by a start date s.

• Abrupt (a harvest, a clear-cut, a construction): c_t = \mathbf{1}_{t \ge s}.
• Transient (a flood, a fire scar that regrows): c_t = \mathbf{1}_{s \le t < s + w} for a few durations w.
• Gradual (a slow degradation): a ramp from s, c_t = \min(1, (t - s)_+ / L).

The factor n_f in the NFA acts as a complexity penalty: a family with more profiles, here the transients with three durations, must explain the data better to be chosen.

6.4 Cloud detection in recurrent satellite imaging

Clouds and their shadows are the main source of false changes. Single-date methods such as Fmask (Zhu and Woodcock, 2012) test spectral properties: clouds are bright, white and cold. Recurrent imaging adds two sources of evidence.

• Temporal evidence. A cloud is a transient bright event: compared with the background built from the other dates, the blue reflectance of a cloudy date is far too high. An a-contrario test with N_{\text{pix}} \times T tests flags it, and a spatial grouping of flagged pixels removes isolated alarms.
• Parallax evidence. Push-broom sensors acquire the spectral bands at slightly different times and angles. Ground objects are superposed in all bands, but objects at altitude, such as clouds, are shifted between bands. Dagobert, Grompone von Gioi, de Franchis, Morel and Hessel (2020) combine luminance and this inter-band parallax.

6.5 Detection and temporal categorisation of changes: a case

Code 6.4 simulates 1000 NDVI series of 100 dates with a seasonal cycle and noise \sigma = 0.04. 850 pixels are stable, and 50 each have an abrupt, a transient and a gradual change of amplitude 0.25. Clouds hit 8 % of the dates.

1. Clouds are screened on the blue band with N_T = 10^5 tests: 8.0 % of the dates are flagged, for 8.0 % truly cloudy.
2. On the remaining dates, Theorem 6.1 tests every profile of the three families, with \varepsilon = 1 for the whole image.

One stable pixel out of 850 is declared changed, consistent with the guarantee of at most one false alarm on average per image. All changes are found and correctly categorised.

Code 6.4 — Cloud screening, change detection and categorisation.

import numpy as np
from scipy.stats import norm, t as student

rng = np.random.default_rng(5)
T, P, npix, sig = 100, 50, 1000, 0.04          # dates, season period, pixels, noise
t = np.arange(T)
kind = np.repeat([0, 1, 2, 3], [850, 50, 50, 50])   # stable, abrupt, transient, gradual
t0 = rng.integers(20, 80, npix)
ndvi = 0.5 + 0.2 * np.sin(2 * np.pi * t / P) + sig * rng.normal(size=(npix, T))
ndvi -= np.where(kind[:, None] == 1, 0.25 * (t >= t0[:, None]), 0)
ndvi -= np.where(kind[:, None] == 2, 0.25 * ((t >= t0[:, None]) & (t < t0[:, None] + 8)), 0)
ndvi -= np.where(kind[:, None] == 3, 0.25 * np.clip((t - t0[:, None]) / 20, 0, 1), 0)
cloud = rng.random((npix, T)) < 0.08
ndvi[cloud] -= 0.5                                  # clouds darken NDVI ...
blue = 0.05 + 0.005 * rng.normal(size=(npix, T)) + 0.3 * cloud   # ... and brighten blue

# 1. cloud screening: a-contrario test of a bright blue value, N_T = npix * T tests
zb = (blue - np.median(blue, 1, keepdims=True)) / 0.005
valid = np.log10(npix * T) + norm.logsf(zb) / np.log(10) > 0

# 2. change families: contrast vectors c (one per row), with their family index
C, fam = [], []
for s in range(1, T):
    C.append(t >= s); fam.append(0)                              # abrupt
    C.append(np.clip((t - s) / 20, 0, 1)); fam.append(2)         # gradual (ramp)
    for w in (4, 8, 16):
        C.append((t >= s) & (t < s + w)); fam.append(1)          # transient
C, fam = np.array(C, float), np.array(fam)
n_f = np.bincount(fam)                                           # tests per family
A = np.stack([np.ones(T), np.sin(2 * np.pi * t / P), np.cos(2 * np.pi * t / P)], 1)

pred = np.zeros(npix, int)                                       # 0 = no change
for i in range(npix):
    v = valid[i]
    Q = np.eye(v.sum()) - A[v] @ np.linalg.pinv(A[v])            # residual projector
    y, Rc = ndvi[i, v], C[:, v] @ Q
    nc = np.maximum(np.linalg.norm(Rc, axis=1), 1e-9)
    proj = Rc @ y / nc                                           # component of y along Qc
    dof = v.sum() - 4                                            # model [A, c] has 4 columns
    s2 = (np.sum((Q @ y) ** 2) - proj ** 2) / dof                # residual variance of [A, c]
    z = proj / np.sqrt(s2)                                       # Student t(dof) under H0
    lnfa = np.log10(npix * 3 * n_f[fam] * 2) + student.logsf(np.abs(z), dof) / np.log(10)
    lnfa[nc < 1e-6] = np.inf
    if lnfa.min() <= 0:
        pred[i] = 1 + fam[np.argmin(lnfa)]                       # 1 abrupt, 2 transient, 3 gradual

print("cloudy dates flagged: %.1f %% (true %.1f %%)" % (100 * (~valid).mean(), 100 * cloud.mean()))
names = ["stable", "abrupt", "transient", "gradual"]
conf = np.zeros((4, 4), int)
np.add.at(conf, (kind, pred), 1)
print("rows = truth, columns = decision (none, abrupt, transient, gradual)")
for k in range(4):
    print(f"{names[k]:9s}", conf[k])

7. Detection of image forgeries

A forged photograph usually looks right: the semantic content gives nothing away. What betrays it are the traces of the camera and of the processing chain, which the forgery breaks locally. This chapter detects such inconsistencies a contrario, where a false alarm means accusing an authentic image.

7.1 Forensic traces

Common forgeries are splicing (a region from another image), copy-move (a region duplicated within the image), removal by inpainting, and, increasingly, regions produced by generative models. The traces used to detect them come from the acquisition pipeline.

• Demosaicking. Each sensor pixel measures one colour of a periodic colour filter array; the two other colours are interpolated.
• Noise. The noise level depends on the intensity through a camera-specific curve.
• JPEG compression. The image is cut into 8 \times 8 blocks on a fixed grid, and the DCT coefficients of each block are quantised.

Each forgery type has too many variants to be learned from examples. The a-contrario approach takes the traces of the image itself as background and detects regions whose traces are too inconsistent to be due to chance.

7.2 JPEG grid detection

JPEG quantises the DCT coefficients of each block with a table Q; decoding applies the inverse DCT and rounds to integers. Many high-frequency coefficients are quantised to zero.

A block that is not aligned with the grid straddles two to four JPEG blocks; in textured regions, its DCT has many coefficients above 0.5. Nikoukhah, Anger, Ehret, Colom, Morel and Grompone von Gioi (2019) build a detector on this contrast.

The 64 pixels of an 8 \times 8 cell share most of their blocks and vote together, which violates the independence of the background model. Code 7.4 therefore keeps one vote per cell; Exercise 14 shows what happens otherwise.

7.3 Demosaicking traces

In a Bayer filter array, the colours repeat with period 2 in each direction, so the array has four possible positions relative to the image. Demosaicking leaves the measured samples and the interpolated ones with different statistics, and the position can be estimated locally. Popescu and Farid (2005) estimate the interpolation correlations with an EM algorithm; Bammey, Grompone von Gioi and Morel (2020) train a small network on the image itself to recognise the position of each pixel in the pattern.

The a-contrario step is the one of Definition 7.2 with four positions instead of 64: a region whose votes favour another position, or no position, is detected with a binomial NFA.

7.4 Noise and copy-move traces

• Noise. Blocks of similar intensity should have similar noise levels. For a block of m residual samples of variance \sigma^2, (m-1)s^2/\sigma^2 follows a \chi^2_{m-1} law, which gives the NFA of an abnormally low or high noise level. Gardella, Musé, Morel and Colom (2021) define a background model of the noise and detect local noise anomalies by their NFA.
• Copy-move. Matching patches within the image produces many spurious matches, since natural images are self-similar. A copy-move is a group of matches with the same displacement; the background model draws displacements uniformly, and the number of coherent matches in a region receives a binomial NFA (Ehret, 2019).

7.5 Evaluation

A forensic tool is judged on two numbers: the localisation quality on forged images, and the number of false detections on authentic ones. The NFA bounds the second by \varepsilon per image, under the background model.

Bammey and coauthors (2022) built a benchmark of forgeries that cannot be detected from the semantic content alone, so that a method is evaluated on the traces it claims to use. Deep forensic detectors are powerful, but most return scores without a false-alarm guarantee (ch. 8).

Code 7.4 — JPEG grid votes and a-contrario detection of a spliced region.

import numpy as np
from scipy.fft import dct, idct
from scipy.ndimage import gaussian_filter
from scipy.stats import binom
from scipy.special import xlogy
from numpy.lib.stride_tricks import sliding_window_view

rng = np.random.default_rng(6)
Q50 = np.array([[16,11,10,16,24,40,51,61],[12,12,14,19,26,58,60,55],
                [14,13,16,24,40,57,69,56],[14,17,22,29,51,87,80,62],
                [18,22,37,56,68,109,103,77],[24,35,55,64,81,104,113,92],
                [49,64,78,87,103,121,120,101],[72,92,95,98,112,100,103,99]])
Q = np.maximum(np.round(Q50 * 0.5), 1)                 # quality 75 luminance table
dct2 = lambda b: dct(dct(b, axis=-1, norm="ortho"), axis=-2, norm="ortho")
idct2 = lambda b: idct(idct(b, axis=-1, norm="ortho"), axis=-2, norm="ortho")

def scene(n):
    u = gaussian_filter(rng.normal(size=(n, n)), 3) * 400 + 128
    return np.clip(u + 8 * rng.normal(size=(n, n)), 0, 255)   # smooth content + texture

def jpeg(u):
    """Quantise 8x8 block DCTs on the grid anchored at (0, 0); decode to integers."""
    n = u.shape[0]
    B = (u - 128).reshape(n // 8, 8, n // 8, 8).transpose(0, 2, 1, 3)
    B = idct2(np.round(dct2(B) / Q) * Q)
    return np.clip(np.round(B.transpose(0, 2, 1, 3).reshape(n, n) + 128), 0, 255)

def votes(v):
    """Each pixel votes for the grid origin of the 8x8 block (among the 64 containing it)
    with the most DCT coefficients rounding to zero; ties give no vote (-1)."""
    n = v.shape[0]
    Z = (np.abs(dct2(sliding_window_view(v - 128, (8, 8)))) < 0.5).sum((-1, -2))
    m = Z.shape[0]
    stack = np.full((64, n, n), -1)
    for a in range(8):
        for b in range(8):                        # pixel at offset (a, b) inside the block
            stack[8 * a + b, a:a + m, b:b + m] = Z
    best = stack.max(0)
    arg = stack.argmax(0)
    x, y = np.mgrid[:n, :n]
    origin = 8 * ((x - arg // 8) % 8) + (y - arg % 8) % 8
    return np.where((stack == best).sum(0) == 1, origin, -1)

n = 256
img = jpeg(scene(n))                                   # authentic image, grid origin (0, 0)
other = jpeg(scene(n + 8))[3:3 + n, 5:5 + n]           # another JPEG, grid origin (5, 3)
forged = img.copy(); forged[96:160, 64:144] = other[96:160, 64:144]

def log10_tail(k, N, p):
    """log10 B(N, k, p), with the Chernoff bound when the exact value underflows."""
    exact = binom.logsf(k - 1, N, p) / np.log(10)
    if np.isfinite(exact):
        return exact
    a = k / N                                          # xlogy(0, .) = 0 covers k = N
    return -N * (xlogy(a, a / p) + xlogy(1 - a, (1 - a) / (1 - p))) / np.log(10)

for name, v in (("authentic", img), ("forged", forged)):
    V = votes(v)[4::8, 4::8]                           # one vote per 8 x 8 cell
    main = np.bincount(V[V >= 0], minlength=64).argmax()
    k, N = np.sum(V == main), np.sum(V >= 0)
    print(f"{name}: main grid origin {divmod(int(main), 8)}, {k}/{N} votes, "
          f"log10 NFA <= {np.log10(64) + log10_tail(k, N, 1 / 64):.0f}")
    W = 4                                              # windows of 4 x 4 cells (32 x 32 pixels)
    m = V.shape[0] - W + 1
    n_tests = m * m * 63
    hits = []
    for i in range(m):
        for j in range(m):
            w = V[i:i + W, j:j + W]
            c = np.bincount(w[w >= 0], minlength=64); c[main] = 0
            g = c.argmax()
            if c[g] and np.log10(n_tests) + log10_tail(c[g], np.sum(w >= 0), 1 / 64) <= 0:
                hits.append((8 * i, 8 * j, divmod(int(g), 8)))
    apart = sum(not (i < 160 and i + 32 > 96 and j < 144 and j + 32 > 64) for i, j, _ in hits)
    print(f"   meaningful foreign-grid windows: {len(hits)} (of {m * m}), disjoint from the"
          f" pasted region: {apart}, origins {sorted({h[2] for h in hits})}")

8. Research directions and open problems

This chapter is for self-study. It surveys questions on which current work focuses, with pointers to papers; the statements about each work summarise its abstract-level claims.

8.1 Dependence and the number of tests

Theorem 3.3 holds whatever the dependence between tests, but each tail is computed under independence inside a test. Real features are correlated: adjacent gradients (Proposition 1.4), grid votes of one cell (ch. 7.2), residuals of a reconstruction. When the dependence is positive, the true tail is typically heavier than the binomial one, and false alarms exceed \varepsilon (Exercise 13).

Myaskouvskey, Gousseau and Lindenbaum (2013) extend the a-contrario decision to dependent features. Two questions remain open in general.

• What is the effective number of independent observations inside a test, for a given image model?
• Conversely, when the tests overlap heavily, Bonferroni counting is pessimistic: Code 4.6 finds 0.20 false segments on average for \varepsilon = 1. Random field theory, used in brain imaging since Worsley and coauthors (1996), counts excursions of smooth Gaussian fields more tightly; its use in a-contrario detection is little explored.

8.2 Learned background models

The reduction of ch. 5.4 separates a reconstruction step, which carries the knowledge of the data, from a calibrated detection step. Davy, Ehret, Morel and Delbracio (2018) apply it with self-similarity in the feature space of pretrained networks. In forensics, Bammey, Grompone von Gioi and Morel (2020) train a network on the very image under analysis.

The open point is the hypothesis of the detection step: the residual of a learned model must be structureless noise on normal data. No general test of this hypothesis exists; when it fails, the NFA is no longer calibrated.

When the background parameters are estimated on the image itself (a mixture model, a noise level), the tests become data-dependent. The guarantee then holds only approximately, a question studied in statistics under the name of selective inference.

8.3 Conformal p-values and the NFA of learned scores

Any anomaly score can be turned into a super-uniform p-value by calibration on normal data.

Bates, Candès, Lei, Romano and Sesia (2023) study these p-values for outlier detection, including their dependence through the shared calibration set.

Their limitation for images is the smallest attainable value: p \ge 1/(n+1), so \mathrm{NFA} \ge N / (n+1). Testing 10^6 pixels at \varepsilon = 1 needs at least 10^6 normal calibration samples. Combining conformal calibration with parametric or extreme-value tails, as in Statistical Learning with Extreme Values, is an open direction.

Code 8.2 — An NFA for an isolation forest. With 5000 calibration samples and 1020 test points, the smallest attainable NFA is 0.204. The detector finds 18 of the 20 anomalies with no false alarm among the 1000 normal points.

import numpy as np
from sklearn.ensemble import IsolationForest

rng = np.random.default_rng(8)
d = 5
train, calib = rng.normal(size=(2000, d)), rng.normal(size=(5000, d))
test = np.vstack([rng.normal(size=(1000, d)),                 # 1000 normal points
                  rng.normal(size=(20, d)) + 2.5])            # 20 anomalies, shifted mean
is_anom = np.r_[np.zeros(1000, bool), np.ones(20, bool)]

model = IsolationForest(n_estimators=200, random_state=0).fit(train)
s_cal = -model.score_samples(calib)                           # large = anomalous
s_test = -model.score_samples(test)

# conformal p-value: super-uniform as soon as calibration and test normals are exchangeable
p = (1 + (s_cal[None, :] >= s_test[:, None]).sum(1)) / (len(s_cal) + 1)
nfa = len(test) * p
det = nfa <= 1
print("smallest attainable NFA: %.3f" % (len(test) / (len(s_cal) + 1)))
print(f"detections: {det[is_anom].sum()}/20 anomalies, {det[~is_anom].sum()} false alarms")

8.4 Forensics in the age of generative models

Images produced or edited by diffusion models are a new kind of forgery. Corvi, Cozzolino, Zingarini, Poggi, Nagano and Verdoliva (2023) studied detectors of images generated by diffusion models. Guillaro, Cozzolino, Sud, Dufour and Verdoliva (2023) proposed TruFor, a network that localises many kinds of forgeries and outputs a reliability map.

Two problems stay open: generalisation to generators unseen in training, and calibrated false-alarm rates on authentic images. Whether the traces of generative pipelines, such as up-sampling layers and auto-encoder decoders, can serve as background models for a-contrario tests is an active question.

8.5 Deep detectors and a-contrario validation

Learned detectors now compete with LSD and with classical change detection: Huang and coauthors (2018) parse wireframes of man-made scenes with a network, and Daudt, Le Saux and Boulch (2018) detect changes between satellite images with fully convolutional siamese networks. Their outputs are scores. Validating their proposals a contrario, by an NFA computed on the image measurements, would give them the false-detection control they lack.

8.6 Open problems

1. Compute the effective number of tests of a detector under a realistic dependence model.
2. Choose the weights of Proposition 3.4 optimally across families of candidates (lengths, scales, change types).
3. Reach small NFAs with finite calibration sets for learned scores.
4. Give calibrated false-alarm guarantees for forensic detection of generated content.
5. Test, on normal data, the hypothesis that a learned residual is structureless.

Exercises

1. ★ Binomial tails and their bounds. For n = 50 observations with p = 1/8, compute \mathcal{B}(50, 20, 1/8) exactly, then its Chernoff and Hoeffding bounds (Proposition 0.2). Which bound would you use to explain a threshold, and which to compute it?

   Solution: the exact tail is 9.3 \times 10^{-7}. The Chernoff bound, with a = 0.4, is e^{-50\, \mathrm{KL}(0.4 \| 0.125)} = 6.5 \times 10^{-6} , and Hoeffding gives e^{-100 \times 0.275^2} = 5.2 \times 10^{-4}.

   Thresholds are computed with the exact tail. The Chernoff bound, within a factor 7 here, explains how they scale; Hoeffding is too loose because it ignores the small variance np(1-p).

2. ★ A discrete p-value. Let X \sim \mathrm{Binomial}(2, 1/2) and S(x) = \mathbb{P}(X \ge x). Compute S at 0, 1 and 2, then \mathbb{P}(S(X) \le t) for t \in \{0.25, 0.3, 0.8\}. Check Proposition 0.7 and say when equality holds.

   Solution: S(0) = 1, S(1) = 3/4, S(2) = 1/4. Then \mathbb{P}(S(X) \le 0.25) = \mathbb{P}(X = 2) = 0.25, with equality; \mathbb{P}(S(X) \le 0.3) = 0.25 < 0.3; and \mathbb{P}(S(X) \le 0.8) = \mathbb{P}(X \ge 1) = 0.75 < 0.8.

   Equality holds exactly at the values taken by S; in between, the p-value is conservative.

3. ★ Gradients of white noise. With the notation of Proposition 1.4, compute the covariance between g_x at pixel (i, j) and g_y at pixel (i, j+1). Then explain why orientations sampled every two pixels are independent, while orientations at adjacent pixels are not.

   Solution: the two windows share two samples.

   • The sample b = u(i, j+1) enters g_x(i, j) with coefficient +1/2 and g_y(i, j+1) with coefficient -1/2.
   • The sample d = u(i+1, j+1) enters both with coefficient +1/2.

   The other samples appear in only one of the two components, so the covariance is \sigma^2 \big(-\tfrac14 + \tfrac14\big) = 0 .

   These two components are uncorrelated, but g_x at adjacent pixels has correlation -1/2, so adjacent orientations are dependent. Windows two pixels apart are disjoint, and the gradients are functions of disjoint sets of independent samples.

4. ★ Bonferroni in numbers. A 1000 \times 1000 image is tested pixel by pixel with Gaussian statistics. Give the Bonferroni thresholds for \varepsilon = 1 and \varepsilon = 0.05, and the expected number of false alarms of a per-pixel threshold at 3\sigma. What does \varepsilon = 1 guarantee on the FWER?

   Solution: z = (\Phi^c)^{-1}(10^{-6}) = 4.75 for \varepsilon = 1 and (\Phi^c)^{-1}(5 \times 10^{-8}) = 5.33 for \varepsilon = 0.05. A 3\sigma threshold gives 10^6 \times 1.35 \times 10^{-3} = 1350 false alarms.

   For \varepsilon = 1, the FWER bound \mathbb{P}(V \ge 1) \le 1 is empty, but \mathbb{E}[V] \le 1 remains a useful guarantee (Proposition 2.7).

5. ★ Gaussian thresholds and minimal contrasts. With N_T = 65\,536 as in Example 3.7, compare z(1) = 4.17 with the bound of Proposition 3.6. Then compute the minimal contrast a_{\min} for discs of radius 2 (13 pixels) and 4 (49 pixels).

   Solution: \sqrt{2 \log 65\,536} = 4.71, an upper bound 0.54 above the exact threshold. The minimal contrasts are 4.17 / \sqrt{13} = 1.16 and 4.17 / 7 = 0.60 noise standard deviations: larger objects can be fainter.

6. ★ Minimal length of an alignment. For an image of 1024 \times 1024 pixels, compute the minimal number of perfectly aligned points of a meaningful alignment for p = 1/16 and p = 1/8 (Proposition 4.2). For LSD on the same image, compute the minimal length of a perfectly aligned rectangle one pixel wide.

   Solution: \log(1024^4) / \log 16 = 10 points and \log(1024^4) / \log 8 = 13.3, so 14 points. For LSD, \log_{10} N_T = 2.5 \log_{10}(2^{20}) + \log_{10} 11 = 16.09 , and 16.09 / \log_{10} 8 = 17.8: 18 pixels, one more than for a 512 \times 512 image.

7. ★ RX thresholds. A hyperspectral image has N = 10^6 pixels and d = 10 bands. Give the RX threshold on D for \varepsilon = 1 with known background statistics, then with statistics estimated in windows of n = 50 pixels (Proposition 5.2).

   Solution: the chi-square threshold is \chi^2_{10} at tail 10^{-6}, that is D \ge 46.9. With estimated statistics, the Fisher law at tail 10^{-6} gives \hat D \ge 96.6 . The second threshold is twice the first: a 50-pixel window estimates a 10 \times 10 covariance poorly.

8. ★ The smallest detectable step. A series of T = 60 dates with constant background has a step after date 20, noise \sigma = 0.03, and a detection threshold |z| \ge 5. Compute \|Qc\| and the smallest step detected half of the time.

   Solution: Qc = c - \bar c, so \|Qc\|^2 = n_a n_b / T with n_b = 20 dates before and n_a = 40 after: \|Qc\| = \sqrt{800/60} = 3.65. The mean statistic is \delta \|Qc\| / \sigma, so \delta_{\min} = \frac{5 \times 0.03}{3.65} = 0.041 .

9. ★★ The matched filter in coloured noise. Let x = as + w with w \sim \mathcal{N}(0, C). Prove that the Neyman–Pearson test thresholds s^\top C^{-1} x, and give its power at level \alpha. Show that among all linear statistics \langle h, x \rangle, the ratio of the mean shift to the standard deviation is maximal for h \propto C^{-1}s.

   Solution: the log-likelihood ratio is a\, s^\top C^{-1} x - \frac{a^2}{2} s^\top C^{-1} s. With d = a \sqrt{s^\top C^{-1} s}, the power is \Phi^c(z_\alpha - d), for example 0.82 for \alpha = 10^{-3} and d = 4.

   For a linear statistic, the ratio is a \langle h, s \rangle / \sqrt{h^\top C h}. Writing \langle h, s \rangle = \langle C^{1/2} h, C^{-1/2} s \rangle, Cauchy–Schwarz bounds it by a \sqrt{s^\top C^{-1} s}, with equality for C^{1/2} h \propto C^{-1/2} s.

10. ★★ Weighted tests by length. Candidates are segments with lengths l = 1, \dots, L, and N_l segments have length l. Show that, for a segment S of length l with k aligned points, \mathrm{NFA}(S) = L\, N_l\, \mathcal{B}(l, k, p) is an NFA. When is it smaller than the unweighted NFA N_T \mathcal{B}(l, k, p), where N_T = \sum_l N_l?

    Solution: apply Proposition 3.4 with weights w_S = 1/(L N_l), whose sum over all segments is \sum_l N_l / (L N_l) = 1.

    The weighted NFA is smaller exactly when L N_l < N_T, that is when length l has fewer candidates than the average length.

    On the grid lines of Code 4.6, a line of n sites holds n - l + 1 segments of length l. Long segments are therefore rarer than short ones, and weighting favours long alignments, whose detection matters most.

11. ★★ Why Hoeffding fails for clusters. In Example 3.5, write the Chernoff threshold as the solution of n\, \mathrm{KL}(k/n \| a) = \log N_T and the Hoeffding threshold explicitly. Explain the gap between 27.8 and 83.6 with the variance of the count.

    Solution: Hoeffding gives k \ge na + \sqrt{\frac{n}{2} \log N_T} = 10 + 73.6 = 83.6. The Chernoff equation, solved numerically, gives k = 27.8, close to the exact 26.

    Hoeffding’s bound holds for any variable in [0, 1] and therefore uses the worst variance 1/4 per observation, 250 in total. The actual variance of the count is na(1-a) = 9.9, so deviations are about five times smaller than Hoeffding assumes.

12. ★★ The bias of local RX. In the setting of Proposition 5.2, prove that x - \bar y \sim \mathcal{N}(0, (1 + 1/n)\Sigma) and is independent of S. Then use the moment of the inverse Wishart matrix, valid for n > d + 2, \mathbb{E}[S^{-1}] = \frac{n-1}{n-d-2}\, \Sigma^{-1} , to compute \mathbb{E}[\hat D], and evaluate it for n = 50 and d = 10.

    Solution: x - \bar y is a linear combination of independent Gaussian vectors with covariance \Sigma + \Sigma / n. The sample mean and covariance of a Gaussian sample are independent, and x is independent of the sample. Conditioning on S, \mathbb{E}[\hat D] = \Big(1 + \frac1n\Big) \mathrm{tr}\big(\mathbb{E}[S^{-1}]\Sigma\big) = \Big(1 + \frac1n\Big) \frac{d\,(n-1)}{n - d - 2} . For n = 50 and d = 10, this is 13.2 instead of d = 10: the estimated distance is inflated, hence the larger threshold.

13. ★★ Alignments on dependent orientations (Python). Modify Code 4.6 to estimate the mean number of meaningful segments on noise for \varepsilon \in \{0.1, 1, 10\}. Then replace the independent orientations by the level-line orientations of the 2 \times 2 gradient of a white noise image, computed (a) at every pixel and (b) at every other pixel. Explain the results.

    Hints and expected result: use the level-line angle, the gradient angle plus \pi/2. Average over at least 200 grids, because one chance alignment produces many overlapping meaningful segments and the counts are heavy-tailed. Typical means for the three values of \varepsilon (they vary by about 30 % with the seed) are:

    • independent orientations: about 0.05, 0.5 and 4 segments, below \varepsilon each time;
    • orientations at every pixel: about 1.5, 8 and 45, well above \varepsilon each time;
    • orientations at every other pixel: about 0, 0.06 and 1.6, below \varepsilon again.

    At every pixel the guarantee fails, because the tail of each test is computed as if adjacent orientations were independent. Along a row, a horizontal segment tests the sign and size of g_y. The computation of Proposition 1.4, point 3, shows that g_y has correlation +1/2 between adjacent windows of a row, and this positive dependence makes the tail heavier.

    With the raw gradient angle, the same segments test g_x, of correlation -1/2, and almost nothing is detected. On a noise image smoothed by a Gaussian of standard deviation 1, a grid typically holds tens of 1-meaningful segments, sometimes more than a hundred: smoothing creates genuine local alignments, which the background model excludes.

14. ★★ Dependent votes in JPEG grid detection (Python). In Code 7.4, use all pixel votes instead of one vote per 8 \times 8 cell, with windows of 32 \times 32 pixels at stride 8. Count the meaningful foreign-grid windows on the authentic image and explain.

    Hints and expected result: false detections appear on the authentic image, voting for origins adjacent to the main grid such as (0, 7). The 64 pixels of a cell share most of their blocks and vote almost identically, so a window of 1024 votes holds about 16 independent ones. The binomial tail with n = 1024 then underestimates the probability of a cluster of foreign votes.

15. ★★★ An effective number of observations for LSD. LSD counts every pixel of a rectangle, although adjacent gradients are dependent. Propose a corrected count, for instance using only a sub-grid of step 2 or an effective number n_{\text{eff}} < n. Implement it on noise images and on synthetic segments, and compare the false detections and the shortest detectable segment with the original count.

    Hints: on noise, the sub-grid count is calibrated (Exercise 13) but needs about twice the length to detect a segment. An intermediate count, such as a binomial on n/c observations, can be calibrated by simulation. Discuss the trade-off with the experiments on noise reported by the LSD authors and with Myaskouvskey, Gousseau and Lindenbaum (2013).

16. ★★★ An NFA from a learned score. With Code 8.2: (a) explain why the NFA stays valid when the isolation forest is badly trained; (b) compute how many calibration samples are needed to test every pixel of a 512 \times 512 image at \varepsilon = 1; (c) propose a way to reach smaller NFAs, for instance by fitting a generalised Pareto tail to the largest calibration scores, and say which guarantee is lost.

    Hints: (a) Proposition 8.1 only needs exchangeability; a bad score loses power, not calibration. (b) N / (n+1) \le 1 requires n \ge 262\,143. (c) A parametric tail extrapolates beyond the calibration range; the guarantee becomes asymptotic and depends on the tail model, which can itself be tested on held-out normal data.

Projects and assessment

The official assessment is continuous. master-mva.com announces 4 homework assignments, including exercises and a report analysing a research paper; the Université Paris-Saclay detail adds experimental reports based on IPOL articles with an online demo, and at least one in-depth analysis of a detection paper. No weighting is displayed.

The subjects below are proposals for that report and for personal training, in the MVA tradition of studying and reimplementing a research paper. The written exam is not part of the official assessment: it shows what a timed test on these notes could look like, and serves as a self-test.

Solo projects

1. [Solo · ≈ 20 h] Paper study: the LSD line segment detector (Grompone von Gioi, Jakubowicz, Morel and Randall, 2010; IPOL 2012). Does a reimplementation keep the false-detection control that the paper claims, and what does the dependence of adjacent gradients cost?
   • Implement: a simplified LSD in NumPy (Algorithm 4.4): the 2 \times 2 gradient, region growing at \tau = 22.5°, the rectangle fit and the NFA of Definition 4.5. Compare it with the reference C code published with the IPOL article.
   • Data: 100 white-noise images of 512 \times 512 pixels; synthetic images of segments of controlled length and contrast; the sample photographs bundled with scikit-learn (sklearn.datasets.load_sample_image).
   • Success criterion: at most one detection per noise image on average; on noiseless synthetic images, a shortest detected segment within two pixels of the bound of ch. 4.3 (17 pixels).
   • Deliverables and evaluation: code, a 4–6 page report in the format of an IPOL experimental report, a 10-minute talk; rubric below. Chapters used: ch. 1, 3, 4.
2. [Solo · ≈ 20 h] Paper study: JPEG grid detection (Nikoukhah, Anger, Ehret, Colom, Morel and Grompone von Gioi, 2019). How small and how compressed can a spliced region be and still be found, with no false accusation of authentic images?
   • Implement: Code 7.4 with real JPEG compression (the Pillow library installed with Matplotlib), one vote per 8 \times 8 cell, and the main-grid and foreign-grid NFAs of Definition 7.2.
   • Data: crops of the scikit-learn sample images, compressed at qualities 60 to 95; splices whose grid origin differs from the host’s, with sizes from 32 \times 32 to 128 \times 128 pixels.
   • Success criterion: fewer than one meaningful foreign window per authentic image on average over 100 crops; a detection rate plotted against quality and region size.
   • Deliverables and evaluation: code, report, talk; rubric below. Chapters used: ch. 0, 3, 7.
3. [Solo · ≈ 15 h] Paper study: stability of multiple testing (Gordon, Glazko, Qiu and Yakovlev, 2007). Does the Bonferroni rule keep its advantage in stability when the tests are as correlated as image tests?
   • Implement: Code 2.8 with three correlation structures: equicorrelation, block correlation, and the spatial correlation of a smoothed Gaussian random field tested pixel by pixel.
   • Prove: \mathbb{E}[V] \le \varepsilon N_0 / N under any dependence (Proposition 2.5), and show by an example that the FDR bound of Theorem 2.6 can fail under negative dependence (two null p-values, coupled so that \{p_1 \le q/2\} and \{p_2 \le q/2\} are disjoint and both lie in (q/2, q] with probability q/2, give an FDR of 3q/2).
   • Success criterion: the paper’s qualitative finding reproduced, a smaller variance of the numbers of true and of total discoveries for Bonferroni, with intervals over 4000 runs.
   • Deliverables and evaluation: code, a 4-page report, a 10-minute talk; rubric below. Chapters used: ch. 0, 2, 3.

Group projects

1. [Group of 2–3 · ≈ 30 h each] Change detection and categorisation in satellite image time series. Can the exact tests of ch. 6 date and classify changes on real series, with a false-alarm level that matches \varepsilon?
   • Implement: temporal cloud screening; the Student test of Theorem 6.1 with yearly harmonics; the families and the categorisation of Definition 6.3; a CUSUM baseline (Page, 1954) whose threshold is tuned on the simulator.
   • Data: the simulator of Code 6.4, with known truth; a public Sentinel-2 Level-2A time series of one tile, downloaded from the Copernicus Data Space Ecosystem, around a documented event (a harvest, a clear-cut, a flood).
   • Success criterion: on the simulator, the number of false changes per 10^6 stable series stays below 10^6 \varepsilon / N_{\text{pix}} for three values of \varepsilon, and the categorisation accuracy is reported against the change amplitude.
   • Split: member 1 builds cloud screening and data loading; member 2 builds the tests and the NFA; member 3 runs the evaluation and the baseline.
   • Deliverables and evaluation: code, a joint 6–8 page report, a 15-minute talk with an individual question per member; rubric below. Chapters used: ch. 3, 6.
2. [Group of 2–3 · ≈ 30 h each] Paper study: reducing anomaly detection in images to detection in noise (Davy, Ehret, Morel and Delbracio, 2018). Does a self-similarity reconstruction leave a residual on which the NFA is calibrated?
   • Implement: the three steps of ch. 5.4 — a non-local patch reconstruction that excludes a neighbourhood of each patch, the normalisation of the residual, and the multiscale two-sided NFA. Compare with global RX on 5 \times 5 patch vectors.
   • Data: synthetic textures (periodic patterns, filtered noise) with planted defects of controlled size and contrast; defect-free crops of the scikit-learn sample images.
   • Success criterion: on defect-free textures, the mean number of detections stays within a factor 2 of \varepsilon; on textures with defects, a detection-rate curve against the defect contrast.
   • Split: reconstruction; detection and calibration (including a test of the whiteness of the residual); evaluation and comparison with RX.
   • Deliverables and evaluation: code, joint report, talk; rubric below. Chapters used: ch. 3, 5.
3. [Group of 3 · ≈ 30 h each] Mosaic consistency as a forensic trace (after Popescu and Farid, 2005, and Bammey, Grompone von Gioi and Morel, 2020). Can a simple interpolation residual locate a region whose colour filter array is shifted, with a binomial NFA over the four positions?
   • Implement: a synthetic camera pipeline (Bayer mosaic at a chosen offset, bilinear demosaicking, optional JPEG compression); a per-block vote for the most likely of the four mosaic positions; the NFA of ch. 7.3.
   • Data: the scikit-learn sample images, mosaicked and demosaicked; splices with a shifted mosaic.
   • Success criterion: fewer than one false region per authentic image on average; the size of the smallest detected splice, before and after JPEG compression at quality 90.
   • Split: pipeline and data; detector and NFA; evaluation, robustness and comparison with the JPEG-grid detector of Code 7.4.
   • Deliverables and evaluation: code, joint report, talk; rubric below. Chapters used: ch. 3, 7.

Written exam (proposal)

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

1. Question (5 points). Define a background model and the number of false alarms of a family of N_T tests. State and prove the fundamental property of the NFA. Explain where independence is used and where it is not, and what a detection means when the background model is wrong.

   Expected:

   • Definitions 3.1 and 3.2, with the property that the expected number of \varepsilon-meaningful candidates is at most \varepsilon (1 pt).
   • The proof of Theorem 3.3: linearity of expectation over the indicators, then the super-uniformity of each S_i(X_i) (2 pts).
   • Independence is needed inside a test, to compute the tail; it is not needed between tests (1 pt).
   • A detection is a local rejection of the background model, not the proof that an object exists (1 pt).

2. Question (5 points). An image has 256 \times 256 pixels, and alignments are tested at precision p = 1/16 with N_T = N^4 segments. (a) Find the minimal length of a perfectly aligned 1-meaningful segment. (b) For a segment of l = 20 independent points, give the number of aligned points that Hoeffding’s bound declares sufficient. (c) The exact threshold is 13 aligned points; explain the gap.

   Expected:

   • 1. N^4 = 2^{32} and p^l = 2^{-4l}, so l \ge 8 points; with 8 points the NFA equals 1 exactly (2 pts).
   • 2. Proposition 3.6 gives k \ge 1.25 + \sqrt{10 \log 2^{32}} = 1.25 + 14.9, so 17 points (1.5 pts).
   • 3. Hoeffding uses the worst variance 1/4 per point, while the true variance p(1-p) \approx 0.06 is four times smaller; the exact tail or the Chernoff bound must be used when p is small (1.5 pts).

3. Question (5 points). (a) Prove that for every real random variable X with survival function S, the variable S(X) is super-uniform. (b) Let w_i > 0 be weights of sum at most 1. Prove that \mathrm{NFA}_w(i, x) = \mathbb{P}_{H_0}(X_i \ge x_i) / w_i is an NFA. (c) Deduce an NFA for K families of candidates of sizes N_1, \dots, N_K.

   Expected:

   • 1. The set where S \le t is a half-line; the closed case gives a probability S(x^*) \le t, and the open case a limit of values at most t (Proposition 0.7) (2 pts).
   • 2. Linearity of expectation and super-uniformity give at most \varepsilon \sum_i w_i \le \varepsilon false detections (2 pts).
   • 3. The weights w_i = 1/(K N_k) give the NFA K N_k\, \mathbb{P}_{H_0}(X_i \ge x_i) for a candidate of family k (1 pt).

4. Question (5 points). A vegetation series has T = 40 dates, a constant background of unknown level, Gaussian noise of unknown standard deviation \sigma, and a possible step after date 10. (a) Give the law of the test statistic of Theorem 6.1 under the background model, with its number of degrees of freedom. (b) Compute \|Qc\|. (c) The image holds 10^5 tests, and the threshold for \varepsilon = 1 is |z| \ge 5.1 . Find the smallest step detected half of the time when \sigma = 0.05.

   Expected:

   • 1. A Student law with T - q - 1 = 38 degrees of freedom, exactly, whatever the level and \sigma, because the projection removes the level and the ratio removes \sigma (2 pts).
   • 2. Qc = c - \bar c, so \|Qc\|^2 = 10 \times 30 / 40 = 7.5 and \|Qc\| = 2.74 (1.5 pts).
   • 3. The statistic has mean about \delta \|Qc\| / \sigma, so \delta_{\min} = 5.1 \times 0.05 / 2.74 = 0.093 (1.5 pts).

Grading rubric (proposal)

References

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

• Lecture notes of the course (listed on master-mva.com; no further detail displayed).
• A. Desolneux, L. Moisan, J.-M. Morel. From Gestalt Theory to Image Analysis: A Probabilistic Approach. Springer, Interdisciplinary Applied Mathematics 34, 2007 (the date given on master-mva.com; the printed volume bears 2008).
• G. Kanizsa. Organization in Vision: Essays on Gestalt Perception. Praeger Publishers, New York, 1979.
• A. Gordon, G. Glazko, X. Qiu, A. Yakovlev. Control of the mean number of false discoveries, Bonferroni and stability of multiple testing. The Annals of Applied Statistics, 1(1):179–190, 2007.
• “Various research papers”: the papers below that the chapters cite.

Textbooks

• E. L. Lehmann, J. P. Romano. Testing Statistical Hypotheses, 3rd edition. Springer, 2005.
• T. W. Anderson. An Introduction to Multivariate Statistical Analysis, 3rd edition. Wiley, 2003.
• S. M. Kay. Fundamentals of Statistical Signal Processing, Volume II: Detection Theory. Prentice Hall, 1998.
• D. G. Lowe. Perceptual Organization and Visual Recognition. Kluwer Academic Publishers, 1985.

Perception and the Helmholtz principle (ch. 1)

• M. Wertheimer. Untersuchungen zur Lehre von der Gestalt II. Psychologische Forschung, 4:301–350, 1923.
• F. Attneave. Some informational aspects of visual perception. Psychological Review, 61(3):183–193, 1954.

Multiple testing (ch. 2)

• Y. Benjamini, Y. Hochberg. Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society, Series B, 57(1):289–300, 1995.
• Y. Benjamini, D. Yekutieli. The control of the false discovery rate in multiple testing under dependency. The Annals of Statistics, 29(4):1165–1188, 2001.

A-contrario detection and geometric structures (ch. 3–4)

• A. Desolneux, L. Moisan, J.-M. Morel. Meaningful alignments. International Journal of Computer Vision, 40(1):7–23, 2000.
• A. Desolneux, L. Moisan, J.-M. Morel. Edge detection by Helmholtz principle. Journal of Mathematical Imaging and Vision, 14(3):271–284, 2001.
• R. Grompone von Gioi, J. Jakubowicz, J.-M. Morel, G. Randall. LSD: a fast line segment detector with a false detection control. IEEE Transactions on Pattern Analysis and Machine Intelligence, 32(4):722–732, 2010.
• R. Grompone von Gioi, J. Jakubowicz, J.-M. Morel, G. Randall. LSD: a line segment detector. Image Processing On Line, 2:35–55, 2012.
• B. Grosjean, L. Moisan. A-contrario detectability of spots in textured backgrounds. Journal of Mathematical Imaging and Vision, 33(3):313–337, 2009.
• L. Moisan, B. Stival. A probabilistic criterion to detect rigid point matches between two images and estimate the fundamental matrix. International Journal of Computer Vision, 57(3):201–218, 2004.
• J. Lezama, J.-M. Morel, G. Randall, R. Grompone von Gioi. A contrario 2D point alignment detection. IEEE Transactions on Pattern Analysis and Machine Intelligence, 37(3):499–512, 2015.
• V. Pătrăucean, P. Gurdjos, R. Grompone von Gioi. Joint a contrario ellipse and line detection. IEEE Transactions on Pattern Analysis and Machine Intelligence, 39(4):788–802, 2017.

Anomaly detection (ch. 5)

• I. S. Reed, X. Yu. Adaptive multiple-band CFAR detection of an optical pattern with unknown spectral distribution. IEEE Transactions on Acoustics, Speech, and Signal Processing, 38(10):1760–1770, 1990.
• A. Davy, T. Ehret, J.-M. Morel, M. Delbracio. Reducing anomaly detection in images to detection in noise. IEEE International Conference on Image Processing (ICIP), 2018.
• T. Ehret, A. Davy, J.-M. Morel, M. Delbracio. Image anomalies: a review and synthesis of detection methods. Journal of Mathematical Imaging and Vision, 61(5):710–743, 2019.

Time series and satellite imaging (ch. 6)

• E. S. Page. Continuous inspection schemes. Biometrika, 41(1–2):100–115, 1954.
• R. Killick, P. Fearnhead, I. A. Eckley. Optimal detection of changepoints with a linear computational cost. Journal of the American Statistical Association, 107(500):1590–1598, 2012.
• J. Verbesselt, R. Hyndman, G. Newnham, D. Culvenor. Detecting trend and seasonal changes in satellite image time series. Remote Sensing of Environment, 114(1):106–115, 2010.
• Z. Zhu, C. E. Woodcock. Object-based cloud and cloud shadow detection in Landsat imagery. Remote Sensing of Environment, 118:83–94, 2012.
• Z. Zhu, C. E. Woodcock. Continuous change detection and classification of land cover using all available Landsat data. Remote Sensing of Environment, 144:152–171, 2014.
• T. Dagobert, R. Grompone von Gioi, C. de Franchis, J.-M. Morel, C. Hessel. Cloud detection by luminance and inter-band parallax analysis for pushbroom satellite imagers. Image Processing On Line, 2020.

Image forensics (ch. 7)

• A. C. Popescu, H. Farid. Exposing digital forgeries in color filter array interpolated images. IEEE Transactions on Signal Processing, 53(10):3948–3959, 2005.
• T. Nikoukhah, J. Anger, T. Ehret, M. Colom, J.-M. Morel, R. Grompone von Gioi. JPEG grid detection based on the number of DCT zeros and its application to automatic and localized forgery detection. IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops (CVPRW), 2019.
• Q. Bammey, R. Grompone von Gioi, J.-M. Morel. An adaptive neural network for unsupervised mosaic consistency analysis in image forensics. IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), 2020.
• M. Gardella, P. Musé, J.-M. Morel, M. Colom. Noisesniffer: a fully automatic image forgery detector based on noise analysis. IEEE International Workshop on Biometrics and Forensics (IWBF), 2021.
• T. Ehret. Robust copy-move forgery detection by false alarms control. arXiv preprint, 2019.
• Q. Bammey, T. Nikoukhah, M. Gardella, R. Grompone von Gioi, M. Colom, J.-M. Morel. Non-semantic evaluation of image forensics tools: methodology and database. IEEE/CVF Winter Conference on Applications of Computer Vision (WACV), 2022.

Research directions (ch. 8)

• N. Myaskouvskey, Y. Gousseau, M. Lindenbaum. Beyond independence: an extension of the a contrario decision procedure. International Journal of Computer Vision, 101(1):22–44, 2013.
• K. J. Worsley, S. Marrett, P. Neelin, A. C. Vandal, K. J. Friston, A. C. Evans. A unified statistical approach for determining significant signals in images of cerebral activation. Human Brain Mapping, 4(1):58–73, 1996.
• S. Bates, E. Candès, L. Lei, Y. Romano, M. Sesia. Testing for outliers with conformal p-values. The Annals of Statistics, 51(1):149–178, 2023.
• R. Corvi, D. Cozzolino, G. Zingarini, G. Poggi, K. Nagano, L. Verdoliva. On the detection of synthetic images generated by diffusion models. IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP), 2023.
• F. Guillaro, D. Cozzolino, A. Sud, N. Dufour, L. Verdoliva. TruFor: leveraging all-round clues for trustworthy image forgery detection and localization. IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), 2023.
• K. Huang, Y. Wang, Z. Zhou, T. Ding, S. Gao, Y. Ma. Learning to parse wireframes in images of man-made environments. IEEE/CVF Conference on Computer Vision and Pattern Recognition (CVPR), 2018.
• R. Caye Daudt, B. Le Saux, A. Boulch. Fully convolutional siamese networks for change detection. IEEE International Conference on Image Processing (ICIP), 2018.