Artificial Intelligence and Computer Vision for Cultural Heritage

Contents

• Introduction and overview
• 0. Review: images, filters, Fourier analysis, least squares and probability
• 1. Technical and chemical imaging of paintings
• 2. Cross-modal image registration
• 3. Conditional expectation in multimodal imaging
• 4. Microstructure: craquelure and canvas
• 5. Text and image embeddings with pretrained networks
• 6. Image-to-image translation with CNNs
• 7. Complements: optics of paint, information and uncertainty
• 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 the basics of image processing and neural networks, and wants to apply them to paintings and drawings, where every image has a physical meaning.

Empirical first: photograph a painting in visible light and in the near infrared, and the two images disagree. In the infrared the paint becomes partly transparent and the artist’s underdrawing appears (ch. 1). Overlay the two images and nothing matches to the pixel, because each camera has its own lens distortion. After a sub-pixel registration (ch. 2), the pixels whose infrared value the colour does not explain point to underdrawing, retouching or a change of mind (ch. 3).

Art history and conservation science now collect gigabytes, sometimes terabytes, of images of a single painting, in several physical modalities. This course teaches the imaging physics behind those images and the computer vision that turns them into evidence: registration, density estimation, segmentation of cracks and canvas weave, pretrained embeddings and image-to-image networks. The course is given by Robert G. Erdmann (University of Amsterdam), as listed on master-mva.com.

Related MVA pages: 3D Computer Vision for camera geometry, Sub-Pixel Image Processing for precise shifts, Remote Sensing Data for multispectral imaging at another scale, and Deep Learning for Image Restoration and Synthesis for image-to-image networks.

What you will be able to do

• Explain what visible photography, infrared reflectography, reflectance spectroscopy, radiography, macro X-ray fluorescence and structured-light scanning measure, and at which depth of a painting.
• Calibrate colour by least squares, unmix reflectance spectra, and recover a height map from phase-shifted fringes.
• Register two modalities: detect, describe and match keypoints, estimate a homography with lens distortion robustly, and refine local offsets to a fraction of a pixel.
• Choose or learn a similarity measure that works when the intensities of two modalities are not related by a monotone map.
• Estimate a low-dimensional density, compute surprisal and conditional surprisal maps, and interpret them.
• Train a CNN that predicts colour from element maps, and use its residuals to locate retouching.
• Segment craquelure with scale-normalised Hessian filters and describe a crack network by its orientation statistics.
• Measure thread densities from a radiograph and test whether two canvases were cut from the same roll.
• Build a text-and-image search engine with CLIP, use semantic directions, and find patch correspondences with DINOv2 features.
• Explain the U-Net and neural style transfer, prove that the Gram-matrix loss is a maximum mean discrepancy, and design an image-to-image network for a conservation question.

Official objectives and syllabus (2026–2027)

The course description on master-mva.com starts from a change in the field. New and inexpensive imaging modalities, fast computation and fast networks let researchers collect gigabytes or terabytes of multimodal images of one artwork. These data can improve what we know about who made it, how it was made, its state of conservation and its relation to other artworks.

The description adds that heritage institutions are digitising millions of works and making them freely available. Researchers can then study central questions of art history, such as attribution, style and iconography, in new ways. Extracting this value requires sophisticated computational tools.

The aims are two-fold: to introduce the rich and little-known world of cultural heritage imaging and its problems, and to introduce a selection of the most important recent developments in computer vision and deep learning that solve them. The Université Paris-Saclay course detail lists the course as “AI and Computer vision for cultural heritage” and gives only its form of assessment.

Every official topic is taught on this page, except the instructor’s own data sets. The extensions beyond the syllabus are the review chapter 0 and the complements of chapter 7: the Kubelka–Munk optics of paint, mutual information, the information bound of contrastive learning and distribution-free uncertainty for reconstructions.

Where this course fits in the programme

Official prerequisites: Python with NumPy (required), and some familiarity with SciPy, OpenCV, scikit-image, pandas and a deep-learning framework (recommended). It is an Opening course of semester 2 (period 2), with 2.5 ECTS and 12 h of lectures.

• Preparation — the refreshers Scientific Python, NumPy and PyTorch (every code block), Linear Algebra, Analysis and Optimization (SVD, least squares) and Probability and Statistics (densities, conditional expectation).
• Semester 1 courses that prepare it — 3D Computer Vision (homographies, calibration), Object Recognition and Computer Vision (local features, CNNs), Deep Learning, Sub-Pixel Image Processing, Mathematical Models for Digital Photography (colour, lenses) and Introduction to Statistical Learning.
• Related courses of semester 2 — Remote Sensing Data (spectral imaging, registration), Point Clouds and 3D Modelling (scanning), Multimodal Explainable AI (CLIP), Representation Learning for Computer Vision and Medical Imaging (self-supervision), Deep Learning for Image Restoration and Synthesis, Generative Models for Images and Kernel Methods for Machine Learning (MMD).
• Afterwards — a heritage-imaging question makes a natural Research Internship.

Course outline

• Review: images, filters, Fourier analysis, least squares and probability — the M1 tools used later, in the form used later (self-study).
• Technical and chemical imaging of paintings — what each of six modalities measures, at which depth, with the physics that links images to materials.
• Cross-modal image registration — keypoints and learned descriptors, homographies and lens distortion, robust estimation, sub-pixel patchwise offsets.
• Conditional expectation in multimodal imaging — densities and surprisal of pixel vectors, and a CNN that predicts colour from chemistry.
• Microstructure: craquelure and canvas — crack segmentation and statistics, thread counting and canvas matching from radiographs.
• Text and image embeddings with pretrained networks — CLIP for search and semantic directions, DINOv2 for correspondences and PCA maps.
• Image-to-image translation with CNNs — the U-Net, Gram matrices and style transfer, and two case studies: ink removal and the Night Watch.
• Complements: optics of paint, information and uncertainty — Kubelka–Munk, mutual information, the InfoNCE bound, conformal intervals (self-study).
• Research directions and open problems — where the field is going and what is still unsolved (self-study).

Session plan

12 h in class = 4 sessions of 3 h. The volume is the Université Paris-Saclay programme table (12 h of lectures, 2.5 ECTS), and it matches the organisation announced on master-mva.com: four 3-hour lectures, each with a considerable part of hands-on work on real data. The split of the six official themes between the four sessions is a proposal of these notes.

Level 1 = foundations (the physics of the modalities, first models), 2 = core methods, 3 = advanced methods and case studies; the level never decreases.

Self-study, outside the 12 h (about 7 h in total): ch. 0 Review, to read before session 1 (≈ 2 h) · ch. 7 Complements: Kubelka–Munk, mutual information, InfoNCE, conformal intervals, after session 4 (≈ 3 h, with Ex. 17) · ch. 8 Research directions and open problems (≈ 2 h, with Ex. 18). The project comes on top (see Projects and assessment).

How to study this module

• Read ch. 0 before session 1, even quickly: chapters 2, 3 and 6 use its theorems without proving them again.
• Before each session, read the chapters of its row; after it, do the exercises of its row. ★ applies a definition, ★★ asks for a short proof or a computation, ★★★ is open-ended.
• Run every code block in a Jupyter notebook, the format of the official assessment. They use NumPy, SciPy, scikit-learn and PyTorch on a CPU, with fixed seeds, and take under 30 seconds each.
• The code works on synthetic stand-ins for museum data. For each one, write down which physical assumption the simulation makes, and how a real image could break it.
• Read the official papers in the order of the sessions: Key.Net, HardNet and AffNet, then MAGSAC++ (sessions 1–2); CLIP (session 3); DINOv2, U-Net and Gatys et al. (session 4).
• Keep Tummers (2011) for the end: it shows how connoisseurs weigh evidence, which is what the computational methods must feed.

0. Review: images, filters, Fourier analysis, least squares and probability

The course combines the physics of imaging, classical image processing and deep learning. This review collects the M1-level tools that later chapters use without proof: linear light, convolution, the Fourier transform, the singular value decomposition, conditional expectation and the vocabulary of convolutional networks.

0.1 Images and linear light

Physical models, such as the attenuation of light by a layer, act on quantities proportional to energy. Camera files are usually gamma-encoded: a stored sRGB value C \in [0, 1] is decoded to linear light by C_{\mathrm{lin}} = \begin{cases} C / 12.92, & C \le 0.04045, \\ \big((C + 0.055)/1.055\big)^{2.4}, & C > 0.04045. \end{cases} Every model of ch. 1 and ch. 3 assumes linear values.

0.2 Convolution, correlation and Gaussian smoothing

Three facts are used constantly.

• The Gaussian is separable: a 2-D smoothing is two 1-D smoothings, one per axis.
• Smoothings compose: G_\sigma * G_\tau = G_{\sqrt{\sigma^2 + \tau^2}}.
• Derivatives commute with convolution, so a smoothed derivative is a convolution with a derivative of the Gaussian: \partial_i (G_\sigma * I) = (\partial_i G_\sigma) * I .

0.3 The discrete Fourier transform

The shift theorem is the basis of phase correlation in ch. 2.5, and the convolution theorem makes large correlations cheap.

0.4 Least squares and the SVD

The ordinary problem \min_x \|Ax - b\|^2 with A of full column rank has the solution of the normal equations, x^\star = (A^\top A)^{-1} A^\top b . Homography estimation (ch. 2.3) uses the first form, colour calibration (ch. 1.1) the second.

0.5 Densities, information and conditional expectation

0.6 The toolbox of convolutional networks

• A convolutional layer applies learned k \times k filters to all input channels, adds a bias and a non-linearity such as ReLU, \max(0, u). Pooling or a stride s > 1 reduces the resolution.
• The receptive field of a unit is the input region it depends on. For layers with kernel sizes k_l and strides s_l, it is r_L = 1 + \sum_{l=1}^{L} (k_l - 1) \prod_{j < l} s_j .
• A softmax with temperature \tau turns scores z_j into probabilities e^{z_j/\tau} / \sum_k e^{z_k/\tau}. The cosine similarity of two vectors is their inner product after normalisation to unit length.
• Training minimises an average loss by stochastic gradient methods such as Adam; automatic differentiation computes the gradients.

The code below checks Theorem 0.4 and Proposition 0.5 numerically.

import numpy as np

rng = np.random.default_rng(0)
n = 64
f, g = rng.standard_normal(n), rng.standard_normal(n)

# Convolution theorem (circular, 1-D): DFT(f * g) = DFT(f) DFT(g)
circ = np.array([sum(f[m] * g[(k - m) % n] for m in range(n)) for k in range(n)])
err = np.abs(np.fft.fft(circ) - np.fft.fft(f) * np.fft.fft(g)).max()
print(f"convolution theorem, max error: {err:.1e}")

# Shift theorem: a circular shift by t multiplies the DFT by exp(-2 i pi k t / n)
t, k = 5, np.arange(n)
lhs = np.fft.fft(np.roll(f, t))
rhs = np.exp(-2j * np.pi * k * t / n) * np.fft.fft(f)
print(f"shift theorem, max error: {np.abs(lhs - rhs).max():.1e}")

# Homogeneous least squares: argmin ||A x|| over unit x = last right singular vector
A = rng.standard_normal((20, 9))
x_true = rng.standard_normal(9)
x_true /= np.linalg.norm(x_true)
A -= np.outer(A @ x_true, x_true)          # x_true is now an exact null vector of A
_, s, Vt = np.linalg.svd(A)
print(f"smallest singular value: {s[-1]:.1e}, |<v_9, x_true>| = {abs(Vt[-1] @ x_true):.6f}")

1. Technical and chemical imaging of paintings

A painting is a layered object, and each imaging modality measures one physical quantity at one depth.

1.1 Visible photography

A colour chart with n patches of known linear values y_k \in \mathbb{R}^3 is photographed with the painting, giving camera values x_k \in \mathbb{R}^3.

1.2 Infrared reflectography

Paint scatters less, and many pigments absorb less, at longer wavelengths, so infrared light crosses layers that are opaque in the visible. Carbon-based drawing materials absorb at all wavelengths.

The contrast grows with T(\lambda), hence with the wavelength.

1.3 Reflectance imaging spectroscopy

A RIS camera records a cube r(x, \lambda) of often hundreds of bands, which identify pigments far better than three colour channels.

It holds for side-by-side areas of pigment, not for pigments mixed in one layer.

import numpy as np
from scipy.optimize import nnls

# Moseley's law for K-alpha lines: E = (3/4) * 13.6 eV * (Z - 1)^2
for name, Z, tabulated in [("Ca", 20, 3.69), ("Fe", 26, 6.40), ("Cu", 29, 8.05)]:
    print(f"{name}: Moseley {0.75 * 13.6e-3 * (Z - 1) ** 2:5.2f} keV, tabulated {tabulated:5.2f} keV")

# Linear unmixing of a synthetic reflectance spectrum (400-1000 nm, 121 bands)
lam = np.linspace(400, 1000, 121)
sig = lambda x, c, w: 1 / (1 + np.exp(-(x - c) / w))
E = np.stack([0.85 + 0.05 * sig(lam, 450, 30),                      # synthetic "white"
              0.05 + 0.80 * sig(lam, 600, 12),                      # synthetic "red"
              0.10 + 0.30 * np.exp(-((lam - 460) / 40) ** 2)        # synthetic "blue"
              + 0.50 * sig(lam, 800, 40)], axis=1)                  # E has shape (121, 3)
a_true = np.array([0.2, 0.5, 0.3])
rng = np.random.default_rng(1)
r = E @ a_true + 0.01 * rng.standard_normal(lam.size)

# Non-negative least squares with a heavily weighted sum-to-one row
w = 100.0
a_hat, _ = nnls(np.vstack([E, w * np.ones(3)]), np.append(r, w))
print("abundances: true", a_true, " estimated", a_hat.round(3))

# Spectral angle between the measured spectrum and each endmember (radians)
sam = np.arccos(E.T @ r / (np.linalg.norm(E, axis=0) * np.linalg.norm(r)))
print("spectral angles:", sam.round(3))

1.4 Radiography

Lead white attenuates X-rays far more than binders, canvas or wood, so a radiograph of an old master is mostly a map of lead-white thickness. Where the ground fills the gaps between the threads, the canvas weave appears (ch. 4.2).

1.5 Macro X-ray fluorescence

An X-ray beam scanned over the painting makes each element emit lines at energies specific to it. One spectrum per position gives one map per element.

Deep layers are seen through the layers above, which absorb low energies most, so element maps also carry depth information. Dik et al. (2008) used them to reveal a head hidden under a painting by Van Gogh.

1.6 Structured-light scanning

A projector casts sinusoidal fringes on the surface and a camera views them from another direction. Relief displaces the fringes.

import numpy as np

# Three-step phase shifting on a synthetic impasto bump (units: mm)
x, y = np.meshgrid(np.linspace(-10, 10, 201), np.linspace(-10, 10, 201))
h = 0.5 * np.exp(-(x**2 + y**2) / 8)                 # height map, 0.5 mm peak
p, theta = 2.0, np.deg2rad(30)                        # fringe period, projection angle
phi0 = 2 * np.pi * x / p                              # phase on the flat reference plane
phi = phi0 + 2 * np.pi * h * np.tan(theta) / p        # phase displaced by the relief
A, B = 0.5, 0.4
rng = np.random.default_rng(0)
I1, I2, I3 = (A + B * np.cos(phi + d) + 0.005 * rng.standard_normal(phi.shape)
              for d in (0, 2 * np.pi / 3, 4 * np.pi / 3))

wrapped = np.arctan2(np.sqrt(3) * (I3 - I2), 2 * I1 - I2 - I3)
dphi = np.angle(np.exp(1j * (wrapped - phi0)))        # phase difference, wrapped to (-pi, pi]
h_hat = dphi * p / (2 * np.pi * np.tan(theta))
print(f"max phase change {dphi.max():.3f} rad (below pi: no unwrapping needed)")
print(f"height RMS error: {1000 * np.sqrt(np.mean((h_hat - h) ** 2)):.1f} micrometres")

2. Cross-modal image registration

Registration maps each pixel of one image to the corresponding point of another. Across modalities intensities differ, so features must be designed or learned to agree. The pipeline has three stages: keypoint matches, a robust global transform, local sub-pixel offsets.

2.1 Keypoints: scale space and the structure tensor

Two small eigenvalues mean a flat region, one large eigenvalue an edge, two large ones a corner. The Harris detector keeps local maxima of \det J - \kappa\,(\operatorname{tr} J)^2, with \kappa \approx 0.05.

Key.Net (Barroso-Laguna et al., 2019) combines such handcrafted derivative filters with learned filters in a multi-scale pyramid, trained on pairs related by random synthetic transforms.

2.2 Descriptors and matching

AffNet (Mishkin et al., 2018) estimates the affine shape of each keypoint region, trained to make descriptors discriminative rather than merely repeatable. HardNet (Mishchuk et al., 2017) maps a 32 \times 32 patch to a unit vector in dimension 128.

A nearest-neighbour match is kept if it is mutual and if the ratio of the first to the second distance is below about 0.8 (Lowe, 2004). Descriptors trained on photographs transfer well between VIS and IRR, which share many edges, less well to element maps.

2.3 Homographies and lens distortion

A painting is nearly planar, so its images are related by homographies. Each match gives two linear equations in the entries of H; with four or more, the direct linear transform (DLT) solves A h = 0 by Proposition 0.5, after normalising the coordinates (Hartley, 1997).

2.4 Robust estimation: RANSAC and MAGSAC++

Many matches between modalities are wrong. RANSAC (Fischler and Bolles, 1981) fits the model to random minimal samples and keeps the one with the most inliers below a threshold.

For w = 0.4, s = 4 and p = 0.99: 178 draws. The fixed threshold is the weak point; MAGSAC++ (Barath et al., 2020) marginalises over the noise scale \sigma instead, and refits by iteratively reweighted least squares.

The weight decreases smoothly, without a threshold; MAGSAC++ applies this construction to each model’s residual distribution, with incomplete gamma functions.

import numpy as np
from scipy.special import erfc

def normalise(x):
    """Similarity T sending the points to zero mean and mean distance sqrt(2)."""
    c = x.mean(0)
    s = np.sqrt(2) / np.linalg.norm(x - c, axis=1).mean()
    return np.array([[s, 0, -s * c[0]], [0, s, -s * c[1]], [0, 0, 1]])

def dlt(x, y, w=None):
    """Weighted normalised DLT: homography H with y ~ H x."""
    w = np.ones(len(x)) if w is None else w
    Tx, Ty = normalise(x), normalise(y)
    xh = np.c_[x, np.ones(len(x))] @ Tx.T
    yh = np.c_[y, np.ones(len(y))] @ Ty.T
    z = np.zeros_like(xh)
    A = np.r_[np.c_[z, -xh, yh[:, 1:2] * xh], np.c_[xh, z, -yh[:, 0:1] * xh]]
    A *= np.sqrt(np.r_[w, w])[:, None]
    H = np.linalg.svd(A)[2][-1].reshape(3, 3)          # Proposition 0.5
    H = np.linalg.inv(Ty) @ H @ Tx
    return H / H[2, 2]

def apply(H, x):
    yh = np.c_[x, np.ones(len(x))] @ H.T
    return yh[:, :2] / yh[:, 2:]

rng = np.random.default_rng(0)
H_true = np.array([[1.02, 0.03, 15.0], [-0.02, 0.98, -8.0], [2e-5, -1e-5, 1.0]])
n, w_in = 300, 0.4
x = rng.uniform(0, 1000, (n, 2))
y = apply(H_true, x) + 0.5 * rng.standard_normal((n, 2))          # inliers: 0.5 px noise
out = rng.random(n) > w_in
y[out] = rng.uniform(0, 1000, (out.sum(), 2))                     # 60 % gross outliers

# RANSAC with the number of draws of Proposition 2.6
N = int(np.ceil(np.log(1 - 0.99) / np.log(1 - w_in**4)))
best, thr = None, 3.0
for _ in range(N):
    idx = rng.choice(n, 4, replace=False)
    H = dlt(x[idx], y[idx])
    inl = np.linalg.norm(apply(H, x) - y, axis=1) < thr
    if best is None or inl.sum() > best.sum():
        best = inl
H_ransac = dlt(x[best], y[best])

# Iteratively reweighted DLT with the sigma-marginalised weight of Proposition 2.7
H_irls = H_ransac
for _ in range(10):
    r = np.linalg.norm(apply(H_irls, x) - y, axis=1)
    H_irls = dlt(x, y, erfc(r / (np.sqrt(2) * thr)))

grid = np.stack(np.meshgrid(np.linspace(0, 1000, 11), np.linspace(0, 1000, 11)), -1).reshape(-1, 2)
for name, H in [("all points", dlt(x, y)), ("RANSAC + refit", H_ransac), ("sigma-weighted IRLS", H_irls)]:
    err = np.linalg.norm(apply(H, grid) - apply(H_true, grid), axis=1)
    print(f"{name:20s}: mean error {err.mean():8.3f} px, max {err.max():8.3f} px")
print(f"RANSAC draws: {N}, inliers kept: {best.sum()} of {(~out).sum()} true inliers")

2.5 Patchwise offsets to sub-pixel precision

After the global transform, local offsets of a few pixels remain (scanner drift, deformed supports, XRF rasters). An offset is estimated per tile, then interpolated into a smooth field.

On a single straight edge \lambda_2 = 0: only the offset across the edge is measured. Tiles are weighted by their smallest eigenvalue.

NCC tolerates gain, offset and sign changes, not a non-monotone relation between modalities. The remedy is a non-linear filter per modality, such as a gradient magnitude, mapping both images to a common representation.

A parabola through the values c_{-1}, c_0, c_{+1} around the peak refines it on each axis: \delta = \frac{c_{-1} - c_{+1}}{2\,(c_{-1} - 2 c_0 + c_{+1})} . Stacking all tiles in one tensor lets a single batched FFT process them in parallel.

import numpy as np
from scipy import ndimage as ndi

rng = np.random.default_rng(0)
n, T = 512, 64                                         # image size, tile size
A = ndi.gaussian_filter(rng.standard_normal((n, n)), 2.0)
A = (A - A.min()) / (A.max() - A.min())                # "modality 1" in [0, 1]

# Smooth displacement field of up to about 3 px, then a non-monotone change of contrast
gy, gx = np.mgrid[0:n, 0:n]
dy = 2.5 * np.sin(2 * np.pi * gx / n) + 0.3
dx = 2.0 * np.cos(2 * np.pi * gy / n) - 0.4
B = np.abs(ndi.map_coordinates(A, [gy - dy, gx - dx], order=3) - 0.5)   # "modality 2"

def tiles(img):
    """Stack of non-overlapping T x T tiles, shape (number of tiles, T, T)."""
    return img.reshape(n // T, T, n // T, T).swapaxes(1, 2).reshape(-1, T, T)

def phase_correlation(a, b):
    """Batched phase correlation (Proposition 2.10) with parabolic sub-pixel refinement."""
    win = np.outer(np.hanning(T), np.hanning(T))
    Fa = np.fft.fft2((a - a.mean((1, 2), keepdims=True)) * win)
    Fb = np.fft.fft2((b - b.mean((1, 2), keepdims=True)) * win)
    R = Fb * np.conj(Fa)
    c = np.real(np.fft.ifft2(R / (np.abs(R) + 1e-12)))   # one correlation surface per tile
    i, j = np.unravel_index(c.reshape(len(c), -1).argmax(1), (T, T))
    k = np.arange(len(c))
    refine = lambda cm, c0, cp: 0.5 * (cm - cp) / (cm - 2 * c0 + cp)
    di = refine(c[k, (i - 1) % T, j], c[k, i, j], c[k, (i + 1) % T, j])
    dj = refine(c[k, i, (j - 1) % T], c[k, i, j], c[k, i, (j + 1) % T])
    wrap = lambda q: (q + T // 2) % T - T // 2
    return np.stack([wrap(i) + di, wrap(j) + dj], 1)

centres = np.arange(n // T) * T + T // 2
true = np.stack([dy[np.ix_(centres, centres)].ravel(), dx[np.ix_(centres, centres)].ravel()], 1)
grad = lambda img: ndi.gaussian_gradient_magnitude(img, 1.0)   # a non-linear common representation
for name, a, b in [("raw intensities", A, B), ("gradient magnitude", grad(A), grad(B))]:
    err = np.linalg.norm(phase_correlation(tiles(a), tiles(b)) - true, axis=1)
    print(f"{name:19s}: median error {np.median(err):6.3f} px, tiles within 0.5 px: {np.mean(err < 0.5):.0%}")

Here the second modality is |A - 0.5|. Raw correlation fails, while the gradient magnitude, unchanged by this map except on the level line A = 0.5, registers every tile.

3. Conditional expectation in multimodal imaging

Once registered, every pixel is a short vector of measurements, such as (R, G, B, IRR). The question “what is unusual here?” becomes a question about densities and conditional expectations.

3.1 Density estimation and surprisal

With d = 4 and 10^8 pixels, the rate n^{-1/2} is excellent: densities of a few channels are easy to estimate, which is why the syllabus speaks of low-dimensional data.

For densities these values can be negative; only comparisons matter. A high s(u, v) flags a rare pixel, perhaps only a rare colour. A high s(v \mid u) flags an infrared value that the colour does not explain, such as an underdrawing or a change of composition. The syllabus illustrates surprisal maps on paintings by Hieronymus Bosch.

import numpy as np
from sklearn.neighbors import KernelDensity
from sklearn.metrics import roc_auc_score

rng = np.random.default_rng(0)
# Synthetic pixels (R, G, B, IRR): three common paints and one rare red paint
means = np.array([[0.85, 0.70, 0.60, 0.80],     # light flesh tone, bright in infrared
                  [0.30, 0.35, 0.50, 0.55],     # blue-grey
                  [0.20, 0.25, 0.15, 0.30],     # dark brown
                  [0.75, 0.15, 0.10, 0.70]])    # rare red
counts = [1500, 1000, 800, 40]
X = np.vstack([m + 0.03 * rng.standard_normal((c, 4)) for m, c in zip(means, counts)])
kind = np.repeat(np.arange(4), counts)
under = (kind == 0) & (rng.random(len(X)) < 0.05)   # carbon underdrawing below the flesh tone
X[under, 3] -= 0.25                                 # darker in IRR only, same colour

kde4 = KernelDensity(bandwidth=0.04).fit(X)         # p(R, G, B, IRR)
kde3 = KernelDensity(bandwidth=0.04).fit(X[:, :3])  # p(R, G, B)
joint = -kde4.score_samples(X)                      # surprisal, in nats
cond = joint + kde3.score_samples(X[:, :3])         # -log p(IRR | RGB)

k = under.sum()                                     # flag as many pixels as there are true ones
for name, s in [("joint", joint), ("conditional", cond)]:
    top = np.argsort(s)[-k:]
    print(f"{name:11s} surprisal: AUC {roc_auc_score(under, s):.3f}, among the {k} most surprising "
          f"pixels {under[top].sum()} are underdrawing and {(kind[top] == 3).sum()} are the rare red")

# Nadaraya-Watson estimate of E[IRR | RGB] and a standardised residual
def nw(x_query, x, y, h=0.04):
    d2 = ((x_query[:, None, :] - x[None, :, :]) ** 2).sum(-1)
    w = np.exp(-d2 / (2 * h * h))
    return (w * y).sum(1) / w.sum(1)

res = X[:, 3] - nw(X[:, :3], X[:, :3], X[:, 3])
z = res / res[~under].std()
print(f"residual z-score: underdrawing {z[under].mean():.1f}, other pixels {z[~under].mean():.2f}")

Joint surprisal mixes the rare red with the underdrawing; conditional surprisal isolates the underdrawing.

3.2 Predicting colour from chemistry

A neighbourhood is needed rather than a single pixel. XRF maps are photon-limited and noisy, light diffuses inside the paint, and lines from deep layers are attenuated. Retouching made with a different colour over the same chemistry, or with materials absent from the input maps, leaves a large residual. The syllabus applies this to the detection of restorations in paintings by Rembrandt.

Train and test on separate blocks of the image: neighbouring pixels are correlated, and a random split of pixels overstates accuracy.

import numpy as np
import torch
from scipy import ndimage as ndi
from sklearn.metrics import roc_auc_score

torch.set_num_threads(1)
torch.manual_seed(0)
rng = np.random.default_rng(0)
S = 48
# Absorption of light (R, G, B) per unit concentration: Hg, Cu, Fe; Pb scatters and brightens
alpha = np.array([[0.1, 2.0, 2.0],      # Hg: vermilion-like red
                  [2.0, 0.6, 0.4],      # Cu: blue-green
                  [0.3, 0.8, 1.6]])     # Fe: earth, yellow-brown

def painting():
    """Element concentrations (Pb, Hg, Cu, Fe) and the visible colour they produce."""
    f = ndi.gaussian_filter(rng.standard_normal((4, S, S)), (0, 4, 4))
    c = np.log1p(np.exp(2 * f / f.std() - 0.5))              # smooth, non-negative maps
    rgb = np.exp(-np.einsum("ehw,ec->chw", c[1:], alpha)) * (0.4 + 0.6 * np.tanh(c[0]))[None]
    rgb = ndi.gaussian_filter(rgb, (0, 1.0, 1.0))           # light diffusion in the paint
    xrf = rng.poisson(20 * c) / 20.0                          # photon-counting noise of MA-XRF
    return xrf.astype(np.float32), rgb.astype(np.float32)

train = [painting() for _ in range(6)]
xrf_t, rgb_t = painting()                                     # a separate test painting
rgb_t[:, 15:26, 27:39] += np.array([0.06, -0.03, -0.04])[:, None, None]   # retouching, chemistry unchanged
patch = np.zeros((S, S), bool); patch[15:26, 27:39] = True

def net(k):
    p = k // 2
    return torch.nn.Sequential(
        torch.nn.Conv2d(4, 24, k, padding=p), torch.nn.ReLU(),
        torch.nn.Conv2d(24, 24, k, padding=p), torch.nn.ReLU(),
        torch.nn.Conv2d(24, 3, 1), torch.nn.Sigmoid())

X = torch.tensor(np.stack([a for a, _ in train])); Y = torch.tensor(np.stack([b for _, b in train]))
for name, k in [("pixel-wise (1x1)", 1), ("CNN (5x5)", 5)]:
    model = net(k)
    opt = torch.optim.Adam(model.parameters(), lr=3e-3)
    for step in range(200):
        opt.zero_grad()
        loss = torch.mean((model(X) - Y) ** 2)            # estimates E[VIS | XRF neighbourhood]
        loss.backward(); opt.step()
    with torch.no_grad():
        pred = model(torch.tensor(xrf_t)[None])[0].numpy()
    res = np.linalg.norm(rgb_t - pred, axis=0)
    print(f"{name:16s}: RMSE outside the retouching {np.sqrt(np.mean(res[~patch] ** 2)):.4f}, "
          f"residual AUC for the retouching {roc_auc_score(patch.ravel(), res.ravel()):.3f}")

4. Microstructure: craquelure and canvas

Some evidence lies in the microstructure of a painting: its cracks, and the weave of its canvas seen in radiographs.

4.1 Craquelure

Cracks are thin dark curves. At scale \sigma, the Hessian of L_\sigma has a large positive eigenvalue across a dark line and a small one along it.

So with \gamma = 3/4 (Lindeberg, 1998), the selected scale equals the crack width. A multi-scale line measure, as in the vesselness filter of Frangi et al. (1998), keeps pixels where the large eigenvalue dominates the small one.

Coherence near 1 means one dominant direction, as on many panels. With crack density, island areas and junction degrees, it feeds a classifier. The syllabus takes its examples from Van Eyck; Cornelis et al. (2013) detect cracks on the Ghent Altarpiece.

import numpy as np
from scipy import ndimage as ndi

rng = np.random.default_rng(0)
n, w = 256, 1.5                                     # image size, crack profile width (px)

def cracks(stretch):
    """Voronoi crack network; stretch > 1 elongates the islands along the vertical axis."""
    seeds = rng.uniform(0, n, (60, 2))
    yy, xx = np.mgrid[0:n, 0:n]
    d = ((yy[..., None] - seeds[:, 0]) / stretch) ** 2 + (xx[..., None] - seeds[:, 1]) ** 2
    lab = d.argmin(-1)
    edge = (np.diff(lab, axis=0, append=lab[-1:]) != 0) | (np.diff(lab, axis=1, append=lab[:, -1:]) != 0)
    img = 0.8 + 0.03 * rng.standard_normal((n, n))                  # paint surface
    img -= 0.5 * np.minimum(1, ndi.gaussian_filter(edge * 1.0, w) * np.sqrt(2 * np.pi) * w)
    return img, edge

def ridge(img, sigmas, gamma=0.75):
    """Scale-normalised Hessian line measure for dark lines; best scale and local orientation."""
    best, theta = np.zeros_like(img), np.zeros_like(img)
    for s in sigmas:
        Hyy = ndi.gaussian_filter(img, s, order=(2, 0)); Hxx = ndi.gaussian_filter(img, s, order=(0, 2))
        Hxy = ndi.gaussian_filter(img, s, order=(1, 1))
        root = np.sqrt(((Hxx - Hyy) / 2) ** 2 + Hxy ** 2)
        lmax, lmin = (Hxx + Hyy) / 2 + root, (Hxx + Hyy) / 2 - root
        r = s ** (2 * gamma) * np.maximum(lmax, 0) * np.exp(-(lmin / (lmax + 1e-9)) ** 2 / 0.5)
        th = 0.5 * np.arctan2(2 * Hxy, Hxx - Hyy) + np.pi / 2    # direction along the line
        upd = r > best
        best[upd], theta[upd] = r[upd], th[upd]
    return best, theta

for name, stretch in [("isotropic network (canvas-like)", 1.0), ("elongated islands (panel-like)", 4.0)]:
    img, edge = cracks(stretch)
    R, th = ridge(img, sigmas=[1.0, 1.5, 2.0, 3.0])
    mask = R > 0.5 * np.percentile(R, 99)
    truth = ndi.binary_dilation(edge, iterations=2)
    precision, recall = (mask & truth).sum() / mask.sum(), (mask & edge).sum() / edge.sum()
    Rbar = np.abs(np.mean(np.exp(2j * th[mask])))               # Definition 4.3
    print(f"{name:32s}: precision {precision:.2f}, recall {recall:.2f}, "
          f"crack density {mask.mean():.3f}, orientation coherence {Rbar:.2f}")

# Scale selection on a 1-D Gaussian line of width w: argmax of sigma^(2 gamma) L_xx(0)
sig = np.linspace(0.3, 6, 2000)
for gamma in (0.75, 1.0):
    resp = sig ** (2 * gamma) * w * (w**2 + sig**2) ** -1.5
    print(f"gamma = {gamma}: selected sigma / w = {sig[resp.argmax()] / w:.3f}")

4.2 Canvas analysis from radiographs

In a radiograph, the lead-white ground that fills the gaps between threads draws the weave (Proposition 1.6). The threads running along a roll of canvas are the warp; the weft threads cross them.

Scanning tiles along the painting gives a warp-density profile; zero-padding and a parabolic fit locate the peak to a fraction of a bin.

Matching therefore maximises the normalised correlation of two warp-density profiles over shifts and flips. The syllabus applies this to works by Poussin, Van Gogh, Rembrandt and Frans Hals; the method is described by Johnson, Johnson and Erdmann (2013).

import numpy as np
from scipy import ndimage as ndi

rng = np.random.default_rng(0)
px = 0.01                                          # pixel size in cm (0.1 mm)

def roll(width_cm=60, density=14.0):
    """Warp thread positions (cm) across the width of a roll: irregular, correlated spacings."""
    k = int(width_cm * density)
    eps = 0.08 * ndi.gaussian_filter1d(rng.standard_normal(k), 8) * np.sqrt(8) + 0.02 * rng.standard_normal(k)
    return np.cumsum((1 + eps) / density)

def radiograph(warp, x0_cm, w_px=1536, h_px=256, weft_density=12.0):
    """X-ray image of a piece cut at x0 across the roll: the ground fills the interstices."""
    x = x0_cm + px * np.arange(w_px)
    u = np.interp(x, warp, np.arange(len(warp)))                      # warp thread phase
    v = np.cumsum((1 + 0.05 * rng.standard_normal(h_px)) * weft_density * px)  # weft phase
    img = 0.2 * np.cos(2 * np.pi * u)[None, :] + 0.2 * np.cos(2 * np.pi * v)[:, None]
    return img + 0.15 * rng.standard_normal((h_px, w_px))

def warp_profile(img, T=64, step=32, pad=512):
    """Warp thread density (threads/cm) along x, from the peak of windowed tile spectra."""
    win = np.outer(np.hanning(T), np.hanning(T))
    f = np.fft.rfftfreq(pad) / px                                      # threads per cm
    band = (f > 8) & (f < 20)
    prof = []
    for x in range(0, img.shape[1] - T + 1, step):
        spec = [np.abs(np.fft.rfft2(img[y:y + T, x:x + T] * win, s=(pad, pad)))[0]
                for y in range(0, img.shape[0] - T + 1, T)]
        s = np.mean(spec, 0)
        i = np.flatnonzero(band)[s[band].argmax()]
        d = 0.5 * (s[i - 1] - s[i + 1]) / (s[i - 1] - 2 * s[i] + s[i + 1])   # parabolic refinement
        prof.append((i + d) * (f[1] - f[0]))
    return np.array(prof)

def best_match(a, b, min_overlap=24):
    """Largest normalised correlation of two density profiles over shifts (in tile steps)."""
    out = []
    for lag in range(-len(b) + min_overlap, len(a) - min_overlap + 1):
        p, q = a[max(lag, 0):], b[max(-lag, 0):]
        m = min(len(p), len(q))
        out.append((np.corrcoef(p[:m], q[:m])[0, 1], lag))
    return max(out)

R1 = roll()
A = warp_profile(radiograph(R1, x0_cm=10.0))
B = warp_profile(radiograph(R1, x0_cm=14.0))       # same roll, overlapping 11 cm of warp threads
print(f"mean warp density of A: {A.mean():.2f} threads/cm (the roll was woven at 14)")
same, lag = best_match(A, B)
null = np.array([best_match(A, warp_profile(radiograph(roll(), x0_cm=10.0)))[0] for _ in range(20)])
print(f"same roll: correlation {same:.2f} at a shift of {lag * 32 * px:.2f} cm (true: 4 cm)")
print(f"20 unrelated rolls: mean {null.mean():.2f}, max {null.max():.2f}")
print(f"z-score of the match against the unrelated rolls: {(same - null.mean()) / null.std():.1f}")

Unrelated canvases reach 0.8 by chance: a match needs a calibrated null.

5. Text and image embeddings with pretrained networks

Pretrained networks map images and texts to vectors whose geometry reflects meaning, turning search and comparison into nearest-neighbour problems.

5.1 CLIP: a joint embedding of images and texts

CLIP (Radford et al., 2021) trained such encoders on 400 million image–text pairs from the internet. A label set becomes a zero-shot classifier: embed prompts such as “a painting of a dog”, and pick the closest.

5.2 Search engine and semantic directions

A search engine embeds the collection once, then returns the images of largest cosine similarity to a text or image query. At millions of images, an approximate nearest-neighbour index such as FAISS (Johnson et al., 2021) replaces the exhaustive scan.

Directions such as “night scene” minus “day scene” rank a collection along a concept without any labelled example.

5.3 High-resolution microstructural comparison

Most CLIP models see 224 \times 224 pixels; a painting has billions. Embed tiles cut at native resolution, with unit embeddings a_1, \dots, a_m for one work and b_1, \dots, b_n for the other. Compare the two sets by the symmetric mean best-match similarity \mathrm{sim}(A, B) = \frac{1}{2m} \sum_{i=1}^{m} \max_{j} \langle a_i, b_j \rangle + \frac{1}{2n} \sum_{j=1}^{n} \max_{i} \langle a_i, b_j \rangle . Brushwork, craquelure or paper texture then become comparable. Micrographs are rare in CLIP’s training data, so trust relative, not absolute, similarities.

The code trains a miniature CLIP from scratch on the scikit-learn digits, with bag-of-words captions such as “a sketch of a seven on a light ground”.

import numpy as np
import torch
import torch.nn.functional as F
from sklearn.datasets import load_digits
from sklearn.metrics import roc_auc_score

torch.set_num_threads(1)
torch.manual_seed(0)
rng = np.random.default_rng(0)
digits = load_digits()
X = digits.data / 16.0
y = digits.target
inv = rng.random(len(X)) < 0.5                    # half the images: light ground, dark ink
X[inv] = 1.0 - X[inv]
names = "zero one two three four five six seven eight nine".split()
templates = ["a drawing of the digit {}", "the number {} written by hand", "a sketch of a {}"]
ground = {False: "on a dark ground", True: "on a light ground"}

def caption(i):
    return templates[rng.integers(3)].format(names[y[i]]) + " " + ground[bool(inv[i])]

vocab = {w: k for k, w in enumerate(sorted({w for t in templates for w in t.split()}
                                           | set(names) | {"on", "a", "dark", "light", "ground"}))}
def tokens(texts):
    """Bag of words, as a normalised count matrix (a stand-in for a text transformer)."""
    M = torch.zeros(len(texts), len(vocab))
    for r, t in enumerate(texts):
        for w in t.split():
            M[r, vocab[w]] += 1
    return M / M.sum(1, keepdim=True)

img_enc = torch.nn.Sequential(torch.nn.Linear(64, 128), torch.nn.ReLU(), torch.nn.Linear(128, 32))
txt_enc = torch.nn.Sequential(torch.nn.Linear(len(vocab), 128), torch.nn.ReLU(), torch.nn.Linear(128, 32))
log_scale = torch.nn.Parameter(torch.tensor(np.log(1 / 0.07), dtype=torch.float32))
params = list(img_enc.parameters()) + list(txt_enc.parameters()) + [log_scale]
opt = torch.optim.Adam(params, lr=2e-3)
train, test = np.arange(1300), np.arange(1300, len(X))
Xt = torch.tensor(X, dtype=torch.float32)

for step in range(400):
    b = rng.choice(train, 256, replace=False)
    zi = F.normalize(img_enc(Xt[b]), dim=1)
    zt = F.normalize(txt_enc(tokens([caption(i) for i in b])), dim=1)
    logits = log_scale.exp().clamp(max=100) * zi @ zt.T          # s_ij of Definition 5.1
    target = torch.arange(len(b))
    loss = (F.cross_entropy(logits, target) + F.cross_entropy(logits.T, target)) / 2
    opt.zero_grad(); loss.backward(); opt.step()

with torch.no_grad():
    zi = F.normalize(img_enc(Xt[test]), dim=1)
    enc = lambda texts: F.normalize(txt_enc(tokens(texts)), dim=1)
    # Zero-shot classification with a prompt ensemble over the two grounds
    zc = F.normalize(sum(enc([f"a drawing of the digit {d} {g}" for d in names]) for g in ground.values()), dim=1)
    acc = ((zi @ zc.T).argmax(1).numpy() == y[test]).mean()
    # Text-to-image search: the 10 nearest test images to one query
    q = enc(["a sketch of a seven on a light ground"])
    top = (zi @ q.T).squeeze(1).argsort(descending=True)[:10].numpy()
    hits = np.mean((y[test][top] == 7) & inv[test][top])
    # Semantic direction (Definition 5.3): light ground minus dark ground
    d = enc(["a drawing of the digit five on a light ground"]) - enc(["a drawing of the digit five on a dark ground"])
    score = (zi @ d.T).squeeze(1).numpy()
print(f"final loss {loss.item():.3f}, temperature {1 / log_scale.exp().item():.3f}")
print(f"zero-shot accuracy on {len(test)} held-out images: {acc:.3f}")
print(f"query 'a sketch of a seven on a light ground': precision at 10 = {hits:.1f}")
print(f"semantic direction 'light ground': AUC {roc_auc_score(inv[test], score):.3f}")

5.4 DINOv2: self-supervised patch embeddings

DINOv2 (Oquab et al., 2024) trains vision transformers without labels on 142 million curated images. A student network matches a teacher, its moving average, across crops (DINO loss) and predicts masked patches (iBOT loss). The models cut images into patches of 14 \times 14 pixels and return one embedding per patch.

Because the embeddings encode parts rather than pixels, correspondences link the same figure in a painting, a copy and a print after it, despite changes of medium.

5.5 Visualisation by PCA

Fit a PCA to the patch embeddings of one or several images, and display the first three components as red, green and blue. Parts with similar meaning receive similar colours across images, and the first component often separates the subject from the background.

The code uses a contrast-invariant stand-in for the patch embeddings, on a scikit-learn sample image and a rotated, tone-inverted copy of it; DINOv2 patch tokens plug into the same pipeline.

import numpy as np
from scipy import ndimage as ndi
from sklearn.datasets import load_sample_image
from sklearn.decomposition import PCA

P = 14                                                   # patch size of the DINOv2 ViT-*/14 models
img1 = load_sample_image("china.jpg").mean(2)[:420, :630] / 255.0
# A "copy": rotated by 8 degrees, scaled by 0.92, with an inverted, non-linear tone curve
a, s = np.deg2rad(8), 0.92
M = np.array([[np.cos(a), -np.sin(a)], [np.sin(a), np.cos(a)]]) / s
off = np.array(img1.shape) / 2 - M @ (np.array(img1.shape) / 2)
img2 = 1.0 - ndi.affine_transform(img1, M, offset=off, order=1, mode="nearest") ** 0.6

def patch_features(img, bins=4):
    """Stand-in for ViT patch tokens: orientation histograms (modulo pi, so contrast-invariant)
    pooled in each 14 x 14 patch at two scales, concatenated over the 3 x 3 neighbouring patches."""
    h, w = img.shape[0] // P, img.shape[1] // P
    chans = []
    for sg in (1.0, 2.0):
        gy, gx = ndi.gaussian_filter(img, sg, order=(1, 0)), ndi.gaussian_filter(img, sg, order=(0, 1))
        mag, th = np.hypot(gx, gy), np.arctan2(gy, gx)
        for k in range(bins):
            ch = mag * np.cos(th - np.pi * k / bins) ** 2              # axial orientation bin
            chans.append(ch[:h * P, :w * P].reshape(h, P, w, P).mean((1, 3)))
    f = np.stack(chans, -1)                                            # (h, w, 2 * bins)
    f = np.concatenate([np.roll(f, (dy, dx), (0, 1)) for dy in (-1, 0, 1) for dx in (-1, 0, 1)], -1)
    f = np.sqrt(f.reshape(h * w, -1))
    f -= f.mean(0)
    return f / np.linalg.norm(f, axis=1, keepdims=True), (h, w)

f1, (h, w) = patch_features(img1)
f2, _ = patch_features(img2)
S = f1 @ f2.T                                            # cosine similarities of all patch pairs
nn12, nn21 = S.argmax(1), S.argmax(0)
mutual = np.flatnonzero(nn21[nn12] == np.arange(len(f1)))    # Definition 5.4
centre = lambda k: np.stack([k // w, k % w], 1) * P + P / 2
truth = (centre(np.arange(len(f1))) - off) @ np.linalg.inv(M).T    # where each patch of img1 lands
seen = np.all((truth >= 0) & (truth < np.array(img2.shape)), axis=1)
ok = np.linalg.norm(centre(nn12) - truth, axis=1) < P              # correct within one patch
print(f"plain nearest neighbours: {ok[seen].mean():.0%} correct ({seen.sum()} visible patches)")
print(f"mutual nearest neighbours: {len(mutual)} pairs, {ok[mutual].mean():.0%} correct")

# PCA of the patch features: three components mapped to RGB give a false-colour map
pca = PCA(n_components=3).fit(f1)
rgb = pca.transform(f1).reshape(h, w, 3)
rgb = (rgb - rgb.min((0, 1))) / (rgb.max((0, 1)) - rgb.min((0, 1)))
print("explained variance ratios:", pca.explained_variance_ratio_.round(3), "| false-colour map:", rgb.shape)

6. Image-to-image translation with CNNs

Many conservation questions take an image and ask for another: the sheet without its ink, the colours of a lost strip.

6.1 The U-Net

The U-Net (Ronneberger et al., 2015) has a contracting path of 3 \times 3 convolutions and 2 \times 2 max pooling, and an expanding path of up-convolutions. Skip connections concatenate each encoder level to the decoder level of the same resolution, so fine detail survives the bottleneck. The original network uses unpadded convolutions and maps a 572 \times 572 input to a 388 \times 388 output; large images are processed tile by tile.

6.2 VGG16, Gram matrices and style transfer

VGG16 (Simonyan and Zisserman, 2015) has 13 convolutional layers with 3 \times 3 kernels; its features, trained on ImageNet, serve as a perceptual measure. Gatys et al. (2016) separated content from style with them.

Style is thus matched as a distribution of local features, wherever they occur (Li et al., 2017). The code checks the identity, then synthesises a texture with random convolutional features, which already capture much of a texture (Ustyuzhaninov et al., 2017).

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

torch.set_num_threads(1)
torch.manual_seed(0)

# (1) Squared Gram distance = M^2 MMD^2 with the kernel k(x, y) = (x . y)^2 (Proposition 6.3)
N, M = 16, 200                                    # channels, positions
Fx, Fy = torch.randn(N, M, dtype=torch.float64), 1.3 * torch.randn(N, M, dtype=torch.float64)
G = lambda F: F @ F.T
k = lambda A, B: (A.T @ B) ** 2                   # kernel between the columns (positions)
lhs = ((G(Fx) - G(Fy)) ** 2).sum()
rhs = k(Fx, Fx).sum() + k(Fy, Fy).sum() - 2 * k(Fx, Fy).sum()
print(f"Gram distance {lhs.item():.6e} = MMD form {rhs.item():.6e}")
print("Gram matrix invariant to shuffling the positions:", torch.allclose(G(Fx), G(Fx[:, torch.randperm(M)])))

# (2) Texture synthesis by Gram matching, with random convolutional features
ex = torch.tensor(load_sample_image("flower.jpg")[100:228:2, 200:328:2] / 255.0,
                  dtype=torch.float32).permute(2, 0, 1)[None]        # 64 x 64 colour exemplar
layers = [torch.nn.Conv2d(3, 32, 5, padding=2), torch.nn.Conv2d(32, 32, 5, padding=2),
          torch.nn.Conv2d(32, 32, 5, padding=2)]
for l in layers:
    l.requires_grad_(False)

def grams(x):
    out = []
    for i, l in enumerate(layers):
        x = torch.relu(l(x if i == 0 else torch.nn.functional.avg_pool2d(x, 2)))
        f = x.flatten(2)[0]
        out.append(f @ f.T / f.shape[1])
    return out

target = grams(ex)
img = torch.rand_like(ex).requires_grad_(True)
opt = torch.optim.Adam([img], lr=0.03)
style = lambda x: sum(((g - t) ** 2).sum() for g, t in zip(grams(x), target))
print(f"style loss at start: {style(img).item():.4f}")
for step in range(250):
    opt.zero_grad()
    loss = style(img)
    loss.backward(); opt.step()
    img.data.clamp_(0, 1)
centred = lambda t: (t - t.mean((2, 3), keepdim=True)).detach().flatten().numpy()
corr = np.corrcoef(centred(img), centred(ex))[0, 1]
print(f"style loss after 250 steps: {style(img).item():.4f}; "
      f"pixel correlation with the exemplar (channel means removed): {corr:.2f}")
print(f"channel means: exemplar {ex.mean((0, 2, 3)).numpy().round(2)}, synthesis {img.detach().mean((0, 2, 3)).numpy().round(2)}")

6.3 Case study: the Magic Eraser

In transmitted light, a sheet shows its watermark, chain lines and laid lines, which date and locate the paper; ink drawn on it hides them. By Proposition 1.6, in the log domain \log I = \log T_{\text{paper}} - d_{\text{ink}}, \qquad d_{\text{ink}} \ge 0 . A U-Net trained on synthetic pairs, clean sheets with simulated strokes, predicts the non-negative d_{\text{ink}}; adding it back erases the ink. In one channel, a bright watermark under ink looks like thinner ink; the shape of strokes, or a coloured ink, resolves the ambiguity. The revealed watermark is compared with repertories such as Briquet’s (1907).

import numpy as np
import torch
import torch.nn as nn
from scipy import ndimage as ndi

torch.set_num_threads(1)
torch.manual_seed(0)
rng = np.random.default_rng(0)
S = 40
tint = np.log([0.95, 0.90, 0.80])[:, None, None]          # colour of the paper in transmitted light
a_ink = np.array([0.6, 1.0, 1.6])[:, None, None]          # optical density of a brown ink per channel

def sheet():
    """Log transmitted light (R, G, B) of laid paper with a watermark, covered by pen strokes."""
    yy, xx = np.mgrid[0:S, 0:S]
    paper = 0.6 + 0.03 * np.sin(2 * np.pi * yy / 2.2)                       # laid lines
    paper += 0.08 * (np.abs((xx + rng.integers(50)) % 50 - 25) < 1)         # a chain line
    cy, cx, r = rng.uniform(12, 28), rng.uniform(12, 28), rng.uniform(5, 10)
    ring = np.abs(np.hypot(yy - cy, xx - cx) - r) < 0.8                     # watermark: thinner paper
    paper += 0.12 * ring + 0.02 * ndi.gaussian_filter(rng.standard_normal((S, S)), 3)
    ink = np.zeros((S, S))
    for _ in range(rng.integers(2, 5)):                                     # quadratic Bezier strokes
        t = np.linspace(0, 1, 200)[:, None]
        p0, p1, p2 = rng.uniform(0, S, (3, 2))
        cv = np.clip(((1 - t) ** 2 * p0 + 2 * t * (1 - t) * p1 + t ** 2 * p2).astype(int), 0, S - 1)
        ink[cv[:, 0], cv[:, 1]] = 1
    ink = ndi.gaussian_filter(ndi.binary_dilation(ink, iterations=rng.integers(1, 3)) * 1.0, 0.7)
    clean = np.log(paper)[None] + tint
    seen = clean - a_ink * ink + 0.01 * rng.standard_normal((3, S, S))     # Beer-Lambert: additive in log
    return seen, clean, ring, ink > 0.3

def unet(cin, c=8):
    """A two-level U-Net whose output adds back a non-negative ink density."""
    block = lambda i, o: nn.Sequential(nn.Conv2d(i, o, 3, padding=1), nn.ReLU(), nn.Conv2d(o, o, 3, padding=1), nn.ReLU())
    class UNet(nn.Module):
        def __init__(self):
            super().__init__()
            self.e1, self.e2, self.mid = block(cin, c), block(c, 2 * c), block(2 * c, 4 * c)
            self.u2, self.d2 = nn.ConvTranspose2d(4 * c, 2 * c, 2, stride=2), block(4 * c, 2 * c)
            self.u1, self.d1 = nn.ConvTranspose2d(2 * c, c, 2, stride=2), block(2 * c, c)
            self.out = nn.Conv2d(c, cin, 1)
            nn.init.zeros_(self.out.weight); nn.init.constant_(self.out.bias, -1.0)
        def forward(self, x):
            e1 = self.e1(x)
            e2 = self.e2(nn.functional.max_pool2d(e1, 2))
            m = self.mid(nn.functional.max_pool2d(e2, 2))
            d2 = self.d2(torch.cat([self.u2(m), e2], 1))                   # skip connection
            d1 = self.d1(torch.cat([self.u1(d2), e1], 1))                  # skip connection
            return x + nn.functional.softplus(self.out(d1))                # ink only darkens: add back >= 0
    return UNet()

train = [sheet() for _ in range(2400)]
test = [sheet() for _ in range(50)]

def visibility(img, b):
    """Watermark contrast-to-clutter ratio on the inked pixels (log units)."""
    on, off = img[b[2] & b[3]], img[~b[2] & b[3]]
    return (on.mean() - off.mean()) / off.std()

model = unet(3)
opt = torch.optim.Adam(model.parameters(), lr=3e-3)
for step in range(150):
    batch = train[16 * step: 16 * step + 16]
    x = torch.tensor(np.stack([b[0] for b in batch]), dtype=torch.float32)
    y = torch.tensor(np.stack([b[1] for b in batch]), dtype=torch.float32)
    loss = nn.functional.mse_loss(model(x), y)
    opt.zero_grad(); loss.backward(); opt.step()
with torch.no_grad():
    out = model(torch.tensor(np.stack([b[0] for b in test]), dtype=torch.float32)).numpy()

keep = [i for i, b in enumerate(test) if (b[2] & b[3]).sum() > 3]      # watermark crosses the ink
for name, imgs in [("inked input", [b[0] for b in test]), ("after erasing", out), ("clean paper", [b[1] for b in test])]:
    g = [imgs[i].mean(0) for i in range(len(test))]                     # grey level, log units
    err = np.mean([np.sqrt(np.mean((g[i] - test[i][1].mean(0))[test[i][3]] ** 2)) for i in keep])
    vis = np.mean([visibility(g[i], test[i]) for i in keep])
    print(f"{name:13s}: RMS error under the ink {err:.3f}, watermark visibility {vis:5.2f}")

After 150 steps on a CPU, the error under the ink drops tenfold, but the watermark is still far less visible than on the clean sheet: closing that gap is Exercise 19.

6.4 Case study: the Night Watch

In 1715 Rembrandt’s Night Watch was cut down on all sides, most of all on the left, to fit a wall of Amsterdam’s town hall. A smaller seventeenth-century copy attributed to Gerrit Lundens, in the National Gallery, London, shows the full composition. The Rijksmuseum showed reconstructed strips with the painting in 2021. A reconstruction of this kind has four steps.

1. Register the copy to the painting on the region both show (ch. 2); the copy’s proportions are not exact, so the warp must be local.
2. On that overlap \Omega_o, learn a translation from the copy’s handling and colours to Rembrandt’s (ch. 6.1–6.2).
3. Apply it on the missing strips \Omega_m.
4. Validate on strips of \Omega_o held out from training, with the errors reported.

The result is a hypothesis about the lost parts, not a restoration: the network extrapolates from \Omega_o to \Omega_m.

7. Complements: optics of paint, information and uncertainty

This self-study chapter deepens four points of the course: why pigments mixed in one layer do not mix linearly, how to compare modalities without learning, why contrastive learning needs large batches, and how to attach honest error bars to a reconstruction.

7.1 Kubelka–Munk: the optics of a paint layer

Proposition 1.4 followed light through the paint once. In a turbid layer, light is absorbed and scattered many times. The two-flux model of Kubelka and Munk (1931) follows a downward flux i(z) and an upward flux j(z) at depth z, with absorption coefficient K and scattering coefficient S: \frac{di}{dz} = -(K + S)\, i + S\, j, \qquad \frac{dj}{dz} = (K + S)\, j - S\, i .

The reflectance of a mixture is therefore not a mixture of reflectances. But the rule is linear in c once q = (K/S)_{\text{mix}} is computed from the measured reflectance: \sum_i c_i\, \big(K_i(\lambda) - q(\lambda)\, S_i(\lambda)\big) = 0 \quad \text{for every } \lambda . With \sum_i c_i = 1, non-negative least squares recovers the concentrations. The code compares the two unmixings for a strong blue in a white.

import numpy as np
from scipy.optimize import nnls

lam = np.linspace(400, 700, 61)
# Absorption K and scattering S (per unit concentration) of two synthetic pigments
K = np.stack([0.02 + 0 * lam,                                     # white: weak absorption
              0.05 + 2.0 / (1 + np.exp(-(lam - 560) / 15))], 1)   # blue: absorbs red light
S = np.stack([1.0 + 0 * lam, 0.3 + 0 * lam], 1)                   # white scatters strongly

def R_inf(k_over_s):
    """Reflectance of an opaque layer (Proposition 7.1)."""
    return 1 + k_over_s - np.sqrt(k_over_s**2 + 2 * k_over_s)

c_true = np.array([0.8, 0.2])                                     # concentrations
r_mix = R_inf((K @ c_true) / (S @ c_true))                        # mixing rule (Definition 7.2)
R_pure = R_inf(K / S)                                             # spectra of the pure pigments

# (a) Linear unmixing of reflectances, as in Definition 1.5
w = 100.0
a_lin, _ = nnls(np.vstack([R_pure, w * np.ones(2)]), np.append(r_mix, w))
# (b) Unmixing in K/S space: sum_i c_i (K_i - q S_i) = 0 with q = (K/S) of the mixture
q = (1 - r_mix) ** 2 / (2 * r_mix)
A = K - q[:, None] * S
c_km, _ = nnls(np.vstack([A, w * np.ones(2)]), np.r_[np.zeros(len(lam)), w])
print("true concentrations   :", c_true)
print("linear unmixing       :", a_lin.round(3))
print("Kubelka-Munk unmixing :", c_km.round(3))
print(f"reflectance at 650 nm: white {R_pure[50, 0]:.3f}, blue {R_pure[50, 1]:.3f}, mixture {r_mix[50]:.3f}")

A fifth of blue pulls the reflectance at 650 nm from 0.82 to 0.38, so linear unmixing reads about half blue. In the infrared, both K and S of most paints decrease, so a layer of given thickness hides less of what lies below: this is the physics behind ch. 1.2.

7.2 Mutual information as a similarity

So mutual information, estimated from the joint histogram of the intensities of two overlapping images, tolerates any one-to-one intensity relation, monotone or not. Registration by maximisation of mutual information over the transform parameters was introduced by Viola and Wells (1997) and Maes et al. (1997). It needs no training data, but its histograms are noisy on small tiles, where the learned filters of Definition 2.11 are more precise.

7.3 Contrastive learning and mutual information

The bound cannot exceed \log N, so estimating large dependencies requires many negatives: CLIP used batches of 32,768 pairs. Duplicates in a batch (Proposition 5.2) raise the floor of the loss without saying anything about the encoders.

7.4 Honest error bars for reconstructions

A reconstruction such as that of ch. 6.4 is validated on held-out strips. Split conformal prediction turns those strips into a guaranteed interval.

The guarantee is marginal, and it rests on exchangeability. Neighbouring pixels are dependent, and a missing strip differs from the overlap; calibrate on whole held-out blocks, and report the interval as conditional on that assumption (Vovk et al., 2005; Lei et al., 2018).

8. Research directions and open problems

This self-study chapter surveys active questions, with entry points in the literature. Each item names a problem that the methods of the course do not yet solve.

8.1 Registration and matching

• Detector-free matchers, such as LoFTR (Sun et al., 2021) and RoMa (Edstedt et al., 2024), and fast learned matchers, such as LightGlue (Lindenberger et al., 2023), have changed matching between photographs.
• Open problem: matching between modalities that share few edges, such as element maps and photographs, at gigapixel scale and with sub-pixel accuracy, with reliable estimates of the local error.

8.2 From spectra and element maps to materials

• Reflectance spectra, element maps and radiographs see the same layers through different physics. Hyperspectral unmixing is reviewed by Bioucas-Dias et al. (2012); the physics of mixtures is that of ch. 7.1.
• Open problems: joint inversion of RIS and MA-XRF into pigment maps with uncertainty, and models of layering, where the order of layers matters as much as their content.

8.3 Microstructure

• Convolutional networks now detect cracks in paintings, as in the work of Sizyakin et al. (2020).
• Open problems: a typology of craquelure learned from data and related to supports and schools, as Bucklow (1997) did by hand; and standards of statistical evidence for canvas matches, since profiles are autocorrelated and searches over shifts are multiple comparisons (Example 4.6).

8.4 Foundation models and digitised collections

• Work by Shen, Efros and Aubry (2019) discovers repeated motifs in art collections with spatially consistent features. Elgammal et al. (2018) studied how networks trained on style labels organise art history. Mensink and van Gemert (2014) built a recognition benchmark from the Rijksmuseum collection.
• Open problems: the domain gap between web photographs and paintings, drawings or micrographs; fine-tuning on small expert-labelled sets; and search interfaces whose semantic directions are checked by art historians, not only by retrieval scores.

8.5 Generative reconstruction

• Inpainting has moved from partial convolutions (Liu et al., 2018) to diffusion models, such as RePaint (Lugmayr et al., 2022).
• Open problems: calibrated uncertainty for reconstructed regions (ch. 7.4), and conventions to disclose what is measured, what is inferred and what is invented in an image shown to the public.

8.6 Attribution and connoisseurship

• Tummers (2011) examines how connoisseurs authenticate paintings by Rembrandt and his contemporaries. Computational studies have proposed statistical signatures of artists: Lyu, Rockmore and Farid (2004) used wavelet statistics, and Johnson et al. (2008) compared several image-processing approaches to attributing paintings to Van Gogh from his brushwork.
• Open problems: small data sets, workshops where several hands painted one work, deliberate forgeries, and the communication of evidence as calibrated probabilities that a connoisseur can weigh against other evidence.

Exercises

1. ★ The size of a multimodal stack. A panel of 80 \times 60 cm is imaged in VIS at 10 µm per pixel (3 channels, 16 bits), in IRR at 50 µm (1 channel, 16 bits), in RIS at 0.25 mm (300 bands, 16 bits) and in MA-XRF at 0.5 mm (spectra of 4096 channels of 32 bits). (a) Compute the size of each image. (b) How many 64 \times 64 tiles does the VIS image hold?

   Solution: (a) VIS: 80\,000 \times 60\,000 \times 6 bytes = 28.8 GB. IRR: 16\,000 \times 12\,000 \times 2 = 0.38 GB. RIS: 3200 \times 2400 \times 600 = 4.6 GB. MA-XRF: 1600 \times 1200 \times 16\,384 = 31.5 GB. The stack holds about 65 GB.

   (b) 4.8 \times 10^9 / 4096 \approx 1.17 million tiles: every per-tile computation must be batched (ch. 2.5).

2. ★ Moseley’s law and a modern pigment. (a) Estimate the K\alpha energies of titanium (Z = 22) and zinc (Z = 30) with Proposition 1.7, and compare with the tabulated 4.51 and 8.64 keV. (b) An MA-XRF scan of a seventeenth-century painting shows a strong titanium signal confined to a small patch. What does it suggest, and what should be checked? (c) How would you tell an arsenic pigment from lead white?

   Solution: (a) 10.2 \text{ eV} \times 21^2 \approx 4.50 keV and 10.2 \text{ eV} \times 29^2 \approx 8.58 keV: within 1 %.

   (b) Titanium white is a pigment of the twentieth century, so a confined signal suggests a later retouching. Check its shape against the photograph, since traces of titanium also occur in some natural earths. (c) As K\alpha and Pb L\alpha coincide near 10.5 keV (Example 1.8); compare the As K\beta (11.73 keV) and Pb L\beta (12.61 keV) maps.

3. ★ Phase shifting by hand. At one pixel the three fringe images of Proposition 1.9 give 0.9, 0.3 and 0.6. (a) Compute \varphi, A and B. (b) The reference phase there is 0, with p = 2 mm and \theta = 30^\circ: compute the height. (c) Which heights can be measured without phase unwrapping?

   Solution: (a) By Proposition 1.9, \varphi = \operatorname{atan2}(\sqrt3 \times 0.3,\ 0.9) = \pi/6 . A is the mean, 0.6, and 3B\cos\varphi = 0.9 gives B \approx 0.346.

   (b) The height is h = \frac{\varphi\, p}{2\pi\tan\theta} = \frac{(\pi/6) \times 2}{2\pi \times 0.577} \approx 0.289 \text{ mm}.

   (c) Without unwrapping, |\Delta\varphi| < \pi, that is |h| < p / (2\tan\theta) \approx 1.73 mm.

4. ★ A RANSAC budget. (a) A quarter of the matches between an X-ray and a photograph are correct. How many draws give a 99 % chance of an all-inlier sample for a homography? (b) And for an affine transform (s = 3)? (c) Why is a fixed inlier threshold risky on a large painting photographed with a wide-angle lens?

   Solution: (a) \log 0.01 / \log(1 - 0.25^4) \approx 1176.6, so 1177 draws (Proposition 2.6). (b) \log 0.01 / \log(1 - 0.25^3) \approx 292.4, so 293 draws.

   (c) Residual lens distortion makes the errors of correct matches grow towards the corners (Example 2.5). A threshold tuned at the centre rejects good matches at the edges, which motivates the σ-marginalised weights of Proposition 2.7.

5. ★ When correlation fails. (a) Prove Proposition 2.9. (b) Let a = (1, 2, 3, 4) and b = (a - 2.5)^2 componentwise. Compute \mathrm{NCC}(a, b). (c) Which representation of a and b makes them correlate?

   Solution: (a) Centring removes \beta; the factor \alpha appears in the numerator and |\alpha| in the denominator.

   (b) The centred vectors are (-1.5, -0.5, 0.5, 1.5) and (1, -1, -1, 1), whose inner product is 0: \mathrm{NCC} = 0 although b is a function of a. (c) A representation unchanged by the map, such as the magnitude of a smoothed gradient on images, as in the code of ch. 2.5, or mutual information (ch. 7.2).

6. ★ Thread density from a spectral peak. A 64 \times 64 tile of a radiograph at 0.1 mm per pixel is zero-padded to 512 points. (a) The spectral peak is at index 72: compute the thread density. (b) Compare the bin width with and without padding. (c) The spectrum at indices 71, 72, 73 is 0.8, 1.0, 0.9: refine the estimate.

   Solution: (a) f = 72/512 cycles per pixel, so \rho = 0.1406 / 0.01 \approx 14.06 threads per cm (Definition 4.4). (b) 0.195 threads per cm with padding, 1.56 without.

   (c) The parabola gives \delta = 0.5 \times (0.8 - 0.9)/(0.8 - 2 + 0.9) \approx 0.167, so \rho \approx 72.17/512/0.01 \approx 14.10 threads per cm.

7. ★ Temperature and retrieval. A query has cosine similarities 0.31, 0.29, 0.25 and 0.10 with four images. (a) Compute the softmax probabilities for \tau = 0.01 and \tau = 0.1. (b) Does the temperature change the ranking returned by a search engine? (c) What does it change in training?

   Solution: (a) For \tau = 0.01: 0.879, 0.119, 0.002, 0.000. For \tau = 0.1: 0.402, 0.329, 0.220, 0.049.

   (b) No: the softmax is increasing in each similarity. (c) In the loss of Definition 5.1, a small \tau concentrates the gradient on the hardest negatives; CLIP learns \tau (the code starts from 0.07).

8. ★ Tiling a gigapixel image. (a) Compute the receptive field of two 3 \times 3 convolutions, a 2 \times 2 pooling, then two 3 \times 3 convolutions (ch. 0.6). (b) A U-Net with two poolings and a receptive field of 93 pixels runs on tiles of 512 pixels. By Proposition 6.1, what core is kept per tile, and how many tiles cover a 40\,000 \times 30\,000 image? (c) What is the computational overhead?

   Solution: (a) r = 1 + 2 + 2 + 1 + 2 \times 2 + 2 \times 2 = 14.

   (b) The margin is \rho = 46, so the core is 512 - 92 = 420 pixels, and the number of tiles is \lceil 40\,000/420 \rceil \times \lceil 30\,000/420 \rceil = 96 \times 72 = 6912 . The step of 420 pixels is a multiple of the total stride s = 4, as Proposition 6.1 requires. With four poolings (s = 16), the step drops to 416 pixels and the count rises to 97 \times 73 = 7081 tiles.

   (c) (512/420)^2 \approx 1.49: about 49 % more computation than one pass over the image.

9. ★ Mutual nearest neighbours. The cosine similarities between three patches of a painting (rows) and three patches of a print after it (columns) are S = \begin{pmatrix} 0.9 & 0.2 & 0.4 \\ 0.8 & 0.3 & 0.1 \\ 0.1 & 0.6 & 0.5 \end{pmatrix}. (a) Find the correspondences of Definition 5.4. (b) Why is the second patch left unmatched, and why is that desirable?

   Solution: (a) Row maxima: 1→1, 2→1, 3→2. Column maxima: 1→1, 2→3, 3→3. The mutual pairs are (1, 1) and (3, 2).

   (b) Patch 2 prefers column 1, which prefers patch 1: two patches compete for one partner. Discarding such one-sided matches removes most errors, as the code of ch. 5.5 shows (61 % against 82 % correct).

10. ★★ The homography of a plane. (a) Write the two linear equations in H given by one match, and explain why four matches in general position determine H. (b) Deduce from Proposition 2.3 that a camera rotating about its centre also gives a homography, whatever the scene. (c) Show that H is affine when K = K', R = I, the plane is n = (0, 0, 1)^\top and the translation is parallel to it.

    Solution: (a) With x' = (u, v, 1) and rows h_1, h_2, h_3 of H: h_1^\top x - u\, h_3^\top x = 0 and h_2^\top x - v\, h_3^\top x = 0. H has 8 degrees of freedom (scale is free), and four points, no three collinear, give 8 independent equations.

    (b) With t = 0, H = K' R K^{-1} does not depend on the plane: every point maps the same way.

    (c) Here H = I + K t n^\top K^{-1} / d . With n = e_3 and t_3 = 0, the third entry of K t is zero, because the last row of K is e_3^\top. So H keeps the last row (0, 0, 1) of the identity: it is affine.

11. ★★ The σ-marginalised weight. (a) Check Proposition 2.7 numerically for r = 1 and \sigma_{\max} = 2 by quadrature. (b) Compute the weight w(r) = \operatorname{erfc}\big(r / (\sqrt2\,\sigma_{\max})\big) at r equal to one, two and three times \sigma_{\max}, and compare with a hard threshold at \sigma_{\max}. (c) Show that p(r) integrates to 1.

    Solution: (a) Both sides equal 0.3867.

    (b) 0.317, 0.046 and 0.003: matches just beyond the scale still count a little, far ones not at all, and no single threshold has to be chosen.

    (c) With u = r/(\sqrt2\,\sigma_{\max}) and \int_0^\infty \operatorname{erfc}(u)\,du = 1/\sqrt\pi, the integral is \frac{1}{\sigma_{\max}} \sqrt{\frac{\pi}{2}} \times \frac{\sqrt2\,\sigma_{\max}}{\sqrt\pi} = 1 .

12. ★★ Dimension and surprisal. (a) Minimise c_1 h^4 + c_2/(n h^d) in h and give the resulting rate (Theorem 3.2). (b) By which factor must n grow to halve the mean squared error for d = 4, and for the d = 300 bands of a RIS cube? What follows for surprisal maps of spectra? (c) Prove Proposition 3.4.

    Solution: (a) h^\star = \big(d\, c_2 / (4 c_1 n)\big)^{1/(d+4)}, and the error decreases as n^{-4/(d+4)}.

    (b) By 2^{(d+4)/4}: a factor 4 for d = 4, and 2^{76} for d = 300. Spectra must first be reduced to a few components (PCA or unmixing, Definition 1.5).

    (c) Given U = u, V is normal with mean m(u), linear in u, and variance s^2; take minus the logarithm of its density.

13. ★★ Classifying crack networks. (Python) In the code of ch. 4.1: (a) compute three features per image, the crack density, the orientation coherence of Definition 4.3 and the logarithm of the mean island area (label the connected components of the complement of the mask with scipy.ndimage.label). (b) Generate 30 networks with stretch 1 and 30 with stretch 2, and estimate the 5-fold cross-validated accuracy of a logistic regression on each feature alone and on all three. (c) Repeat with stretch 1.5, and say which feature carries the classification.

    Hints and expected result: (a)–(b) With stretch 2, coherence alone separates the classes (accuracy about 1.0, mean coherence about 0.12 against 0.52). Density reaches about 0.97 and island area about 0.7.

    (c) With stretch 1.5, coherence still gives about 0.98, while density falls to about 0.78 and island area to near chance (about 0.57). Real crack types differ in more than anisotropy: width, sharpness, junctions, support and restorations vary from painting to painting. Validate on whole held-out paintings, as Exercise 18 explains.

14. ★★ How significant is a canvas match? (Python) In the code of ch. 4.2: (a) apply Fisher’s transform of Example 4.6 to the same-roll correlation, over its n = 34 overlapping tiles at the best shift, and to the largest unrelated one, with the minimum overlap n = 24. (b) Estimate the lag-one autocorrelation \rho_1 of the profiles and the effective sample size. (c) Replace the 20 unrelated rolls by 200, and report an empirical p-value.

    Hints and expected result: (a) z \approx 1.832 \times \sqrt{31} \approx 10.2 for r = 0.95, but also z \approx 5.2 for the chance correlation of 0.81: the naive test would call both significant.

    (b) The profiles have \rho_1 \approx 0.8. For two independent AR(1)-like profiles with lag-one autocorrelations \rho_1 and \rho_1', Bartlett’s approximation gives n_{\text{eff}} \approx n\, \frac{1 - \rho_1 \rho_1'}{1 + \rho_1 \rho_1'} , about 8 of the 34 tiles here.

    (c) With k exceedances among N unrelated rolls, the empirical p-value is (1 + k)/(1 + N); with none in 200, report p \le 1/201 \approx 0.005.

15. ★★ The floor of a contrastive loss. (Python) In the code of ch. 5.3, the image cannot know which of the three templates its caption used. (a) With 20 classes (digit and ground) in batches of 256, estimate the floor of Proposition 5.2 and compare with the printed final loss. (b) Prove the analogue of Proposition 5.2 for the text-to-image term. (c) Modify the loss so that all captions of the same class count as positives, and observe the new final loss.

    Hints and expected result: (a) Each image has on average 1 + 255/20 \approx 13.75 indistinguishable captions in the batch, its own included, so the floor is about \log 13.75 \approx 2.6; the code prints 2.601.

    (b) Images of one class with the same caption are interchangeable for the text encoder, so the same argument applies by columns.

    (c) The answer depends on the form of the loss, with p_{ij} the softmax probabilities of each row and P(i) the captions of its class. With -\log \sum_{j \in P(i)} p_{ij}, the probability mass on all positives, the floor is 0 and the loss approaches it (about 0.002 after 400 steps). With the average of -\log p_{ij} over P(i), the floor stays \log m (about 2.59). The retrieval metrics barely change in either case.

16. ★★ What the Gram matrix sees. (a) Prove \langle a a^\top, b b^\top \rangle_F = (a^\top b)^2 and deduce Proposition 6.3. (b) For a single feature channel, show that the style loss only compares the mean of f_m^2. (c) Show that the Gram matrix of centred features is M times their covariance, and relate this to matching channel means and variances, as adaptive instance normalisation does (Huang and Belongie, 2017).

    Solution: (a) For vectors a and b, \langle a a^\top, b b^\top \rangle_F = \operatorname{tr}(a a^\top b b^\top) = (a^\top b)^2 ; expand the squared norm of \sum_m f_m f_m^\top - \sum_m p_m p_m^\top with it.

    (b) For N = 1, G = \sum_m f_m^2, a scalar: two images with the same energy have the same style, whatever their arrangement.

    (c) The Gram matrix of centred features is \sum_m (f_m - \bar f)(f_m - \bar f)^\top = M\, \widehat{\mathrm{Cov}} . AdaIN matches only the diagonal (means and variances per channel), a cheaper and weaker version of Gram matching.

17. ★★ Kubelka–Munk unmixing. (Python) (a) From Proposition 7.1, compute K/S for R_\infty = 0.82, 0.38 and 0.064, and check against the code of ch. 7.1. (b) Explain why linear unmixing overestimates the blue. (c) Add 1 % noise and a third pigment to the code, and compare the errors of the two unmixings over 100 noisy draws.

    Hints and expected result: (a) 0.0198, 0.506 and 6.84; the blue at 650 nm has K/S = 2.045/0.3 \approx 6.8, and the mixture 0.425/0.86 \approx 0.49.

    (b) The mixture’s reflectance is far below the average of the pure reflectances, since absorption adds up in K/S, not in R; linear unmixing explains the darkness with too much blue. (c) The Kubelka–Munk unmixing stays close to the truth; the linear one keeps a bias that noise does not remove.

18. ★★★ Calibrated evidence for an attribution. A classifier trained on tiles of authenticated paintings by an artist and by contemporaries gives scores to the tiles of a disputed painting. (a) Explain why the accuracy on random held-out tiles overstates the evidence. (b) Propose a validation and a calibration that turn the scores into a likelihood ratio for the whole painting. (c) Discuss how this ratio could be weighed with other evidence (Tummers, 2011).

    Hints: (a) Tiles of one painting share support, ground, condition and photography: the classifier can learn the painting, not the hand. (b) Leave one painting out; aggregate the tile scores per painting; fit the score distributions under both hypotheses on held-out paintings, and report, with its uncertainty, the likelihood ratio \mathrm{LR} = \frac{p(\text{score} \mid \text{artist})}{p(\text{score} \mid \text{other})} .

    (c) In Bayesian terms, posterior odds are the likelihood ratio times the prior odds formed by provenance, documents and connoisseurship. Discuss workshop participation, the small number of works, and the risk of deliberately deceptive works.

19. ★★★ Closing the gap of the Magic Eraser. (Python) In the code of ch. 6.3, the error under the ink drops tenfold, yet the watermark visibility stays far below that of the clean sheet. (a) List the causes you can test. (b) Train longer, and weight the loss on inked pixels; report the visibility over several seeds. (c) Propose a measurement that would make the separation of ink and watermark easier on real sheets.

    Hints and expected result: (a) The short training, errors at the edges of strokes, and the ambiguity between thin ink and thin paper.

    (b) In a test run, 600 steps instead of 150 raised the visibility from 1.4 to about 1.9, with an error of 0.054 under the ink; weighting the inked pixels five times gave about 2.1 and 0.046. The clean sheet is at 4.6. Report means and spreads over at least five seeds.

    (c) A second image in which the ink is more transparent, for instance in the infrared, where iron-gall inks are known to fade, or a raking-light image that shows the texture of the paper; then feed all channels to the network.

Projects and assessment

The official assessment is by projects, submitted as code and documented and demonstrated in Jupyter notebooks (master-mva.com). Université Paris-Saclay describes it as “Evaluation Continue non Intégrale”; no weights are displayed. The subjects below follow that format: a notebook that runs, explains and demonstrates, with a short report and a talk. Each reproduces or extends published work, in the MVA tradition. No written exam is announced; the exam proposal is for self-assessment.

Solo projects

1. [Solo · ≈ 20 h] Paper study: MAGSAC++ (Barath et al., 2020) against tuned RANSAC. Does marginalising over the noise scale beat a tuned threshold when the noise of correct matches varies across the image?
   • Implement the normalised DLT, RANSAC, the simplified weight of Proposition 2.7 and the MAGSAC++ weights of the paper, with incomplete gamma functions.
   • Data: synthetic matches with radial distortion and noise growing towards the corners; optionally, image pairs of the public HPatches benchmark (Balntas et al., 2017).
   • Success criterion: grid reprojection error within 10 % of an oracle fit on the true inliers, over 20 seeds and all thresholds tried. Deliverables: notebook, 4–6 page report, 10-minute talk. Chapters 0 and 2.
2. [Solo · ≈ 20 h] Paper study: style as a distribution (Gatys et al., 2016; Li et al., 2017). Which feature statistics reproduce the brushwork of a painting?
   • Synthesise textures from crops of public-domain paintings (open-access images of the Rijksmuseum or the Metropolitan Museum of Art) with pretrained VGG16 features from torchvision.
   • Compare Gram matching on several layer sets, matching of channel means and variances only, and random features (Ustyuzhaninov et al., 2017); verify Proposition 6.3 numerically.
   • Success criterion: a figure per method, and an MMD on held-out VGG features ranking the methods. Deliverables: notebook, report, talk. Chapters 5.3 and 6.2.
3. [Solo · ≈ 20 h] Paper study: weave analysis from radiographs (Johnson, Johnson and Erdmann, 2013), on simulated canvases. How long must the common part of two canvases be for a reliable match?
   • Extend the simulator of ch. 4.2 with thread-thickness variations, a slight skew of the weave and stronger noise.
   • Implement thread counting by tile spectra, warp matching over shifts and flips, and a null from 500 unrelated rolls.
   • Success criterion: ROC curves of same-roll detection against the overlap length, and the overlap needed for a 90 % detection rate at a 1 % false-match rate. Deliverables: notebook, report, talk. Chapters 1.4 and 4.2.

Group projects

1. [Group of 2–3 · ≈ 30 h each] A multimodal registration pipeline, end to end. Simulate a stack from an open-access image of a painting: an IRR-like channel with a synthetic underdrawing, element-like maps with Poisson noise, each modality with its own homography, lens distortion and smooth local warp.
   • Member 1: keypoints, descriptors (handcrafted, or HardNet through the Kornia library of Riba et al., 2020), matching, robust homography with distortion (ch. 2.1–2.4).
   • Member 2: patchwise offsets with tile stacks, structure-tensor weights and learned filters (Definition 2.11).
   • Member 3: conditional surprisal on the registered stack (ch. 3.1), as a function of the residual misregistration.
   • Success criterion: a median residual below 0.3 px, and a curve of underdrawing detection (AUC) against registration error. Deliverables: notebooks, a 6–8 page report, a 20-minute talk. Chapters 1–3.
2. [Group of 2–3 · ≈ 30 h each] Paper study: CLIP (Radford et al., 2021) and DINOv2 (Oquab et al., 2024) on an open-access collection. Download a few thousand public-domain images with their metadata from the Metropolitan Museum of Art or the Rijksmuseum, and run the pretrained models on a free GPU tier.
   • Member 1: a search engine with an approximate nearest-neighbour index; precision at 10 on queries built from the metadata.
   • Member 2: semantic directions (Definition 5.3) tested against dates and techniques of the metadata, with an analysis of failures.
   • Member 3: DINOv2 correspondences (Definition 5.4) between works and prints or copies after them, and PCA maps.
   • Success criterion: precision at 10 above a keyword baseline, and a hand-checked set of correspondences with its error rate. Deliverables: notebooks, report, talk. Chapter 5.
3. [Group of 2–3 · ≈ 30 h each] Paper study: U-Net (Ronneberger et al., 2015) for ink removal, after the Magic Eraser case. Can a watermark hidden under a drawing be identified after erasing the ink?
   • Member 1: a simulator of sheets in transmitted light (laid and chain lines, 100 synthetic watermark designs, inks with spectral signatures).
   • Member 2: the U-Net with a non-negative ink density, trained with the ablations of Exercise 19.
   • Member 3: watermark retrieval from the synthetic repertory, and conformal error bars (Proposition 7.6) on held-out sheets.
   • Success criterion: top-1 retrieval of the right watermark under the ink above 80 %. Deliverables: notebooks, report, talk. Chapters 1, 6 and 7.

Written exam (proposal)

Two hours, one handwritten A4 sheet and a calculator allowed; four questions of 5 points.

1. Question (5 points). For visible photography, infrared reflectography, radiography and MA-XRF, state the physical quantity measured, the depth reached and one question about a painting that each modality helps answer. Then explain, with a model, why an underdrawing appears in the infrared. Expected: the table of ch. 1 (reflectance at the surface; reflectance below paint; X-ray transmission dominated by lead white; element fluorescence attenuated with depth), and Proposition 1.4 with a transmittance that grows with the wavelength over a carbon-based drawing.

2. Question (5 points). (a) A third of the matches are correct: how many RANSAC draws give 99 % confidence for a homography? (b) With f = 6000 px and k_1 = -0.01, how far is a point at r = 0.5 displaced? (c) A 128-pixel tile at 0.05 mm per pixel, zero-padded to 1024 points, has its spectral peak at index 72: give the thread density.

   Expected: (a) \log 0.01 / \log(1 - 3^{-4}) \approx 370.7, so 371 draws. (b) 0.01 \times 0.5^3 \times 6000 = 7.5 px. (c) 72/1024 = 0.0703 cycles per pixel, divided by 0.005 cm: about 14.1 threads per cm.

3. Question (5 points). Prove Proposition 2.10 from the shift theorem. Then derive the normal equations of Proposition 2.8, and explain why a tile containing a single straight edge cannot be registered in both directions. Expected: the cross-power spectrum reduces to a pure phase whose inverse DFT is a Kronecker delta; linearisation of A(x - \delta), the system J\delta = -\sum \nabla A\,(B - A), and a rank-one structure tensor when all gradients are parallel.

4. Question (5 points). Prove Proposition 6.3, and deduce that two images whose features are the same up to a permutation of positions have the same style loss. What does this imply for the transfer of style to a painting with a given composition? Expected: \langle a a^\top, b b^\top\rangle_F = (a^\top b)^2, expansion into the three kernel sums, invariance of \sum_m f_m f_m^\top to permutations; style is a distribution of local features, so composition must come from the content loss.

Grading rubric (proposal)

The official grade comes from the projects, in the form of code and notebooks. The weights below are a proposal for grading one project.

References

Official sources and readings of the course

• master-mva.com, course page “Artificial Intelligence and Computer Vision for Cultural Heritage”, 2026–2027: https://www.master-mva.com/cours/ai-and-computer-vision-for-cultural-heritage/.
• Université Paris-Saclay, M2 Mathematics, Vision, Learning, programme table and course detail “AI and Computer vision for cultural heritage” (ECTS, hourly volume, form of assessment).
• A. Tummers, The Eye of the Connoisseur: Authenticating Paintings by Rembrandt and His Contemporaries, J. Paul Getty Museum, 2011 (additional reading named by the course).
• A. Barroso-Laguna, E. Riba, D. Ponsa, K. Mikolajczyk, “Key.Net: keypoint detection by handcrafted and learned CNN filters”, ICCV, 2019.
• D. Barath, J. Noskova, M. Ivashechkin, J. Matas, “MAGSAC++, a fast, reliable and accurate robust estimator”, CVPR, 2020.
• L. A. Gatys, A. S. Ecker, M. Bethge, “Image style transfer using convolutional neural networks”, CVPR, 2016.
• O. Ronneberger, P. Fischer, T. Brox, “U-Net: convolutional networks for biomedical image segmentation”, MICCAI, 2015.
• A. Radford, J. W. Kim, C. Hallacy, A. Ramesh, G. Goh, S. Agarwal, G. Sastry, A. Askell, P. Mishkin, J. Clark, G. Krueger, I. Sutskever, “Learning transferable visual models from natural language supervision”, ICML, 2021.
• M. Oquab, T. Darcet, T. Moutakanni et al., “DINOv2: learning robust visual features without supervision”, Transactions on Machine Learning Research, 2024 (arXiv:2304.07193, 2023).
• D. Mishkin, F. Radenović, J. Matas, “Repeatability is not enough: learning affine regions via discriminability”, ECCV, 2018 (AffNet, named in the syllabus).
• A. Mishchuk, D. Mishkin, F. Radenović, J. Matas, “Working hard to know your neighbor’s margins: local descriptor learning loss”, NeurIPS, 2017 (HardNet, named in the syllabus).
• K. Simonyan, A. Zisserman, “Very deep convolutional networks for large-scale image recognition”, ICLR, 2015 (VGG16, named in the syllabus).

Textbooks

• R. Hartley, A. Zisserman, Multiple View Geometry in Computer Vision, 2nd ed., Cambridge University Press, 2004.
• R. Szeliski, Computer Vision: Algorithms and Applications, 2nd ed., Springer, 2022.
• T. Lindeberg, Scale-Space Theory in Computer Vision, Kluwer Academic Publishers, 1994.
• B. W. Silverman, Density Estimation for Statistics and Data Analysis, Chapman and Hall, 1986.
• C. M. Bishop, Pattern Recognition and Machine Learning, Springer, 2006.
• I. Goodfellow, Y. Bengio, A. Courville, Deep Learning, MIT Press, 2016.
• V. Vovk, A. Gammerman, G. Shafer, Algorithmic Learning in a Random World, Springer, 2005.

Imaging of paintings and works on paper (ch. 1, 3, 4, 6, 7)

• P. Kubelka, F. Munk, “Ein Beitrag zur Optik der Farbanstriche”, Zeitschrift für technische Physik 12, 593–601, 1931.
• J. R. J. van Asperen de Boer, “Reflectography of paintings using an infrared vidicon television system”, Studies in Conservation 14(3), 96–118, 1969.
• J. K. Delaney, J. G. Zeibel, M. Thoury, R. Littleton, M. Palmer, K. M. Morales, E. R. de la Rie, A. Hoenigswald, “Visible and infrared imaging spectroscopy of Picasso’s Harlequin Musician: mapping and identification of artist materials in situ”, Applied Spectroscopy 64(6), 584–594, 2010.
• J. M. Bioucas-Dias, A. Plaza, N. Dobigeon, M. Parente, Q. Du, P. Gader, J. Chanussot, “Hyperspectral unmixing overview: geometrical, statistical, and sparse regression-based approaches”, IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing 5(2), 354–379, 2012.
• J. Dik, K. Janssens, G. Van der Snickt, L. van der Loeff, K. Rickers, M. Cotte, “Visualization of a lost painting by Vincent van Gogh using synchrotron radiation based X-ray fluorescence elemental mapping”, Analytical Chemistry 80(16), 6436–6442, 2008.
• M. Alfeld, K. Janssens, J. Dik, W. de Nolf, G. van der Snickt, “Optimization of mobile scanning macro-XRF systems for the in situ investigation of historical paintings”, Journal of Analytical Atomic Spectrometry 26, 899–909, 2011.
• J. Geng, “Structured-light 3D surface imaging: a tutorial”, Advances in Optics and Photonics 3(2), 128–160, 2011.
• C. Zuo, S. Feng, L. Huang, T. Tao, W. Yin, Q. Chen, “Phase shifting algorithms for fringe projection profilometry: a review”, Optics and Lasers in Engineering 109, 23–59, 2018.
• S. Bucklow, “The description of craquelure patterns”, Studies in Conservation 42(3), 129–140, 1997.
• B. Cornelis, T. Ružić, E. Gezels, A. Dooms, A. Pižurica, L. Platiša, J. Cornelis, M. Martens, M. De Mey, I. Daubechies, “Crack detection and inpainting for virtual restoration of paintings: the case of the Ghent Altarpiece”, Signal Processing 93(3), 605–619, 2013.
• D. H. Johnson, C. R. Johnson Jr., R. G. Erdmann, “Weave analysis of paintings on canvas from radiographs”, Signal Processing 93(3), 527–540, 2013.
• C. R. Johnson Jr., E. Hendriks, I. J. Berezhnoy, E. Brevdo, S. M. Hughes, I. Daubechies, J. Li, E. Postma, J. Z. Wang, “Image processing for artist identification”, IEEE Signal Processing Magazine 25(4), 37–48, 2008.
• C.-M. Briquet, Les Filigranes: dictionnaire historique des marques du papier dès leur apparition vers 1282 jusqu’en 1600, Geneva, 1907.

Registration and local features (ch. 2, 4, 7)

• C. Harris, M. Stephens, “A combined corner and edge detector”, Alvey Vision Conference, 1988.
• D. G. Lowe, “Distinctive image features from scale-invariant keypoints”, International Journal of Computer Vision 60(2), 91–110, 2004.
• T. Lindeberg, “Edge detection and ridge detection with automatic scale selection”, International Journal of Computer Vision 30(2), 1998.
• A. F. Frangi, W. J. Niessen, K. L. Vincken, M. A. Viergever, “Multiscale vessel enhancement filtering”, MICCAI, 1998.
• M. A. Fischler, R. C. Bolles, “Random sample consensus: a paradigm for model fitting with applications to image analysis and automated cartography”, Communications of the ACM 24(6), 381–395, 1981.
• D. Barath, J. Matas, J. Noskova, “MAGSAC: marginalizing sample consensus”, CVPR, 2019.
• R. I. Hartley, “In defense of the eight-point algorithm”, IEEE Transactions on Pattern Analysis and Machine Intelligence 19(6), 580–593, 1997.
• Z. Zhang, “A flexible new technique for camera calibration”, IEEE Transactions on Pattern Analysis and Machine Intelligence 22(11), 1330–1334, 2000.
• B. D. Lucas, T. Kanade, “An iterative image registration technique with an application to stereo vision”, IJCAI, 1981.
• C. D. Kuglin, D. C. Hines, “The phase correlation image alignment method”, Proceedings of the IEEE International Conference on Cybernetics and Society, 1975.
• M. Guizar-Sicairos, S. T. Thurman, J. R. Fienup, “Efficient subpixel image registration algorithms”, Optics Letters 33(2), 156–158, 2008.
• P. Viola, W. M. Wells III, “Alignment by maximization of mutual information”, International Journal of Computer Vision 24(2), 137–154, 1997.
• F. Maes, A. Collignon, D. Vandermeulen, G. Marchal, P. Suetens, “Multimodality image registration by maximization of mutual information”, IEEE Transactions on Medical Imaging 16(2), 187–198, 1997.
• V. Balntas, K. Lenc, A. Vedaldi, K. Mikolajczyk, “HPatches: a benchmark and evaluation of handcrafted and learned local descriptors”, CVPR, 2017.
• E. Riba, D. Mishkin, D. Ponsa, E. Rublee, G. Bradski, “Kornia: an open source differentiable computer vision library for PyTorch”, WACV, 2020.

Embeddings, style and image-to-image networks (ch. 5–7)

• A. van den Oord, Y. Li, O. Vinyals, “Representation learning with contrastive predictive coding”, arXiv:1807.03748, 2018.
• B. Poole, S. Ozair, A. van den Oord, A. A. Alemi, G. Tucker, “On variational bounds of mutual information”, ICML, 2019.
• M. Caron, H. Touvron, I. Misra, H. Jégou, J. Mairal, P. Bojanowski, A. Joulin, “Emerging properties in self-supervised vision transformers”, ICCV, 2021.
• A. Dosovitskiy et al., “An image is worth 16x16 words: transformers for image recognition at scale”, ICLR, 2021.
• J. Johnson, M. Douze, H. Jégou, “Billion-scale similarity search with GPUs”, IEEE Transactions on Big Data 7(3), 535–547, 2021.
• L. A. Gatys, A. S. Ecker, M. Bethge, “A neural algorithm of artistic style”, arXiv:1508.06576, 2015.
• Y. Li, N. Wang, J. Liu, X. Hou, “Demystifying neural style transfer”, IJCAI, 2017.
• I. Ustyuzhaninov, W. Brendel, L. A. Gatys, M. Bethge, “What does it take to generate natural textures?”, ICLR, 2017.
• X. Huang, S. Belongie, “Arbitrary style transfer in real-time with adaptive instance normalization”, ICCV, 2017.
• J. Johnson, A. Alahi, L. Fei-Fei, “Perceptual losses for real-time style transfer and super-resolution”, ECCV, 2016.
• P. Isola, J.-Y. Zhu, T. Zhou, A. A. Efros, “Image-to-image translation with conditional adversarial networks”, CVPR, 2017.
• J. Lei, M. G’Sell, A. Rinaldo, R. J. Tibshirani, L. Wasserman, “Distribution-free predictive inference for regression”, Journal of the American Statistical Association 113(523), 1094–1111, 2018.

Research directions (ch. 8)

• J. Sun, Z. Shen, Y. Wang, H. Bao, X. Zhou, “LoFTR: detector-free local feature matching with transformers”, CVPR, 2021.
• P. Lindenberger, P.-E. Sarlin, M. Pollefeys, “LightGlue: local feature matching at light speed”, ICCV, 2023.
• J. Edstedt, Q. Sun, G. Bökman, M. Wadenbäck, M. Felsberg, “RoMa: robust dense feature matching”, CVPR, 2024.
• R. Sizyakin, B. Cornelis, L. Meeus, H. Dubois, M. Martens, V. Voronin, A. Pižurica, “Crack detection in paintings using convolutional neural networks”, IEEE Access 8, 2020.
• X. Shen, A. A. Efros, M. Aubry, “Discovering visual patterns in art collections with spatially-consistent feature learning”, CVPR, 2019.
• A. Elgammal, B. Liu, D. Kim, M. Elhoseiny, M. Mazzone, “The shape of art history in the eyes of the machine”, AAAI, 2018.
• T. Mensink, J. van Gemert, “The Rijksmuseum Challenge: museum-centered visual recognition”, ICMR, 2014.
• G. Liu, F. A. Reda, K. J. Shih, T.-C. Wang, A. Tao, B. Catanzaro, “Image inpainting for irregular holes using partial convolutions”, ECCV, 2018.
• A. Lugmayr, M. Danelljan, A. Romero, F. Yu, R. Timofte, L. Van Gool, “RePaint: inpainting using denoising diffusion probabilistic models”, CVPR, 2022.
• S. Lyu, D. Rockmore, H. Farid, “A digital technique for art authentication”, Proceedings of the National Academy of Sciences 101(49), 17006–17010, 2004.