| SOCR ≫ | BPAD1 Website ≫ | BPAD GitHub ≫ |
Biomedical Physics with Applications to Disease (2nd Edition)
BPAD2 Chapter 7: General Medical Image Processing, Quantitative Measurement and Uncertainty
Magdalena I. Ivanova, Ivo D. Dinov (BPAD Team + SOCR)
August 06, 2026
How to use this chapter
This chapter is the computational hinge of the book. The preceding chapters derived the physics of image formation for X-ray/CT, MRI, ultrasound, and nuclear medicine. This chapter asks what happens after an image exists: how an array of reconstructed numbers becomes a calibrated, validated measurement with a defensible uncertainty, and how that measurement feeds a model.
Everything in the main teaching thread is implemented in base R on a synthetic phantom whose ground truth we control, so every claim can be checked numerically rather than asserted. Real-image demonstrations (EBImage, NIfTI volumes) are optional and clearly flagged.
Learning objectives
After working through this chapter, a student should be able to:
- Describe the linear imaging forward model \(g = \mathcal{H}f + \eta\) and identify, for a given modality, what \(f\), \(\mathcal{H}\), and \(\eta\) physically represent.
- Convert between image (voxel) coordinates and physical (patient) coordinates using pixel spacing, slice thickness, and the affine geometry matrix; compute areas and volumes with correct units.
- Predict and verify how a linear filter changes both the noise level (\(\sigma_{\text{out}} = \sigma\lVert w\rVert_2\)) and the effective resolution (FWHMs add in quadrature).
- Explain and demonstrate the equivalence of spatial convolution and frequency-domain multiplication, and use it to reason about PSF, MTF, Nyquist sampling, aliasing, and Gibbs ringing.
- Implement and evaluate thresholding, region-growing, clustering, and morphological segmentation, and score the result with overlap (Dice/Jaccard) and boundary (HD95, ASSD) metrics.
- Propagate uncertainty through an image-derived measurement using the GUM framework, and distinguish accuracy (bias) from precision (repeatability) using Monte-Carlo, Bland–Altman, repeatability coefficient, and ICC analyses.
- Critique a quantitative-imaging claim: identify which pipeline step dominates the error budget and what evidence would be needed to trust the number.
Notation
| Symbol | Meaning |
|---|---|
| \(f(x,y,z)\) | true (unknown) object property being imaged |
| \(g\), \(I\) | measured / reconstructed image |
| \(\mathcal{H}\), \(h\) | imaging operator; point spread function (PSF) |
| \(\eta\) | noise, with standard deviation \(\sigma\) |
| \(\Delta x,\Delta y,\Delta z\) | voxel spacing (mm) |
| \(N\) | number of voxels in a region |
| \(w\) | convolution kernel (filter weights) |
| \(F[u,v]\) | 2D discrete Fourier transform of the image |
| \(f_x,f_y\) | spatial frequency (cycles mm\(^{-1}\)) |
| \(f_{\text{Ny}}\) | Nyquist frequency, \(1/(2\Delta x)\) |
| \(M\) | binary mask (segmentation) |
| \(u_c(\cdot)\) | combined standard uncertainty |
| \(\lVert w \rVert_2\) | \(\sqrt{\sum_i w_i^2}\) |
Instructor note. A two-week module maps naturally onto this chapter: week 1 runs from What is a medical image? through Fourier analysis (representation, inspection, preprocessing, filtering, frequency domain), and week 2 from Registration through the Capstone (registration, segmentation, measurement and uncertainty, radiomics, visualization). The “Checkpoint” boxes are designed as 2-minute in-class polls; the Learning activities section supplies laboratory exercises; Answers to selected problems gives worked solutions for the numerical problems.
What is a medical image?
A medical image is a numerical measurement of biological structure or function. Each image is produced by an imaging system that detects a physical signal, reconstructs that signal into a spatial map, and stores the result as numerical data. The picture we see on the screen is a display of those numbers, not the numbers themselves.
The physical source of the signal differs across modalities: computed tomography (CT) measures X-ray attenuation, magnetic resonance imaging (MRI) measures nuclear magnetic resonance signals, positron emission tomography (PET) measures radiotracer activity, ultrasound measures reflected acoustic waves, and microscopy usually measures transmitted or emitted light. Once acquired and reconstructed, however, most biomedical images are processed using a common set of computational steps: import, visualization, denoising, enhancement, filtering, thresholding, segmentation, registration, feature extraction, quantitative modeling, and AI-assisted interpretation.
This shared computational core is the subject of the present chapter. The modality-specific physics of image formation was developed in the preceding chapters. Here we study what happens after an image exists, regardless of how it was made.
| Question | Where it is answered |
|---|---|
| How does a medical image become data? | Image import and representation |
| What information is stored in an image? | Image import and representation (pixels, voxels, metadata, geometry) |
| What steps turn an image into a reliable measurement? | Workflow through Segmentation |
| How is that measurement’s uncertainty quantified? | From segmentation to measurement |
| How do measurements become predictions? | Feature extraction, radiomics, and modeling |
Clinical example. A CT scan looks like an anatomical photograph, but every voxel carries a numerical value (a Hounsfield unit) related to local X-ray attenuation. Those values can be used to measure lesion size, tissue density, calcium burden, and treatment response, quantities a visual reading alone cannot supply.
The imaging forward model
Almost every statement in this chapter can be traced back to one equation. Write \(f\) for the true spatial distribution of the physical property being imaged (attenuation coefficient, spin density, tracer concentration, acoustic impedance, fluorophore concentration) and \(g\) for the reconstructed image. Then
\[\boxed{\;g \;=\; \mathcal{S}\big[\,b\cdot(h * f)\,\big] \;+\; \eta\;}\]
where
- \(h\) is the point spread function (PSF): the system’s response to a point object. Convolution with \(h\) is blur, and it is irreversible in the presence of noise.
- \(b\) is a slowly varying sensitivity / bias field (MRI receive-coil profile, microscope illumination, CT bowtie/scatter residual, PET attenuation and normalization residual).
- \(\mathcal{S}\) is sampling onto a discrete voxel grid, which imposes the Nyquist limit and creates partial-volume mixing.
- \(\eta\) is noise, whose statistics are modality-specific (Poisson, Gaussian, Rician, speckle).
Key idea. Image processing is applied inverse problem solving. Every operation in this chapter either estimates a nuisance term (\(b\), \(\eta\)), partially inverts an operator (\(h\), \(\mathcal{S}\)), or summarizes \(f\) into a measurement. Nothing we do can create information that \(\mathcal{H}\) destroyed: once a spatial frequency is gone from \(g\), no filter recovers it, and any method that appears to do so is supplying a prior, not data.
From signal to image
A medical image typically passes through several stages before it can be analyzed:
- A physical interaction occurs in the body (attenuation, relaxation, emission, scattering, fluorescence).
- A detector records a raw signal (projections, k-space lines, coincidence events, echoes, photon counts).
- The imaging system reconstructs or processes the signal.
- The reconstructed image is stored as an array of numbers plus metadata.
- The image is displayed using a grayscale or color mapping.
- The image is analyzed quantitatively.
Workflow idea. Acquisition creates the raw information. Image processing organizes, cleans, aligns, segments, and measures that information. The two are distinct: better processing cannot recover information the acquisition never captured.
Qualitative interpretation versus quantitative analysis
Radiologists and clinicians often interpret images qualitatively, by visual pattern recognition. Quantitative analysis goes further, extracting reproducible numbers from the image data: tumor diameter or volume, cell counts, vessel diameter, mean intensity inside a region of interest, texture, shape, and change over time.
Key idea. Quantitative imaging requires a reproducible path from image data to numerical measurement. Every step on that path, import, normalization, filtering, segmentation, can change the final number, and every step therefore belongs in the error budget.
Why quantitative analysis is difficult
A measurement extracted from an image can be affected by many factors that have nothing to do with the underlying biology: acquisition settings, scanner calibration, patient or sample motion, noise, the reconstruction algorithm, contrast timing, voxel size, intensity normalization, the segmentation method, and observer variability.
Common pitfall. Two images may look similar but differ in acquisition settings or voxel size. Quantitative measurements should not be compared across images unless those images have been properly calibrated and processed with a documented, consistent pipeline.
Images, volumes, and tensors
It is useful, and increasingly necessary in the era of machine learning, to think of a digital image as a tensor: a multi-dimensional array of numbers indexed by spatial (and possibly temporal and channel) coordinates.
- A grayscale 2D image is a rank-2 tensor (a matrix), \(I(x,y)\).
- A 3D volume (structural MRI, CT, PET) is a rank-3 tensor, \(I(x,y,z)\).
- A multi-channel image (RGB microscopy, multi-sequence MRI) adds a channel index, \(I(x,y,c)\).
- A dynamic series (fMRI, dynamic-contrast MRI, dynamic PET, time-lapse microscopy) adds time, giving a rank-4 tensor, \(I(x,y,z,t)\).
In some modalities the value stored at each voxel is itself a tensor. In diffusion tensor imaging (DTI), each voxel holds a \(3\times 3\) symmetric positive-definite diffusion tensor \(D\) describing the directional mobility of water. Its eigenvalues \(\lambda_1\ge\lambda_2\ge\lambda_3\) give the two most widely used scalar maps,
\[\text{MD} = \frac{\lambda_1+\lambda_2+\lambda_3}{3}, \qquad \text{FA} = \sqrt{\tfrac{3}{2}}\, \frac{\sqrt{\sum_i (\lambda_i - \overline{\lambda})^2}}{\sqrt{\sum_i \lambda_i^2}},\]
with \(\text{FA}\in[0,1]\): 0 for isotropic diffusion (e.g. free water, \(\lambda_1 =\lambda_2=\lambda_3\)) and approaching 1 for a perfectly linear diffusion profile. This is the literal sense in which medical imaging takes us “from biological structure and function to images and tensors.”
## a single voxel's diffusion tensor (units mm^2/s), as fitted from DWI data
D <- matrix(c(1.70, 0.20, 0.10,
0.20, 0.50, 0.05,
0.10, 0.05, 0.40), 3, 3) * 1e-3
eig <- eigen(D, symmetric = TRUE)
lam <- eig$values
MD <- mean(lam)
FA <- sqrt(1.5 * sum((lam - MD)^2) / sum(lam^2))
data.frame(lambda1 = lam[1], lambda2 = lam[2], lambda3 = lam[3],
MD_mm2_per_s = MD, FA = FA)## lambda1 lambda2 lambda3 MD_mm2_per_s FA
## 1 0.001741 0.00048 0.000379 0.0008667 0.7123## [1] 0.984 0.162 0.079Important consequences of the volumetric and temporal (3D/4D) view:
- filtering and segmentation can use 3D neighborhoods rather than only 2D ones;
- registration may be required across time points, not just across subjects;
- temporal filtering and time-series modeling become available for 4D data;
- voxel-wise time courses can be extracted from a 4D dataset;
- visualization may require orthogonal slice views, maximum-intensity projections, surface rendering, or volume rendering.
Take-home. Every image-processing method in this chapter ultimately operates on the numbers stored at pixel or voxel locations of a tensor. The image looks visual, but mathematically it is an indexed array, and the operations we apply are tensor operations.
Checkpoint 1. In the forward model \(g=\mathcal{S}[b\cdot(h*f)]+\eta\), which term explains each of the following? (a) a bright rim on the surface coil side of a knee MRI; (b) a 4 mm lesion appearing 6 mm wide; (c) grainy texture that changes between two back-to-back scans; (d) a stair-step appearance of a vessel running obliquely through the slice. Answers: (a) \(b\); (b) \(h\); (c) \(\eta\); (d) \(\mathcal{S}\).
The general image-processing workflow
Medical image processing converts acquired images into reliable quantitative data. A general workflow can be summarized as a pipeline:
\[\text{Import} \rightarrow \text{Inspect} \rightarrow \text{Preprocess} \rightarrow \text{Filter} \rightarrow \text{Register} \rightarrow \text{Segment} \rightarrow \text{Measure} \rightarrow \text{Model} \rightarrow \text{Validate}\]
Not every study uses every step, and the order can vary, but most pipelines contain some version of the following:
- Data import and representation, load image arrays together with metadata: pixel spacing, voxel size, modality, acquisition parameters, orientation, and intensity calibration.
- Inspection and visualization, display slices, projections, color maps, and histograms to check quality, contrast, artifacts, and intensity distributions.
- Preprocessing, normalize intensities, correct artifacts and shading, reduce background variation, and prepare the image for analysis.
- Denoising and enhancement, improve visibility of relevant structures while preserving biologically meaningful boundaries.
- Registration, align images across time points, subjects, or modalities.
- Thresholding and segmentation, identify structures of interest.
- Feature extraction and radiomics, convert regions or whole images into quantitative measurements.
- Modeling and AI, use features or whole images for classification, prediction, diagnosis, prognosis, or treatment planning.
- Validation and interpretation, evaluate whether results are reproducible, biologically plausible, and clinically meaningful.
Workflow principle. Each step affects the final quantitative result, and the errors compose. If step \(k\) contributes an independent relative uncertainty \(u_k\), the pipeline’s relative uncertainty is at least \(\sqrt{\sum_k u_k^2}\), and considerably more if the steps are correlated (as smoothing and thresholding are). The final measurement is only as reliable as the weakest step in the pipeline.
A recurring distinction: display vs. quantitative operations. Windowing, color maps, gamma correction, and interpolation-for-display can make an image easier to read without preserving quantitative intensity values. Normalization, denoising, registration, segmentation, and feature extraction change the data used for measurement. Quantitative analysis should be performed on well-documented, appropriately corrected data, not on display-optimized images.
Image import and representation
The first step is to load the image into the analysis environment. Import may involve reading DICOM, NIfTI, TIFF, PNG, JPEG, or proprietary formats; preserving metadata; checking dimensions and bit depth; confirming orientation; and loading a single image, a stack, or a time series.
Common pitfall. Exporting a medical image as a screenshot or lossy JPEG can destroy quantitative information (it rescales intensities, applies a display window, and may add compression artifacts). Quantitative analysis must start from the original image format, with the original intensity calibration intact.
A reproducible synthetic phantom built from the forward model
Throughout this chapter we use a synthetic phantom constructed explicitly from the forward model of the previous section, so that every degradation we later try to undo has a known cause and a known magnitude:
\[g \;=\; \underbrace{b}_{\text{bias field}}\cdot \big(\underbrace{h}_{\text{Gaussian PSF}} * \underbrace{f}_{\text{piecewise-constant truth}}\big) \;+\; \underbrace{\eta}_{\mathcal{N}(0,\sigma^2)} .\]
Using a synthetic image makes every example fully reproducible without external files and, crucially, lets us compare measurements against a known ground truth.
Geometric convention. We adopt the display convention used by every figure in this chapter: the first matrix index is the row and increases downward on screen; the second index is the column and increases to the right. The continuous coordinates \((x,y)\) are attached accordingly, with \(y\) increasing upward in the usual mathematical sense. Keeping this explicit avoids the single most common bug in image code, a silently transposed or flipped array.
n <- 160 # matrix size (pixels)
FOV <- 160 # field of view (mm)
dx <- FOV / n # in-plane spacing, mm/pixel (columns, x)
dy <- FOV / n # in-plane spacing, mm/pixel (rows, y)
xs <- seq(-1, 1, length.out = n) # x: left -> right (columns)
ys <- seq( 1, -1, length.out = n) # y: top -> bottom (rows)
X <- matrix(rep(xs, each = n), n, n) # X[i,j] = xs[j]
Y <- matrix(rep(ys, times = n), n, n) # Y[i,j] = ys[i]
## rotated-ellipse indicator, evaluated on the coordinate matrices
ellipse <- function(X, Y, x0, y0, a, b, angle = 0) {
xr <- cos(angle) * (X - x0) + sin(angle) * (Y - y0)
yr <- -sin(angle) * (X - x0) + cos(angle) * (Y - y0)
(xr / a)^2 + (yr / b)^2 <= 1
}
## ---- ground-truth compartments -------------------------------------------
brain_mask <- ellipse(X, Y, 0.00, 0.00, 0.82, 0.95)
tumor_mask <- ellipse(X, Y, 0.34, -0.34, 0.20, 0.14, angle = 0.5)
vessel_mask <- abs(Y - 0.35 * sin(5 * X)) < 0.022 & abs(X) < 0.70 & brain_mask
lesion_mask <- ellipse(X, Y, -0.42, 0.34, 0.06, 0.06)
par_mask <- brain_mask & !tumor_mask & !vessel_mask & !lesion_mask
## ---- f : the true, noise-free, infinitely sharp object -------------------
img_true <- matrix(0.05, n, n) # background (air / holder)
img_true[brain_mask] <- 0.40 # parenchyma
img_true[vessel_mask] <- 0.62 # vessel
img_true[lesion_mask] <- 0.55 # small lesion
img_true[tumor_mask] <- 0.80 # tumour (brightest)
## ---- b : multiplicative bias / sensitivity field --------------------------
bias_field <- 1 + 0.18 * X + 0.10 * Y # ranges over [0.72, 1.28]
## ---- h : Gaussian PSF, sigma in pixels -----------------------------------
psf_sigma <- 1.2 # FWHM = 2.355 * 1.2 = 2.83 px = 2.83 mm
sigma_noise <- 0.04 # additive Gaussian noise SD## --- padding, used by every neighborhood operator in this chapter ---------
pad_image <- function(img, pr, pc, mode = c("replicate", "zero", "reflect")) {
mode <- match.arg(mode); nr <- nrow(img); nc <- ncol(img)
if (mode == "zero") {
out <- matrix(0, nr + 2 * pr, nc + 2 * pc)
out[(pr + 1):(pr + nr), (pc + 1):(pc + nc)] <- img
return(out)
}
ri <- if (mode == "replicate") c(rep(1, pr), 1:nr, rep(nr, pr))
else c(if (pr) (pr + 1):2 else integer(0), 1:nr,
if (pr) (nr - 1):(nr - pr) else integer(0))
ci <- if (mode == "replicate") c(rep(1, pc), 1:nc, rep(nc, pc))
else c(if (pc) (pc + 1):2 else integer(0), 1:nc,
if (pc) (nc - 1):(nc - pc) else integer(0))
img[ri, ci, drop = FALSE]
}
## --- 2D convolution, vectorised over kernel entries (fast and exact) -------
conv2 <- function(img, kernel, mode = "replicate") {
kr <- nrow(kernel); kc <- ncol(kernel)
pr <- (kr - 1) %/% 2; pc <- (kc - 1) %/% 2
P <- pad_image(img, pr, pc, mode)
nr <- nrow(img); nc <- ncol(img)
out <- matrix(0, nr, nc)
for (u in 1:kr) for (v in 1:kc) {
w <- kernel[kr - u + 1L, kc - v + 1L] # flip the kernel: true convolution
if (w != 0) out <- out + w * P[u:(u + nr - 1), v:(v + nc - 1)]
}
out
}
## --- separable Gaussian blur (O(k) instead of O(k^2) per pixel) ------------
gauss_blur <- function(img, sigma, mode = "replicate") {
if (sigma <= 0) return(img)
h <- ceiling(3 * sigma); ax <- -h:h
g <- exp(-ax^2 / (2 * sigma^2)); g <- g / sum(g)
conv2(conv2(img, matrix(g, nrow = 1), mode), matrix(g, ncol = 1), mode)
}
normalize01 <- function(x) { r <- range(x, na.rm = TRUE); (x - r[1]) / diff(r) }
## --- display helper: fixed orientation, optional fixed intensity window ----
## zlim = NULL -> auto-scale each panel (convenient, but NOT comparable)
## zlim = c(a,b) -> identical display window for every panel (comparable)
show_img <- function(m, main = "", col = gray.colors(256), zlim = NULL) {
if (is.logical(m)) m <- m * 1
if (is.null(zlim)) zlim <- range(m, finite = TRUE)
m <- pmin(pmax(m, zlim[1]), zlim[2])
image(x = 1:ncol(m), y = 1:nrow(m), z = t(m[nrow(m):1, , drop = FALSE]),
col = col, zlim = zlim, axes = FALSE, xlab = "", ylab = "",
main = main, asp = 1, useRaster = TRUE)
}img_blur <- gauss_blur(img_true, psf_sigma) # h * f
img_noisy <- bias_field * img_blur +
matrix(rnorm(n * n, 0, sigma_noise), n, n) # b(h*f) + eta
par(mfrow = c(1, 3), mar = c(1, 1, 2.5, 1))
show_img(img_true, "f: true object", zlim = c(0, 1))
show_img(img_blur, "h*f: after the PSF", zlim = c(0, 1))
show_img(img_noisy, "g: observed (bias + noise)", zlim = c(0, 1))
Note that all three panels share the same display window, zlim = c(0,1), so
they are visually comparable. Auto-scaling each panel separately, the default
in most plotting functions, hides exactly the intensity changes we are trying
to study.
data.frame(
compartment = c("background", "parenchyma", "vessel", "small lesion", "tumour"),
true_value = c(0.05, 0.40, 0.62, 0.55, 0.80),
pixels = c(sum(!brain_mask), sum(par_mask), sum(vessel_mask),
sum(lesion_mask), sum(tumor_mask)),
area_mm2 = c(sum(!brain_mask), sum(par_mask), sum(vessel_mask),
sum(lesion_mask), sum(tumor_mask)) * dx * dy
)## compartment true_value pixels area_mm2
## 1 background 0.05 10144 10144
## 2 parenchyma 0.40 14435 14435
## 3 vessel 0.62 396 396
## 4 small lesion 0.55 72 72
## 5 tumour 0.80 553 553The clean image img_true is the signal we wish we could measure; img_noisy
is the observation we actually get. Much of image processing is the attempt to
recover information about the former from the latter.
Images as numerical arrays
A grayscale image is a matrix \(I(x,y)\), where \(x\) and \(y\) are spatial coordinates and \(I\) is the measured intensity. A small image is just a table of numbers:
small_img <- matrix(c(
12, 15, 17, 18,
10, 14, 20, 22,
11, 16, 25, 30,
9, 13, 19, 24
), nrow = 4, byrow = TRUE)
small_img## [,1] [,2] [,3] [,4]
## [1,] 12 15 17 18
## [2,] 10 14 20 22
## [3,] 11 16 25 30
## [4,] 9 13 19 24(We deliberately name this small_img so it does not overwrite the
\(160\times160\) phantom used elsewhere.) A 3D volume is \(I(x,y,z)\) and a dynamic
series is \(I(x,y,z,t)\); the same indexing logic extends naturally.
Pixels and voxels
A pixel is the smallest addressable element of a 2D image; it has a location and an intensity value. A voxel is its 3D analogue. For a 3D image, the physical size of a voxel is set by the field of view (FOV) and the matrix size:
\[\Delta x = \frac{\text{FOV}_x}{N_x}, \qquad \Delta y = \frac{\text{FOV}_y}{N_y}, \qquad \Delta z = \text{slice thickness (or through-plane spacing)},\]
so the voxel volume is \(V_{\text{voxel}} = \Delta x \,\Delta y\, \Delta z\).
Physics: the resolution–noise–dose triangle. In every quantum-limited modality the number of detected events in a voxel scales with the voxel volume and the acquisition time, \(N \propto \Delta x\,\Delta y\,\Delta z\, T\). For Poisson counting statistics \(\sigma_N = \sqrt{N}\), so
\[\text{SNR} \;=\; \frac{N}{\sqrt N} \;=\; \sqrt{N} \;\propto\; \sqrt{\Delta x\,\Delta y\,\Delta z\; T}\,.\]
Three consequences follow immediately and are worth memorizing:
- Halving the linear voxel size in all three directions divides the voxel volume by 8 and the SNR by \(\sqrt{8}\approx 2.83\).
- Halving the noise at fixed resolution requires four times the counts, hence four times the dose (CT, PET) or four times the scan time (MRI averages, \(\text{SNR}\propto\sqrt{N_{\text{avg}}}\)).
- Therefore resolution is never “free”: it is purchased with dose, time, or noise. Choosing a voxel size is a clinical trade-off decision.
## Poisson counting: SNR = sqrt(N) with N proportional to voxel volume
vox_mm <- c(0.5, 0.7, 1.0, 1.5, 2.0, 3.0)
counts <- 2000 * (vox_mm / 1.0)^3 # 2000 counts at 1 mm isotropic
snr <- sqrt(counts)
plot(vox_mm, snr, type = "b", pch = 19, log = "xy", col = "darkred",
xlab = "isotropic voxel size (mm)", ylab = "SNR (Poisson limit)",
main = expression(SNR %prop% sqrt(Delta*x*Delta*y*Delta*z)))
grid()
data.frame(voxel_mm = vox_mm, counts = round(counts), SNR = round(snr, 1),
dose_multiplier_for_SNR20 = round((20 / snr)^2, 2))## voxel_mm counts SNR dose_multiplier_for_SNR20
## 1 0.5 250 15.8 1.60
## 2 0.7 686 26.2 0.58
## 3 1.0 2000 44.7 0.20
## 4 1.5 6750 82.2 0.06
## 5 2.0 16000 126.5 0.03
## 6 3.0 54000 232.4 0.01Clinical example. A tumor-volume measurement depends on voxel size. Counting 1,000 voxels is meaningless unless we know the physical volume of each voxel — and comparing 1,000 voxels on a 1 mm scan against 1,000 voxels on a 3 mm scan is comparing 1 cm\(^3\) against 27 cm\(^3\).
Image geometry: from voxel indices to patient coordinates
Pixel spacing alone is not enough to place a voxel in the patient. Medical image formats store a full affine geometry: a \(4\times4\) matrix mapping voxel indices \((i,j,k)\) to physical coordinates \((x,y,z)\) in millimeters,
\[\begin{pmatrix} x\\ y\\ z\\ 1\end{pmatrix} = \underbrace{\begin{pmatrix} R\,S & t \\ 0 & 1\end{pmatrix}}_{A} \begin{pmatrix} i\\ j\\ k\\ 1\end{pmatrix}, \qquad S=\operatorname{diag}(\Delta x,\Delta y,\Delta z),\]
where the columns of \(R\) are the direction cosines of the image axes (DICOM tag Image Orientation (Patient)) and \(t\) is the position of voxel \((0,0,0)\) (Image Position (Patient)). Two facts make this practically essential:
- the physical voxel volume is \(|\det(RS)| = \Delta x\Delta y\Delta z\) only when \(R\) is orthonormal, for oblique or sheared acquisitions you must use the determinant, not the product of spacings;
- \(\det(A) < 0\) signals a left–right flip between the stored array and the physical frame. DICOM uses LPS (+x = patient left, +y = posterior, +z = superior); NIfTI conventionally uses RAS. Mixing them silently mirrors the patient, an error that has reached the operating room.
## An oblique axial acquisition: 0.7 x 0.7 mm in-plane, 3 mm slices,
## rotated 10 degrees about the z axis, first voxel at (-90, -110, 40) mm.
theta <- 10 * pi / 180
Rdir <- matrix(c(cos(theta), sin(theta), 0,
-sin(theta), cos(theta), 0,
0, 0, 1), 3, 3)
S <- diag(c(0.7, 0.7, 3.0))
tvec <- c(-90, -110, 40)
A <- rbind(cbind(Rdir %*% S, tvec), c(0, 0, 0, 1))
dimnames(A) <- list(c("x","y","z",""), c("i","j","k","1"))
round(A, 3)## i j k 1
## x 0.689 -0.122 0 -90
## y 0.122 0.689 0 -110
## z 0.000 0.000 3 40
## 0.000 0.000 0 1voxel_to_mm <- function(A, i, j, k) as.vector(A %*% c(i, j, k, 1))[1:3]
round(voxel_to_mm(A, 128, 128, 10), 2) # physical position of a voxel## [1] -17.32 -6.20 70.00## voxel volume: use the determinant, never the naive product, for oblique data
det_vol <- abs(det(Rdir %*% S))
naive_vol <- prod(diag(S))
c(det_volume_mm3 = det_vol, naive_product_mm3 = naive_vol,
left_right_flip = det(A) < 0)## det_volume_mm3 naive_product_mm3 left_right_flip
## 1.47 1.47 0.00Common pitfall. Resampling an image without carrying its affine forward, or writing a mask with the identity affine, decouples the labels from the anatomy. Masks must always be stored with the same geometry as the image they were drawn on, and any tool that reports “volume = voxel count \(\times\) spacing” should be checked against \(\left|\det(RS)\right|\).
Metadata and intensity calibration
Metadata are data about the image. Typical fields include modality; patient/study identifiers; acquisition date and time; image dimensions; pixel spacing; slice thickness; orientation; scanner settings; sequence parameters; reconstruction parameters; and calibration information, carried in DICOM headers and, for research volumes, in NIfTI headers.
Raw stored values are almost never the physical quantity. DICOM stores a linear map from stored value \(SV\) to output units,
\[\text{value} = m\cdot SV + c,\]
with \(m\) = Rescale Slope (0028,1053) and \(c\) = Rescale Intercept (0028,1052). For CT the output unit is the Hounsfield unit, defined by water and air:
\[\text{HU} = 1000\,\frac{\mu - \mu_{\text{water}}}{\mu_{\text{water}} - \mu_{\text{air}}} \;\approx\; 1000\,\frac{\mu - \mu_{\text{water}}}{\mu_{\text{water}}},\]
so water is 0 HU and air is \(-1000\) HU by construction. For PET the standard semi-quantitative unit is the standardized uptake value
\[\text{SUV}_{\text{bw}} = \frac{C_{\text{tissue}}(t)\ [\text{kBq/mL}]} {A_{\text{inj}}\,e^{-\lambda t}\ [\text{kBq}] \,/\, W\ [\text{g}]},\]
with \(A_{\text{inj}}\) the injected activity, \(\lambda=\ln 2/T_{1/2}\) the decay constant, and \(W\) the body weight; SUV \(\approx 1\) means uptake equal to a uniform whole-body distribution.
## --- CT: linear attenuation coefficients at ~70 keV -> Hounsfield units ----
mu_water <- 0.0193 # 1/mm
mu <- c(air = 0.0000, lung = 0.0002, fat = 0.0177, water = 0.0193,
soft_tissue = 0.0208, trabecular_bone = 0.0270, cortical_bone = 0.0480)
HU <- 1000 * (mu - mu_water) / mu_water
round(HU)## air lung fat water soft_tissue
## -1000 -990 -83 0 78
## trabecular_bone cortical_bone
## 399 1487## --- DICOM rescale: stored values are unsigned; HU needs slope/intercept ---
stored_values <- c(0L, 24L, 1024L, 1100L, 2500L)
slope <- 1; intercept <- -1024
data.frame(stored = stored_values, HU = slope * stored_values + intercept)## stored HU
## 1 0 -1024
## 2 24 -1000
## 3 1024 0
## 4 1100 76
## 5 2500 1476## --- PET: decay-corrected SUV ---------------------------------------------
suv <- function(conc_kBq_per_mL, dose_MBq, weight_kg, minutes_post_inj,
half_life_min = 109.77) { # F-18
dose_kBq_dec <- dose_MBq * 1000 * exp(-log(2) * minutes_post_inj / half_life_min)
conc_kBq_per_mL / (dose_kBq_dec / (weight_kg * 1000))
}
round(suv(conc_kBq_per_mL = 18.5, dose_MBq = 370,
weight_kg = 75, minutes_post_inj = 60), 2)## [1] 5.48Important caution. Intensity values do not mean the same thing across modalities. CT intensity is calibrated in Hounsfield units and is comparable across scanners; MRI intensity is usually relative and depends strongly on sequence, coil, and vendor scaling; SUV is calibrated but sensitive to uptake time, blood glucose, and partial-volume effects. The mathematical operations in this chapter are general, but their biological interpretation is modality-specific.
| Modality | Meaning of intensity | Calibrated? |
|---|---|---|
| CT | Hounsfield units, linear in attenuation coefficient | Yes (water/air anchored) |
| MRI | signal from sequence, PD, \(T_1\), \(T_2\), \(T_2^*\), diffusion, flow | Usually no (arbitrary units) |
| Quantitative MRI (\(T_1\), \(T_2\), ADC maps) | physical relaxation/diffusion parameter | Yes |
| Ultrasound (B-mode) | log-compressed echo amplitude, gain-dependent | No |
| PET / SPECT | radiotracer activity concentration, often as SUV | Yes (with corrections) |
| Fluorescence microscopy | photons from fluorophore or autofluorescence | Rarely |
Point spread function, resolution, and blur
The point spread function (PSF) describes how an imaging system represents an ideal point object. Instead of a perfect point, the system produces a blurred spot whose width sets the spatial resolution. For a linear, shift-invariant (LSI) system,
\[I_{\text{obs}}(x,y) = (h * I_{\text{true}})(x,y) + \eta(x,y).\]
Resolution is quoted as the full width at half maximum (FWHM) of the PSF. For a Gaussian PSF of standard deviation \(\sigma_h\),
\[\text{FWHM} = 2\sqrt{2\ln 2}\;\sigma_h \approx 2.3548\,\sigma_h .\]
Two independent blurs combine (for Gaussians, exactly; otherwise approximately) in quadrature:
\[\text{FWHM}_{\text{total}}^2 = \text{FWHM}_1^2 + \text{FWHM}_2^2 .\]
This is why a 4 mm intrinsic detector resolution and a 4 mm reconstruction filter give \(\sqrt{32}=5.7\) mm, not 8 mm, and why post-smoothing a good scanner quickly throws away the resolution it paid for.
fwhm_of_sigma <- function(s) 2 * sqrt(2 * log(2)) * s
c(psf_sigma_px = psf_sigma,
psf_FWHM_px = fwhm_of_sigma(psf_sigma),
psf_FWHM_mm = fwhm_of_sigma(psf_sigma) * dx)## psf_sigma_px psf_FWHM_px psf_FWHM_mm
## 1.200 2.826 2.826## quadrature addition, verified numerically on our phantom
s1 <- 1.2; s2 <- 1.6
lhs <- gauss_blur(gauss_blur(img_true, s1), s2)
rhs <- gauss_blur(img_true, sqrt(s1^2 + s2^2))
c(max_abs_difference = max(abs(lhs - rhs)),
FWHM_combined_mm = fwhm_of_sigma(sqrt(s1^2 + s2^2)) * dx)## max_abs_difference FWHM_combined_mm
## 0.0002972 4.7096401The PSF links image formation, spatial resolution, blur, convolution, and the frequency-domain analysis developed later, where its normalized Fourier magnitude, the modulation transfer function (MTF), quantifies how much contrast survives at each spatial frequency.
Partial-volume effects and the recovery coefficient
A voxel may contain a mixture of tissues; its recorded intensity is then a weighted average rather than a pure-tissue value. This partial-volume effect (PVE) is strongest when the object is small relative to the PSF, when a boundary passes through the voxel, when structures are thin, curved, or oblique relative to the grid, and after resampling or smoothing.
PVE is quantified by the contrast recovery coefficient
\[\text{CRC} = \frac{\bar I_{\text{measured}} - I_{\text{background}}} {I_{\text{true}} - I_{\text{background}}},\]
which approaches 1 for large objects and falls steeply once the object diameter drops below about \(2\times\)FWHM. This single curve explains why small lesions look less bright than they are (PET SUV underestimation), why small structures are systematically under-measured in intensity and over-measured in size, and why lesion size must always accompany a reported concentration.
## disks of increasing radius, blurred by the same PSF, measured in the true ROI
rc_tab <- t(sapply(c(1.5, 2, 3, 4, 6, 8, 12, 20), function(R) {
S <- 121; cc <- 61
d <- (row(matrix(0, S, S)) - cc)^2 + (col(matrix(0, S, S)) - cc)^2 <= R^2
b <- gauss_blur(d * 1, psf_sigma)
c(radius_px = R,
diam_over_FWHM = 2 * R / fwhm_of_sigma(psf_sigma),
peak_CRC = max(b),
mean_CRC = mean(b[d]))
}))
plot(rc_tab[, "diam_over_FWHM"], rc_tab[, "mean_CRC"], type = "b", pch = 19,
ylim = c(0, 1.05), col = "darkorange3",
xlab = "object diameter / PSF FWHM", ylab = "contrast recovery coefficient",
main = "Partial-volume effect: small objects lose contrast")
lines(rc_tab[, "diam_over_FWHM"], rc_tab[, "peak_CRC"], type = "b", pch = 1,
col = "steelblue", lty = 2)
abline(h = 1, lty = 3); abline(v = 2, col = "red", lty = 2)
legend("bottomright", c("mean over true ROI", "peak value"),
col = c("darkorange3", "steelblue"), lty = c(1, 2), pch = c(19, 1), bty = "n")
## radius_px diam_over_FWHM peak_CRC mean_CRC
## [1,] 1.5 1.062 0.644 0.491
## [2,] 2.0 1.416 0.754 0.543
## [3,] 3.0 2.123 0.957 0.684
## [4,] 4.0 2.831 0.996 0.758
## [5,] 6.0 4.247 1.000 0.841
## [6,] 8.0 5.662 1.000 0.879
## [7,] 12.0 8.493 1.000 0.920
## [8,] 20.0 14.155 1.000 0.953## the same effect inside our phantom: measured vs true contrast
bg <- 0.40 # parenchyma
crc <- function(mask, true_val)
(mean(img_blur[mask]) - bg) / (true_val - bg)
data.frame(
object = c("tumour", "small lesion", "vessel"),
true_value = c(0.80, 0.55, 0.62),
measured_mean= round(c(mean(img_blur[tumor_mask]), mean(img_blur[lesion_mask]),
mean(img_blur[vessel_mask])), 3),
CRC = round(c(crc(tumor_mask, 0.80), crc(lesion_mask, 0.55),
crc(vessel_mask, 0.62)), 3)
)## object true_value measured_mean CRC
## 1 tumour 0.80 0.771 0.927
## 2 small lesion 0.55 0.521 0.804
## 3 vessel 0.62 0.530 0.591From image to measurement. Segmentation supplies the number of pixels or voxels in a structure. Metadata supplies their physical size. The PSF supplies the bias correction. All three are required for a real-world measurement, and only the first two are usually reported.
Checkpoint 2. A PET scanner has 5 mm FWHM resolution. A spherical lesion of 6 mm diameter has a true SUV of 8.0. Roughly what SUV\(_{\text{mean}}\) will be measured, and in which direction is the error? Answer: diameter/FWHM = 1.2, so the CRC is well below 1 (≈0.4–0.5 from the curve above); the measured SUV will be roughly 3–4, a severe underestimate. Reporting it without a lesion-size caveat would be misleading.
Section summary. Metadata describe size, spacing, orientation, acquisition, and calibration; the affine matrix places voxels in the patient; intensity values are modality-specific and often require a rescale map; and the PSF sets both the resolution and the partial-volume bias of every subsequent measurement.
Real images with EBImage (optional)
The base-R pipeline above needs no special packages. For working with real
images we use the Bioconductor package EBImage, which reads, writes,
displays, and processes JPEG/PNG/TIFF images and exposes them as ordinary R
arrays. (If EBImage is unavailable, the jpeg/png/tiff CRAN packages provide
basic reading; RNifti or oro.nifti read NIfTI, and oro.dicom reads DICOM.)
## Sketch of a real import, with the two things that must never be skipped:
## (1) recover physical units, (2) recover geometry.
library(oro.dicom); library(oro.nifti)
dcm <- readDICOM("study/series01")
hdr <- dcm$hdr[[1]]
slope <- as.numeric(extractHeader(dcm$hdr, "RescaleSlope"))[1]
intercept <- as.numeric(extractHeader(dcm$hdr, "RescaleIntercept"))[1]
spacing <- as.numeric(strsplit(extractHeader(dcm$hdr, "PixelSpacing",
numeric = FALSE)[1], " ")[[1]])
vol_HU <- slope * dcm$img[[1]] + intercept # stored value -> HU
nii <- readNIfTI("subject01_T1.nii.gz", reorient = FALSE)
pixdim(nii)[2:4] # voxel spacing (mm)
nii@qform_code; nii@sform_code # which affine is authoritativelibrary(EBImage)
f <- system.file("images", "sample.png", package = "EBImage")
img <- readImage(f)
## The Image class extends base 'array'; pixel data live in the .Data slot.
print(img, short = TRUE)## Image
## colorMode : Grayscale
## storage.mode : double
## dim : 768 512
## frames.total : 1
## frames.render: 1## [1] 768 512## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.4471 0.4627 0.4784 0.4980 0.5137 0.5294
## [2,] 0.4510 0.4627 0.4784 0.4824 0.5059 0.5216
## [3,] 0.4627 0.4667 0.4824 0.4980 0.5137 0.5137## [1] 0 1EBImage::display(img, method = "raster")
text(x = 20, y = 20, labels = "Sample image", adj = c(0, 1),
col = "orange", cex = 1.6)
An EBImage Image carries a colorMode slot. For a grayscale image the third
and higher array dimensions index separate frames (z-positions, time points,
replicates); for a Color image the third dimension holds color channels (R, G,
B) and the fourth indexes frames. Multi-frame data are common in microscopy:
## Image
## colorMode : Grayscale
## storage.mode : double
## dim : 510 510 4
## frames.total : 4
## frames.render: 4
##
## imageData(object)[1:5,1:6,1]
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.06275 0.07451 0.07059 0.08235 0.10588 0.09804
## [2,] 0.06275 0.05882 0.07843 0.09020 0.09020 0.10588
## [3,] 0.06667 0.06667 0.08235 0.07843 0.09412 0.09412
## [4,] 0.06667 0.06667 0.07059 0.08627 0.08627 0.09804
## [5,] 0.05882 0.06667 0.07059 0.08235 0.09412 0.10588## [1] 4
Visualization and inspection
Visualization converts numerical arrays into interpretable displays, and inspection is the disciplined habit of looking before measuring. Most catastrophic analysis errors (wrong orientation, corrupted slice, mis-scaled intensities) are obvious on inspection and invisible in a summary statistic.
Key idea. Inspection prevents automated analysis from producing precise but meaningless numbers.
Display operations versus quantitative operations
It bears repeating, because it underlies everything below. Display transformations, windowing/leveling, gamma correction, pseudocolor, display-only interpolation, brightness/contrast, change how an image looks without changing the stored values. Quantitative transformations — normalization, denoising, registration/resampling, segmentation, feature extraction, bias correction, change the data used for measurement. Changing the display window can make a lesion look far more conspicuous without altering a single voxel value; that is appropriate for reading but must never be mistaken for a measurement.
The auto-scaling trap. Almost every plotting function rescales each panel to
its own min/max. Two images that differ by a factor of three in absolute
intensity then look identical. Every comparison figure in this chapter
therefore passes an explicit zlim. If a figure in a paper compares “before”
and “after” without stating the display window, the comparison is
uninterpretable.
half <- img_noisy * 0.5 # a genuinely dimmer image
par(mfrow = c(1, 3), mar = c(1, 1, 2.5, 1))
show_img(img_noisy, "Original (auto-scaled)")
show_img(half, "Half intensity (auto-scaled)")
show_img(half, "Half intensity (zlim = 0..1)", zlim = c(0, 1))
## mean_original mean_half
## 0.2741 0.1370Grayscale, pseudocolor, and lookup tables
A grayscale lookup table maps intensities to shades of gray; a pseudocolor map assigns colors to intensities to emphasize small differences. Pseudocolor can reveal subtle gradients but can also mislead, because the human visual system perceives some color transitions as sharper than the underlying data. A perceptually uniform map (viridis, cividis, magma) keeps perceived differences proportional to data differences; classic rainbow/jet maps do not, and are known to manufacture apparent edges at the cyan and yellow bands.
ramp <- matrix(rep(seq(0, 1, length.out = n), each = 40), 40, n)
par(mfrow = c(1, 3), mar = c(1, 1, 2.5, 1))
show_img(img_noisy, "Grayscale", zlim = c(0, 1))
show_img(img_noisy, "Pseudocolor (heat)", zlim = c(0, 1), col = heat.colors(256))
show_img(img_noisy, "Pseudocolor (rainbow)", zlim = c(0, 1), col = rainbow(256, end = 0.7))
par(mfrow = c(1, 1))
par(mfrow = c(3, 1), mar = c(1, 1, 2, 1))
show_img(ramp, "Linear ramp in grayscale (uniform)", zlim = c(0, 1))
show_img(ramp, "Linear ramp in heat", zlim = c(0, 1), col = heat.colors(256))
show_img(ramp, "Linear ramp in rainbow (banding!)", zlim = c(0, 1), col = rainbow(256, end = 0.7))
Teaching prompt. The three ramps encode exactly the same linear data. Where does the rainbow version appear to have edges? Those apparent boundaries are artifacts of the colormap, not of the data, which is precisely the failure mode when a rainbow map is used on a perfusion or ADC map.
Windowing and contrast adjustment
Windowing (window/level) maps a chosen intensity interval \([L - W/2,\; L + W/2]\), centered at level \(L\) with width \(W\), onto the full display range, clipping values outside the window:
\[I_{\text{disp}} = \operatorname{clip}\!\left(\frac{I - (L - W/2)}{W},\,0,\,1\right).\]
It is the standard way radiologists inspect CT, where different windows reveal different structures from the same underlying Hounsfield data.
window_level <- function(img, level, width)
pmin(pmax((img - (level - width / 2)) / width, 0), 1)
## a synthetic CT slice in Hounsfield units built from our phantom compartments
ct_hu <- matrix(-1000, n, n) # air
ct_hu[brain_mask] <- 40 # soft tissue / brain
ct_hu[vessel_mask] <- 260 # contrast-filled vessel
ct_hu[lesion_mask] <- 70
ct_hu[tumor_mask] <- 900 # calcified / dense lesion
ct_hu <- gauss_blur(ct_hu, psf_sigma) + matrix(rnorm(n * n, 0, 12), n, n)
par(mfrow = c(1, 3), mar = c(1, 1, 2.5, 1))
show_img(window_level(ct_hu, level = 40, width = 400), "Soft-tissue window (L 40 / W 400)", zlim = c(0, 1))
show_img(window_level(ct_hu, level = 300, width = 1500), "Bone window (L 300 / W 1500)", zlim = c(0, 1))
show_img(window_level(ct_hu, level = -600, width = 1500), "Lung window (L -600 / W 1500)", zlim = c(0, 1))
All three panels display identical data. The window changes only which range of Hounsfield units is mapped onto the visible gray levels, the stored HU values and any measurement made from them are unchanged.
Gamma correction
Gamma correction is a nonlinear display transform, \(I_{\text{out}} = I_{\text{in}}^{\gamma}\) (on a \([0,1]\) scale). For \(\gamma < 1\) dark regions brighten; for \(\gamma > 1\) bright structures are emphasized.
img01 <- normalize01(img_noisy)
par(mfrow = c(1, 3), mar = c(1, 1, 2.5, 1))
show_img(img01, "Original", zlim = c(0, 1))
show_img(img01^0.5, expression(gamma == 0.5), zlim = c(0, 1))
show_img(img01^2.0, expression(gamma == 2), zlim = c(0, 1))
Gamma correction changes display contrast but not the underlying biological signal; display enhancement must not be confused with quantitative measurement. Note that a gamma transform is monotone, so it preserves the rank order of voxels: a threshold applied after gamma is equivalent to a different threshold applied before it. Any non-monotone display transform, by contrast, destroys information irreversibly.
Histograms and intensity distributions
The intensity histogram summarizes how often each value occurs and is the single most useful inspection tool: it reveals dynamic range, clipping, multimodality (which motivates thresholding), and bias-field skew.
par(mfrow = c(1, 2))
hist(as.vector(img_noisy), breaks = 60, col = "gray80", border = "gray50",
main = "Whole-image histogram", xlab = "Image intensity")
abline(v = c(0.05, 0.40, 0.55, 0.62, 0.80), col = "red", lty = 2)
hist(img_noisy[brain_mask], breaks = 60, col = "steelblue", border = "gray50",
main = "Histogram inside the head mask", xlab = "Image intensity")
abline(v = c(0.40, 0.55, 0.62, 0.80), col = "red", lty = 2)
The whole-image histogram is dominated by background; the masked histogram reveals the tissue classes. Always histogram inside a mask when the background occupies most of the field of view.
Noise statistics: what kind of noise is it?
The statistical character of noise is set by the physics of detection, and choosing the right filter or estimator depends on getting this right.
Physics of the four common noise models.
- Poisson (shot) noise, photon or event counting in CT, PET, SPECT, and photon-limited microscopy. \(\operatorname{Var}(N)=\mathbb{E}[N]\), so noise is signal-dependent and \(\text{SNR}=\sqrt{N}\). The Anscombe transform \(A(N)=2\sqrt{N+3/8}\) maps it to approximately unit-variance Gaussian noise, so that ordinary Gaussian denoisers become valid.
- Gaussian noise, a good approximation for high-count regimes and for thermal/electronic noise; also the limit of Poisson as \(N\) grows.
- Rician noise, magnitude MRI. If the complex signal has amplitude \(A\) and each channel carries Gaussian noise of SD \(\sigma\), the magnitude \(m\) follows \(p(m)=\frac{m}{\sigma^2}\exp\!\big(-\frac{m^2+A^2}{2\sigma^2}\big) I_0\!\big(\frac{mA}{\sigma^2}\big)\). Two consequences matter: in signal-free regions the magnitude is Rayleigh with mean \(\sigma\sqrt{\pi/2}\approx1.253\sigma\) and SD \(\approx 0.655\sigma\) (hence the classic correction \(\hat\sigma = \text{SD}_{\text{background}}/0.655\)), and low-SNR voxels are biased upward, a noise floor that inflates ADC and \(T_2\) estimates.
- Speckle, coherent imaging (ultrasound, OCT). Fully developed speckle is multiplicative, with SNR \(\approx 1.91\) independent of gain; it is texture, not additive noise, and is better handled in the log domain.
set.seed(11)
lam <- 25
pois <- rpois(2e5, lam)
ansc <- 2 * sqrt(pois + 3/8)
sig <- 1
ric <- sqrt((0 + rnorm(2e5, 0, sig))^2 + rnorm(2e5, 0, sig)^2) # A = 0
spk <- rexp(2e5, 1) # intensity speckle
par(mfrow = c(2, 2))
hist(pois, breaks = 40, col = "gray80", main = "Poisson (lambda = 25)", xlab = "counts")
hist(ansc, breaks = 40, col = "gray80",
main = sprintf("Anscombe transformed (SD = %.2f)", sd(ansc)), xlab = "2*sqrt(N+3/8)")
hist(ric, breaks = 40, col = "gray80",
main = sprintf("Rayleigh (signal-free MRI): mean %.2f sigma", mean(ric)/sig),
xlab = "magnitude")
hist(spk, breaks = 40, col = "gray80",
main = sprintf("Speckle intensity: SNR = %.2f", mean(spk)/sd(spk)), xlab = "intensity")
par(mfrow = c(1, 1))
data.frame(
model = c("Poisson", "Anscombe(Poisson)", "Rayleigh (A=0)", "Speckle amplitude"),
mean = round(c(mean(pois), mean(ansc), mean(ric), mean(sqrt(spk))), 3),
sd = round(c(sd(pois), sd(ansc), sd(ric), sd(sqrt(spk))), 3),
theory = c("Var = mean", "Var ~ 1", "mean = 1.253 sigma, sd = 0.655 sigma",
"SNR ~ 1.91 (intensity)")
)## model mean sd theory
## 1 Poisson 25.007 4.998 Var = mean
## 2 Anscombe(Poisson) 10.026 0.999 Var ~ 1
## 3 Rayleigh (A=0) 1.254 0.654 mean = 1.253 sigma, sd = 0.655 sigma
## 4 Speckle amplitude 0.885 0.464 SNR ~ 1.91 (intensity)Rician noise floor. Because \(\mathbb{E}[m] \ge A\) always, magnitude MRI never averages to zero in background. Fitting a mono-exponential decay \(S(TE)=S_0e^{-TE/T_2}\) through low-SNR points therefore overestimates \(T_2\); the standard remedies are to fit with a Rician likelihood, subtract the noise floor in quadrature, or truncate the fit above an SNR threshold.
Image quality metrics
Before applying complex processing, ask whether the image has enough quality to support the intended analysis.
- Signal-to-noise ratio (SNR), how strongly signal rises above noise.
- Contrast-to-noise ratio (CNR), how well two regions can be distinguished relative to noise.
- Spatial resolution, smallest separation between distinguishable structures (see PSF/MTF).
- Temporal resolution, time between repeated measurements in 4D imaging.
- Dynamic range / bit depth, range of representable intensity values.
- Partial-volume effects, mixing of tissues within a voxel.
The two workhorse definitions are
\[\text{SNR} = \frac{\mu_{\text{signal}}}{\sigma_{\text{noise}}}, \qquad \text{CNR}_{AB} = \frac{|\mu_A - \mu_B|}{\sigma_{\text{noise}}} .\]
The Rose criterion. For a uniform disk of diameter \(d\) against a background, detectability by a human observer requires roughly
\[\text{SNR}_{\text{Rose}} \;=\; \text{CNR}\cdot\sqrt{\tfrac{\pi}{4}}\,\frac{d}{\Delta x} \;\gtrsim\; 3\text{–}5 .\]
Because the object’s area enters, a lesion twice as wide is detectable at half the contrast. This is the quantitative reason why low-contrast detectability, not high-contrast resolution, drives CT dose protocols.
Estimating the noise level correctly
A naive noise estimate, the standard deviation of a “flat” background region — is biased whenever that region also contains shading, structure, or a reconstruction-dependent noise texture. A robust alternative applies a unit-norm high-pass kernel (which annihilates any locally linear signal) and takes the median absolute deviation of the response:
\[\hat\sigma = \operatorname{MAD}\!\big(k * I\big), \qquad k = \frac{1}{\sqrt{20}}\begin{bmatrix}0&1&0\\1&-4&1\\0&1&0\end{bmatrix}, \quad \lVert k\rVert_2 = 1 .\]
lap_unit <- matrix(c(0, 1, 0, 1, -4, 1, 0, 1, 0), 3, 3, byrow = TRUE) / sqrt(20)
sigma_hat_bg <- sd(img_noisy[!brain_mask]) # naive
sigma_hat_mad <- mad(conv2(img_noisy, lap_unit)[par_mask]) # robust
data.frame(
true_sigma = sigma_noise,
background_SD = round(sigma_hat_bg, 4),
bias_of_background = round(100 * (sigma_hat_bg / sigma_noise - 1), 1),
MAD_highpass = round(sigma_hat_mad, 4),
bias_of_MAD = round(100 * (sigma_hat_mad / sigma_noise - 1), 1)
)## true_sigma background_SD bias_of_background MAD_highpass bias_of_MAD
## 1 0.04 0.0484 21 0.0395 -1.3The background estimate is inflated because the bias field varies across the background; the high-pass MAD estimate is essentially unbiased. This matters: every SNR, CNR, and uncertainty statement downstream inherits the error in \(\hat\sigma\).
mu_tumor <- mean(img_noisy[tumor_mask])
mu_paren <- mean(img_noisy[par_mask])
mu_les <- mean(img_noisy[lesion_mask])
sg <- sigma_hat_mad
## Rose criterion for the two lesions
rose <- function(mu_obj, mu_bg, sigma, diam_px)
abs(mu_obj - mu_bg) / sigma * sqrt(pi / 4) * diam_px
data.frame(
region = c("tumour", "small lesion"),
mean_intensity= round(c(mu_tumor, mu_les), 3),
SNR = round(c(mu_tumor, mu_les) / sg, 1),
CNR_vs_paren = round(abs(c(mu_tumor, mu_les) - mu_paren) / sg, 1),
diameter_px = round(c(2 * sqrt(sum(tumor_mask) / pi), 2 * sqrt(sum(lesion_mask) / pi)), 1),
Rose_SNR = round(rose(c(mu_tumor, mu_les), mu_paren, sg,
c(2 * sqrt(sum(tumor_mask) / pi),
2 * sqrt(sum(lesion_mask) / pi))), 0)
)## region mean_intensity SNR CNR_vs_paren diameter_px Rose_SNR
## 1 tumour 0.792 20.1 10.0 26.5 234
## 2 small lesion 0.504 12.8 2.7 9.6 23Both objects clear the Rose threshold comfortably, which is why they are easy to see, and a useful reminder that a quantitative detectability index can be computed rather than guessed.
Common artifacts
Artifacts are systematic errors or distortions. Unlike random noise, they often have recognizable spatial patterns and, crucially, they do not average out over repeats. Their physical causes are modality-specific, but several categories recur:
| Artifact | Typical physical cause | Where it hurts |
|---|---|---|
| Motion (bulk, respiratory, cardiac) | object moves during encoding | blurring, ghosting, false enhancement |
| Aliasing / wrap-around | sampling below Nyquist | anatomy folded into the FOV |
| Gibbs / truncation ringing | finite k-space or projection sampling | false rims, over/undershoot at edges |
| Beam hardening / streaks | polychromatic X-rays, dense objects, sparse views | HU inaccuracy near bone and metal |
| Intensity non-uniformity | coil sensitivity, illumination, scatter | biased intensity measurements |
| Geometric distortion | \(B_0\) inhomogeneity, gradient nonlinearity | mis-registration, volume error |
| Susceptibility / chemical shift | local field perturbation | signal dropout, spatial displacement |
| Processing artifacts | over-smoothing, interpolation, thresholding | fabricated or erased structure |
A central skill in image analysis is deciding whether an observed pattern is biological or artifactual. Two artifacts deserve explicit demonstration because they are direct consequences of the sampling theory developed later.
Aliasing
Sampling a signal below twice its highest frequency folds that content to a lower apparent frequency. The classic demonstration is a radial chirp, whose local frequency increases with radius.
rr <- (row(matrix(0, n, n)) - n/2)^2 + (col(matrix(0, n, n)) - n/2)^2
chirp <- 0.5 + 0.5 * cos(2 * pi * rr / (2.2 * n))
fac <- 4
sub <- function(m) m[seq(1, nrow(m), by = fac), seq(1, ncol(m), by = fac)]
par(mfrow = c(1, 3), mar = c(1, 1, 2.5, 1))
show_img(chirp, "Original chirp (fully sampled)", zlim = c(0, 1))
show_img(sub(chirp), "Sub-sampled 4x: ALIASED", zlim = c(0, 1))
show_img(sub(gauss_blur(chirp, fac*0.42)),"Anti-alias filtered, then 4x", zlim = c(0, 1))
The middle panel shows spurious low-frequency rings that exist nowhere in the object, they are aliases of frequencies above the new Nyquist limit. Low-pass filtering before decimation (right panel) sacrifices genuine detail but does not fabricate structure.
Medical-imaging caution. Aliasing can fabricate patterns that look anatomical. Always low-pass filter before downsampling, and remember that a “resolution” quoted as matrix size is meaningless without the corresponding anti-alias/reconstruction filter.
Gibbs (truncation) ringing
Truncating k-space, equivalently, convolving with a sinc PSF, produces oscillations near sharp edges whose relative overshoot converges to the Gibbs constant, approximately 8.95 % of the step height, independently of how much data are kept. More data narrows the ringing but does not reduce its amplitude.
trunc1d <- function(v, frac) {
N <- length(v); F <- fft(v)
k <- c(0:(N/2), -(N/2 - 1):-1)
Re(fft(F * (abs(k) <= floor(frac * N / 2)), inverse = TRUE)) / N
}
N1 <- 256; step <- c(rep(0, N1/2), rep(1, N1/2))
ov <- sapply(c(0.10, 0.25, 0.50), function(f) max(trunc1d(step, f)) - 1)
par(mfrow = c(1, 2))
plot(step, type = "l", lwd = 2, col = "gray50", ylim = c(-0.15, 1.2),
xlab = "sample", ylab = "value", main = "Gibbs ringing at a step edge",
xlim = c(100, 156))
for (i in seq_along(c(0.10, 0.25, 0.50)))
lines(trunc1d(step, c(0.10, 0.25, 0.50)[i]), col = c("firebrick","darkorange","steelblue")[i], lwd = 1.5)
abline(h = c(1, 1.0895), lty = c(1, 2), col = c("gray70", "red"))
legend("topleft", c("10% of k-space","25%","50%"),
col = c("firebrick","darkorange","steelblue"), lwd = 1.5, bty = "n", cex = 0.8)
show_img(kspace_trunc <- {
F <- fft(img_true); k1 <- c(0:(n/2), -(n/2-1):-1)
keep <- outer(abs(k1) <= n * 0.12, abs(k1) <= n * 0.12, "&")
Re(fft(F * keep, inverse = TRUE)) / (n * n)
}, "2D: 24% of k-space retained", zlim = c(0, 1))
par(mfrow = c(1, 1))
data.frame(kspace_fraction = c(0.10, 0.25, 0.50),
overshoot = round(ov, 4), Gibbs_constant = 0.0895)## kspace_fraction overshoot Gibbs_constant
## 1 0.10 0.0909 0.0895
## 2 0.25 0.0844 0.0895
## 3 0.50 0.0683 0.0895Clinical relevance. Gibbs ringing produces the classic false “syrinx” line in sagittal cervical-cord MRI and dark/bright rims at fat–water and CSF–cord boundaries. It is not pathology; it is the Fourier transform of a rectangle.
A practical inspection checklist
Inspect before you measure. (1) Does the image load with the expected dimensions and bit depth? (2) Is the geometry right, spacing, slice thickness, orientation, and no left–right flip? (3) Are the intensities in plausible physical units (HU, SUV, ADC)? (4) Is there motion, aliasing, ringing, dropout, or metal artifact? (5) Is the intended structure actually visible, with an adequate Rose SNR? (6) Is the noise level what the protocol predicts? Only after all six pass should automated measurement begin.
Checkpoint 3. A colleague reports “the mean signal in the lesion doubled after treatment” from two conventional \(T_1\)-weighted MRI scans acquired on different days. Name three reasons, all covered so far, why this number may be meaningless. Answer: MRI intensity is uncalibrated (arbitrary units); the receive-coil bias field differs between sessions and positions; and if the lesion shrank below \(\sim 2\times\)FWHM, partial-volume recovery changes the measured mean even when the true concentration is constant.
Preprocessing, normalization, and geometric transformation
Preprocessing prepares the image for segmentation and measurement. It should improve reliability without changing the biological meaning of the image. Common steps include denoising, background/bias correction, artifact reduction, intensity normalization, contrast adjustment, orientation correction, resampling/cropping, and masking.
Important caution. Preprocessing can improve analysis, but it can also create artifacts: over-smoothing erases small lesions, aggressive normalization distorts intensity relationships, and every resampling step blurs. Always document preprocessing choices so that measurements remain reproducible, and apply identical preprocessing to every image that will be compared.
Arithmetic manipulations
Because images are numeric arrays, ordinary arithmetic gives useful transforms: a
negative image (max(img) - img), brightness shifts (addition), contrast
scaling (multiplication), and gamma (exponentiation). Subtraction of two
registered images is the basis of digital subtraction angiography and of
contrast-enhancement maps.
par(mfrow = c(1, 4), mar = c(1, 1, 2.5, 1))
show_img(img01, "Original", zlim = c(0, 1))
show_img(1 - img01, "Negative", zlim = c(0, 1))
show_img(pmin(img01 + 0.2, 1), "Brighter (+0.2)", zlim = c(0, 1))
show_img(pmin(img01 * 1.6, 1), "Higher contrast (x1.6)",zlim = c(0, 1))
## "pre-contrast" and "post-contrast" images differing only in the vessel
pre <- bias_field * gauss_blur(replace(img_true, vessel_mask, 0.40), psf_sigma) +
matrix(rnorm(n * n, 0, sigma_noise), n, n)
post <- img_noisy
par(mfrow = c(1, 3), mar = c(1, 1, 2.5, 1))
show_img(pre, "Pre-contrast", zlim = c(0, 1))
show_img(post, "Post-contrast", zlim = c(0, 1))
show_img(post - pre, "Subtraction: vessel isolated", zlim = c(-0.1, 0.3))
Subtraction cancels everything common to the two images, including the bias field, but doubles the noise variance: if each image has noise SD \(\sigma\), the difference has SD \(\sigma\sqrt2\). Subtraction images are therefore intrinsically noisier and are exquisitely sensitive to misregistration.
Spatial transformations, interpolation, and resampling
Spatial transformations, translation, rotation, reflection, scaling, and the general affine map, reposition image content. In homogeneous coordinates a 2D affine map is
\[\begin{pmatrix} x'\\ y'\\ 1\end{pmatrix} = \begin{pmatrix} a_{11} & a_{12} & t_x\\ a_{21} & a_{22} & t_y\\ 0&0&1 \end{pmatrix}\begin{pmatrix} x\\ y\\ 1\end{pmatrix},\]
with 6 free parameters in 2D (12 in 3D). A rigid map restricts \(A\) to a rotation (\(A^{\top}A=I\), \(\det A = +1\)): 3 parameters in 2D, 6 in 3D. Adding an isotropic scale gives a similarity transform (4 / 7 parameters).
Transformations are applied by inverse mapping: for each output voxel we compute the corresponding input location and interpolate there. (Forward mapping would leave holes.) Common interpolation kernels:
| Method | Kernel | behavior | Correct use |
|---|---|---|---|
| Nearest neighbor | box | no new values introduced | masks, label maps, categorical data |
| Linear / bi- / tri-linear | triangle | mild low-pass, no overshoot | continuous grayscale intensities |
| Cubic B-spline / Lanczos | higher order | sharper, but can overshoot and ring | high-quality resampling of smooth data |
Practical rule. Use nearest-neighbor interpolation for segmentation masks and categorical labels, linear interpolation of the labels 1 and 3 produces the nonexistent label 2. Use linear or spline interpolation for grayscale images. Above all, compose all transformations and resample once: interpolation is a low-pass filter, and repeated resampling compounds its blur irreversibly.
## Bilinear sampling at arbitrary (row, col) source locations -----------------
warp_bilinear <- function(img, i_src, j_src, bg = 0) {
nr <- nrow(img); nc <- ncol(img)
i0 <- floor(i_src); j0 <- floor(j_src)
a <- i_src - i0; b <- j_src - j0
grab <- function(ii, jj) {
ok <- ii >= 1 & ii <= nr & jj >= 1 & jj <= nc
v <- rep(bg, length(ii))
v[ok] <- img[cbind(ii[ok], jj[ok])]
v
}
matrix((1-a)*(1-b)*grab(i0, j0) + (1-a)*b*grab(i0, j0+1) +
a *(1-b)*grab(i0+1, j0) + a *b*grab(i0+1, j0+1), nr, nc)
}
warp_nearest <- function(img, i_src, j_src, bg = 0) {
nr <- nrow(img); nc <- ncol(img)
ii <- round(i_src); jj <- round(j_src)
ok <- ii >= 1 & ii <= nr & jj >= 1 & jj <= nc
v <- rep(bg, length(ii)); v[ok] <- img[cbind(ii[ok], jj[ok])]
matrix(v, nr, nc)
}
## index grids reused by every geometric operation
.row_idx <- function(nr, nc) matrix(rep(1:nr, times = nc), nr, nc)
.col_idx <- function(nr, nc) matrix(rep(1:nc, each = nr), nr, nc)
translate_image <- function(img, dx = 0, dy = 0, bg = 0, interp = "linear") {
nr <- nrow(img); nc <- ncol(img)
f <- if (interp == "nearest") warp_nearest else warp_bilinear
f(img, .row_idx(nr, nc) + dy, .col_idx(nr, nc) - dx, bg) # dx right, dy up
}
rotate_image <- function(img, angle_deg, bg = 0, interp = "linear") {
nr <- nrow(img); nc <- ncol(img)
cy <- (nr + 1) / 2; cx <- (nc + 1) / 2; th <- angle_deg * pi / 180
xs_ <- .col_idx(nr, nc) - cx; ys_ <- cy - .row_idx(nr, nc)
xr <- cos(th) * xs_ + sin(th) * ys_ # inverse map
yr <- -sin(th) * xs_ + cos(th) * ys_
f <- if (interp == "nearest") warp_nearest else warp_bilinear
f(img, cy - yr, xr + cx, bg)
}par(mfrow = c(1, 4), mar = c(1, 1, 2.5, 1))
show_img(img01, "Original", zlim = c(0, 1))
show_img(translate_image(img01, 20, -12), "Translated (+20, -12)", zlim = c(0, 1))
show_img(rotate_image(img01, 25), "Rotated 25 deg", zlim = c(0, 1))
show_img(rotate_image(img01, 25, interp = "nearest"),
"Rotated 25 deg (NN)", zlim = c(0, 1))
Interpolation is a filter: repeated resampling destroys resolution
Each bilinear resampling convolves the image with a triangular kernel. Applying \(k\) rotations that sum to a full turn should return the original image exactly; in practice the error grows and then saturates as the image loses high-frequency content.
z <- img01; err <- numeric(9)
for (k in 1:9) { z <- rotate_image(z, 40); err[k] <- sqrt(mean((z - img01)^2)) }
## after 9 x 40 = 360 degrees the geometry is back to the start
par(mfrow = c(1, 2))
plot(1:9, err, type = "b", pch = 19, col = "firebrick",
xlab = "number of 40-degree rotations", ylab = "RMS difference from original",
main = "Resampling error accumulates")
show_img(z, "After 9 rotations (net 360 deg)", zlim = c(0, 1))
par(mfrow = c(1, 1))
## how much extra Gaussian blur is that, expressed in pixels?
sig_grid <- seq(0, 2, by = 0.05)
fit_sig <- sig_grid[which.min(sapply(sig_grid,
function(s) mean((gauss_blur(img01, s) - z)^2)))]
c(rms_after_full_turn = round(err[9], 4),
equivalent_Gaussian_sigma_px = fit_sig,
equivalent_FWHM_mm = round(fwhm_of_sigma(fit_sig) * dx, 2))## rms_after_full_turn equivalent_Gaussian_sigma_px
## 0.0713 1.3000
## equivalent_FWHM_mm
## 3.0600The final image is geometrically identical to the original but visibly softer: nine interpolations have cost roughly as much resolution as an extra smoothing step. In a longitudinal study this is the difference between measuring atrophy and measuring your own interpolation.
Intensity normalization and histogram methods
Images from different times, scanners, or settings may not share a common intensity scale. Some modalities are calibrated (CT in HU, quantitative MRI maps, PET SUV); many produce relative values (conventional MRI, most microscopy).
| Method | Definition | Robustness | When appropriate |
|---|---|---|---|
| min–max | \((I-\min)/(\max-\min)\) | poor, one hot voxel rescales everything | clean, artifact-free data |
| percentile (\(p_1\)–\(p_{99}\)) | clip then rescale | good | general display and ML input |
| z-score in a mask | \((I-\mu_M)/\sigma_M\) | good | conventional MRI across subjects |
| reference-tissue scaling | divide by \(\mu\) in a stable tissue (white matter, muscle, liver, cerebellum) | good | SUV ratios, MRI “white-stripe” |
| histogram matching | quantile map onto a template | strong | multi-site harmonization |
| physical calibration | HU, SUV, \(T_1\), ADC | best | when the physics permits it |
minmax01 <- function(x) { r <- range(x, na.rm = TRUE); (x - r[1]) / diff(r) }
pct01 <- function(x, p = c(0.01, 0.99)) {
q <- quantile(x, p, na.rm = TRUE); pmin(pmax((x - q[1]) / diff(q), 0), 1)
}
zscore_in <- function(x, mask) (x - mean(x[mask])) / sd(x[mask])
## inject a single "hot" voxel, as a dead detector element or spike would
img_spike <- img_noisy; img_spike[10, 10] <- 8
par(mfrow = c(1, 2))
hist(minmax01(img_spike), breaks = 60, col = "gray80",
main = "min-max after one spike: data crushed", xlab = "normalized intensity")
hist(pct01(img_spike), breaks = 60, col = "steelblue",
main = "1st-99th percentile normalization", xlab = "normalized intensity")
par(mfrow = c(1, 1))
data.frame(
method = c("min-max", "percentile 1-99", "z-score in head mask"),
parenchyma_mean = round(c(mean(minmax01(img_spike)[par_mask]),
mean(pct01(img_spike)[par_mask]),
mean(zscore_in(img_spike, brain_mask)[par_mask])), 3),
tumour_mean = round(c(mean(minmax01(img_spike)[tumor_mask]),
mean(pct01(img_spike)[tumor_mask]),
mean(zscore_in(img_spike, brain_mask)[tumor_mask])), 3)
)## method parenchyma_mean tumour_mean
## 1 min-max 0.060 0.109
## 2 percentile 1-99 0.509 0.959
## 3 z-score in head mask -0.184 3.875A single outlier voxel compresses the entire min–max normalized image into a narrow band near zero. Percentile normalization is unaffected. This is not a contrived example: hot pixels, metal, and truncation spikes are routine.
Histogram matching
Histogram matching (histogram equalization) maps the quantiles of one image onto those of a reference, forcing the two to share an intensity distribution. It is a strong harmonization tool and a correspondingly strong assumption: it presumes the two images contain the same tissue proportions, so it can erase real biological differences (a large tumour changes the histogram legitimately).
hist_match <- function(src, ref, mask = NULL, nq = 512) {
s <- if (is.null(mask)) as.vector(src) else src[mask]
r <- if (is.null(mask)) as.vector(ref) else ref[mask]
p <- seq(0, 1, length.out = nq)
qs <- quantile(s, p, names = FALSE); qr <- quantile(r, p, names = FALSE)
out <- src
out[] <- approx(qs, qr, xout = as.vector(src), rule = 2)$y
out
}
## a "different scanner": 1.4x gain, +0.05 offset, slightly more noise
img_siteB <- 1.4 * img_noisy + 0.05 + matrix(rnorm(n * n, 0, 0.01), n, n)
matched <- hist_match(img_siteB, img_noisy, brain_mask)
plot(density(img_noisy[brain_mask]), lwd = 2, col = "black", ylim = c(0, 9),
main = "Histogram matching harmonizes two 'sites'", xlab = "intensity")
lines(density(img_siteB[brain_mask]), lwd = 2, col = "firebrick")
lines(density(matched[brain_mask]), lwd = 2, col = "steelblue", lty = 2)
legend("topright", c("site A (reference)", "site B (raw)", "site B (matched)"),
col = c("black", "firebrick", "steelblue"), lty = c(1, 1, 2), lwd = 2, bty = "n")
Bias, background, and illumination correction
A slowly varying multiplicative field, the MRI receive-coil bias field, uneven microscope illumination, residual CT scatter, corrupts every intensity-based analysis. Our phantom contains a known multiplicative field, so we can correct it and measure exactly how well we did.
Two strategies are common, and the difference between them is instructive.
1. Low-pass estimation. Estimate \(b\) by heavily smoothing the image, then divide. This is fast and assumption-light, but the smoothing window must be much larger than any real structure, otherwise the estimate absorbs the anatomy.
2. Parametric fitting in the log domain. Since \(\log(b\cdot u)=\log b+\log u\), a multiplicative field becomes additive under a logarithm. Fitting a smooth low-order polynomial (or B-spline) surface to \(\log I\) inside a tissue mask, with the brightest and darkest percentiles trimmed, gives a well-conditioned estimate. This is the conceptual core of the standard N3/N4 algorithms.
## --- Otsu's threshold (derived in the Segmentation section; used here) -----
otsu_threshold <- function(x, nbins = 256) {
v <- as.vector(x); v <- v[is.finite(v)]
br <- seq(min(v), max(v), length.out = nbins + 1)
h <- hist(v, breaks = br, plot = FALSE)
p <- h$counts / sum(h$counts); mids <- h$mids
w0 <- cumsum(p); mu <- cumsum(p * mids); muT <- mu[length(mu)]
sb <- (muT * w0 - mu)^2 / (w0 * (1 - w0))
sb[!is.finite(sb)] <- -Inf
mids[which.max(sb)]
}
## --- binary morphology with an explicit structuring element ---------------
struct_elem <- function(r, shape = c("disc", "square")) {
shape <- match.arg(shape); ax <- -r:r
if (shape == "square") matrix(TRUE, 2*r+1, 2*r+1)
else outer(ax, ax, function(a, b) a^2 + b^2) <= r^2 + 1e-9
}
bin_dilate <- function(m, se = struct_elem(1)) {
r <- (nrow(se) - 1) %/% 2; nr <- nrow(m); nc <- ncol(m)
P <- pad_image(m * 1, r, r, "zero"); out <- matrix(0, nr, nc)
for (u in 1:nrow(se)) for (v in 1:ncol(se)) if (se[u, v])
out <- pmax(out, P[u:(u + nr - 1), v:(v + nc - 1)])
out > 0
}
bin_erode <- function(m, se = struct_elem(1)) {
r <- (nrow(se) - 1) %/% 2; nr <- nrow(m); nc <- ncol(m)
P <- pad_image(m * 1, r, r, "replicate"); out <- matrix(1, nr, nc)
for (u in 1:nrow(se)) for (v in 1:ncol(se)) if (se[u, v])
out <- pmin(out, P[u:(u + nr - 1), v:(v + nc - 1)])
out > 0
}
bin_open <- function(m, se = struct_elem(1)) bin_dilate(bin_erode(m, se), se)
bin_close <- function(m, se = struct_elem(1)) bin_erode(bin_dilate(m, se), se)
## --- integral-image box mean: O(1) per pixel regardless of window size -----
box_blur <- function(img, k, mode = "replicate") {
p <- (k - 1) %/% 2
P <- pad_image(img, p, p, mode)
S <- t(apply(apply(P, 2, cumsum), 1, cumsum)) # 2D cumulative sum
S <- rbind(0, cbind(0, S))
i1 <- 1:nrow(img); j1 <- 1:ncol(img)
(S[i1 + k, j1 + k] - S[i1, j1 + k] - S[i1 + k, j1] + S[i1, j1]) / (k * k)
}## tissue mask: Otsu, then morphological opening to drop isolated noise specks
head_mask <- bin_open(img_noisy > otsu_threshold(img_noisy), struct_elem(2))
data.frame(threshold = round(otsu_threshold(img_noisy), 3),
head_pixels = sum(head_mask),
truth_pixels = sum(brain_mask),
Dice_vs_truth = round(2 * sum(head_mask & brain_mask) /
(sum(head_mask) + sum(brain_mask)), 4))## threshold head_pixels truth_pixels Dice_vs_truth
## 1 0.235 15347 15456 0.9951## ---- strategy 1: low-pass estimate ---------------------------------------
b_lowpass <- box_blur(img_noisy, 41)
b_lowpass <- b_lowpass / mean(b_lowpass[head_mask])
img_bc_lp <- img_noisy / pmax(b_lowpass, 0.2)
## ---- strategy 2: log-domain polynomial fit inside the tissue mask --------
estimate_bias <- function(img, mask, degree = 2, trim = 0.02) {
## normalised, mask-centerd coordinates for every pixel
u_all <- (as.vector(col(img)) - mean(col(img)[mask])) / nrow(img)
v_all <- (as.vector(row(img)) - mean(row(img)[mask])) / nrow(img)
## fit on interior tissue only: positive, inside the mask, extremes trimmed
idx <- which(as.vector(mask) & as.vector(img) > 0)
y <- log(as.vector(img)[idx])
q <- quantile(y, c(trim, 1 - trim))
keep<- idx[y > q[1] & y < q[2]]
Xfit <- poly(u_all[keep], v_all[keep], degree = degree, raw = TRUE)
beta <- qr.solve(cbind(1, Xfit), log(as.vector(img)[keep]))
Xall <- poly(u_all, v_all, degree = degree, raw = TRUE)
bhat <- matrix(exp(cbind(1, Xall) %*% beta), nrow(img), ncol(img))
bhat / exp(mean(log(bhat[mask]))) # unit geometric mean in mask
}
b_poly <- estimate_bias(img_noisy, head_mask, degree = 2)
img_bc <- img_noisy / b_poly
b_true <- bias_field / exp(mean(log(bias_field[head_mask])))
par(mfrow = c(2, 3), mar = c(1, 1, 2.5, 1))
show_img(img_noisy, "Observed (biased)", zlim = c(0, 1))
show_img(b_true, "TRUE bias field", zlim = c(0.7, 1.3))
show_img(b_poly, "Estimated (log-polynomial)", zlim = c(0.7, 1.3))
show_img(img_bc, "Corrected (polynomial)", zlim = c(0, 1))
show_img(b_lowpass, "Estimated (41-px low-pass)", zlim = c(0.7, 1.3))
show_img(img_bc_lp, "Corrected (low-pass)", zlim = c(0, 1))
left <- col(img_true) <= n / 2; right <- !left
## interior parenchyma only: erode away voxels contaminated by PSF spill-in
par_core <- bin_erode(par_mask, struct_elem(3))
row_of <- function(x) c(
parenchyma_SD = sd(x[par_core]),
left_mean = mean(x[par_core & left]),
right_mean = mean(x[par_core & right]),
LR_difference_pct = 100 * (mean(x[par_core & right]) / mean(x[par_core & left]) - 1))
res <- round(rbind(
observed = row_of(img_noisy),
lowpass_corrected = row_of(img_bc_lp),
polynomial_corrected= row_of(img_bc)), 4)
res## parenchyma_SD left_mean right_mean LR_difference_pct
## observed 0.0531 0.3773 0.4243 12.4689
## lowpass_corrected 0.0892 0.4176 0.4090 -2.0673
## polynomial_corrected 0.0421 0.3996 0.3987 -0.2256c(noise_floor_SD = sigma_noise,
bias_field_correlation_lowpass = round(cor(b_lowpass[head_mask], b_true[head_mask]), 3),
bias_field_correlation_poly = round(cor(b_poly[head_mask], b_true[head_mask]), 3))## noise_floor_SD bias_field_correlation_lowpass
## 0.040 0.355
## bias_field_correlation_poly
## 0.938The numbers tell the story. In the observed image the parenchyma, which is uniform by construction, shows a spurious left–right intensity difference of over 10 % and a standard deviation well above the noise floor. The polynomial correction removes essentially all of it, driving the parenchyma SD down to the true noise level \(\sigma\). The low-pass estimate is contaminated by the anatomy it was meant to ignore.
Limitations of low-pass bias estimation. A background estimated by smoothing also absorbs large, smooth real structures. The smoothing window must be larger than the objects of interest but smaller than the scale of the artifact, a judgement call that must be checked against the original image. Parametric (polynomial or B-spline) fitting inside a tissue mask, with extremes trimmed, is both more stable and easier to report.
Bias correction is not free. Dividing by an estimated field also divides the noise, making it spatially non-stationary: after correction, \(\sigma\) varies across the image as \(\sigma/\hat b(x)\). Any later statistic that assumes homoscedastic noise (a global threshold, a \(t\)-map) must account for this.
Know your data: modality-specific caution
The same operation means different things in different modalities. Thresholding is a single mathematical operation, but the thresholded quantity is Hounsfield units in CT, sequence-dependent intensity or quantitative maps in MRI, standardized uptake values in PET, echo intensity in ultrasound, and fluorescence intensity in microscopy.
| Modality | Physical signal / contrast | Example quantitative unit | Example use |
|---|---|---|---|
| X-ray radiography | X-ray attenuation (projection) | detector counts, gray levels | fracture visibility, lung opacity |
| CT | X-ray attenuation in 3D | Hounsfield units (HU) | tissue density, lesion volume, calcium score |
| MRI | proton signal, relaxation, diffusion, flow | \(T_1\), \(T_2\), \(T_2^*\), ADC, PD | brain volume, edema, tumour contrast |
| Ultrasound | acoustic reflection / scattering | relative intensity, velocity (m/s) | organ motion, vessel flow, elastography |
| PET / SPECT | radiotracer emission | SUV, binding potential, \(K_i\) | metabolism, perfusion, receptor binding |
| Optical microscopy | absorption, fluorescence, scattering | intensity, photon counts | cell count, morphology, localization |
Key concept. The mathematical operations are shared across modalities, but the biological interpretation of a processed image depends on the physics of the signal, and so does the legitimacy of a given preprocessing step. Intensity normalization is mandatory for conventional MRI and forbidden for CT, where it would destroy the very calibration that makes HU meaningful.
Noise models and the denoising–enhancement trade-off
Denoising reduces random variation; enhancement improves the visibility of structures. The central challenge is to reduce noise without erasing small structures, and the correct denoiser depends on the noise model established earlier: Gaussian noise calls for linear or bilateral smoothing, Poisson noise for a variance-stabilizing transform first, impulse noise for a rank filter, and Rician noise for a bias-aware estimator.
A useful rule. Denoising should reduce noise more than it reduces biologically meaningful structure. Choose and report the method, and always check its effect on the downstream measurement that actually matters, not on how the image looks.
Checkpoint 4. You must resample a \(T_1\) MRI and its manually drawn tumour mask from \(0.9\times0.9\times3\) mm to \(1\times1\times1\) mm isotropic, then rotate into atlas space. What interpolation do you use for each, and how many times do you resample? Answer: linear (or B-spline) for the image, nearest-neighbor for the mask; compose the resample and the rotation into a single affine and apply it once to each.
Image filtering
Medical images contain noise, artifacts, low contrast, and blurred boundaries. Filtering prepares images for visualization, segmentation, registration, and measurement by reducing noise, suppressing outliers, enhancing edges, emphasizing structures of a given size, or improving downstream segmentation.
Important caution. Filtering can improve analysis but it also changes the image. Choose a filter based on the measurement goal, quantify what it does to both noise and resolution, and apply the same pipeline consistently across every image that will be compared.
Why convolution matters
Medical images contain structure at many spatial scales, broad tissue regions, organ boundaries, small vessels, and noise in CT; cells, nuclei, membranes, and subcellular puncta in microscopy. Convolution is the core operation behind most linear filters: by changing the convolution kernel we can smooth noise, sharpen edges, detect boundaries, or emphasize structures of a certain size or direction. And image blur is itself a convolution (with the PSF), so understanding convolution explains both how images degrade and how we repair them.
Neighborhood operations and weighted averaging
The intuitive starting point is local averaging: replace each pixel by the average of a window around it. For a 2D image \(f(x,y)\) and a \((2a+1)\times(2a+1)\) window,
\[f'(x,y) = \frac{1}{N}\sum_{s=-a}^{a}\sum_{t=-a}^{a} f(x+s,\,y+t), \qquad N=(2a+1)^2 .\]
This reduces pixel-to-pixel noise but blurs edges and small structures. The general version uses a matrix of weights \(w(s,t)\), the kernel, giving the convolution
\[(w * f)(x,y) = \sum_{s}\sum_{t} w(s,t)\, f(x-s,\,y-t).\]
One- and two-dimensional convolution
For a 1D signal \(f[n]\) and kernel \(h[n]\), the discrete convolution is
\[(f * h)[n] = \sum_{m} f[m]\,h[n-m].\]
For a 2D image \(I[i,j]\) and kernel \(K[u,v]\),
\[(I * K)[i,j] = \sum_{u}\sum_{v} I[i-u,\,j-v]\,K[u,v].\]
The kernel slides across the image; at each location the overlapping intensities are multiplied by kernel weights and summed.
Two structural properties.
Linearity: \(w * (c_1 f_1 + c_2 f_2) = c_1 (w*f_1) + c_2 (w*f_2)\).
Shift invariance: if \(f'(x)=f(x-x_0)\) then \((w*f')(x)=(w*f)(x-x_0)\).
Together these define a linear shift-invariant (LSI) system. Every LSI system is a convolution, and, as the convolution theorem will show, every convolution is a multiplication in the frequency domain. This is why the PSF/MTF pair completely characterizes an LSI imaging system.
Correlation vs. convolution. Strict convolution flips the kernel before applying it; correlation does not. For symmetric kernels (mean, Gaussian, Laplacian) the two coincide. For asymmetric directional kernels (Sobel, Prewitt) the convention flips the sign of the response, so check which one your software implements. Deep-learning “convolution” layers are in fact correlations.
Our conv2() (defined with the phantom) flips the kernel and therefore performs
a true convolution. It is vectorized over kernel entries rather than pixels,
which makes it exact and fast enough for interactive use.
## sanity check against an explicit, obviously-correct double loop
naive_conv2 <- function(img, kernel) {
nr <- nrow(img); nc <- ncol(img); kr <- nrow(kernel); kc <- ncol(kernel)
pr <- (kr - 1) %/% 2; pc <- (kc - 1) %/% 2
P <- pad_image(img, pr, pc, "replicate"); kf <- kernel[kr:1, kc:1]
out <- matrix(0, nr, nc)
for (i in 1:nr) for (j in 1:nc)
out[i, j] <- sum(P[i:(i + kr - 1), j:(j + kc - 1)] * kf)
out
}
Ktest <- matrix(c(-1, 0, 1, -2, 0, 2, -1, 0, 1), 3, 3, byrow = TRUE)
c(max_abs_difference = max(abs(conv2(img01, Ktest) - naive_conv2(img01, Ktest))))## max_abs_difference
## 2.22e-16Mean filtering and the noise-propagation law
A mean filter uses equal weights; the \(3\times 3\) version is
\[K = \frac{1}{9}\begin{bmatrix}1&1&1\\1&1&1\\1&1&1\end{bmatrix}.\]
Its entries sum to 1 so the filter preserves the overall intensity level (a flat region maps to itself, i.e. the filter has unit DC gain).
How much noise does a linear filter remove? If the input noise is white with variance \(\sigma^2\), then the output of convolution with kernel \(w\) has variance
\[\sigma_{\text{out}}^2 = \sigma^2\sum_{u,v} w(u,v)^2 = \sigma^2\,\lVert w\rVert_2^2 .\]
For a \(k\times k\) mean filter \(\lVert w\rVert_2^2 = 1/k^2\), so \(\sigma_{\text{out}} = \sigma/k\), noise falls with the linear size of the window, not its area. For a Gaussian of width \(\sigma_g\) pixels, \(\lVert w\rVert_2^2 \to 1/(4\pi\sigma_g^2)\), so \(\sigma_{\text{out}} \approx \sigma/(2\sigma_g\sqrt{\pi})\).
The catch: the output noise is no longer white. It is correlated over the kernel’s support, which is exactly why naive standard errors computed on smoothed images are badly wrong (see the uncertainty section).
set.seed(202)
pure <- matrix(rnorm(300 * 300, 0, sigma_noise), 300, 300) # noise only
core <- 40:260 # avoid edge effects
tab <- rbind(
t(sapply(c(3, 5, 7, 9), function(k) {
Kk <- matrix(1, k, k) / k^2
c(filter = k, empirical_SD = sd(conv2(pure, Kk)[core, core]),
predicted = sigma_noise * sqrt(sum(Kk^2)),
simple_rule = sigma_noise / k)
})))
round(as.data.frame(tab), 5)## filter empirical_SD predicted simple_rule
## 1 3 0.01325 0.01333 0.01333
## 2 5 0.00788 0.00800 0.00800
## 3 7 0.00560 0.00571 0.00571
## 4 9 0.00436 0.00444 0.00444## Gaussian kernels
gk <- function(s) { h <- ceiling(3*s); g <- exp(-(-h:h)^2/(2*s^2)); g <- g/sum(g); outer(g, g) }
round(as.data.frame(t(sapply(c(0.8, 1.2, 2.0, 3.0), function(s) {
K <- gk(s)
c(sigma_g = s, empirical_SD = sd(gauss_blur(pure, s)[core, core]),
predicted = sigma_noise * sqrt(sum(K^2)),
asymptotic = sigma_noise / (2 * s * sqrt(pi)))
}))), 5)## sigma_g empirical_SD predicted asymptotic
## 1 0.8 0.01407 0.01416 0.01410
## 2 1.2 0.00931 0.00941 0.00940
## 3 2.0 0.00554 0.00565 0.00564
## 4 3.0 0.00368 0.00377 0.00376Prediction and experiment agree to within a percent, which means the noise consequence of any linear filter can be computed before running it.
img_mean5 <- conv2(img_bc, matrix(1, 5, 5) / 25)
par(mfrow = c(1, 3), mar = c(1, 1, 2.5, 1))
show_img(img_bc, "Bias-corrected input", zlim = c(0, 1))
show_img(img_mean5,"Mean filtered (5x5)", zlim = c(0, 1))
show_img(img_bc - img_mean5, "Removed component", zlim = c(-0.15, 0.15))
The “removed” panel is the diagnostic that matters: it should look like noise. If it shows anatomy, edges, the vessel, lesion rims, the filter is deleting signal, not noise.
Common pitfall. A large mean filter makes an image look smoother but can erase small lesions, thin vessels, cell boundaries, or fine anatomy. Always inspect the difference image.
Gaussian smoothing
A Gaussian kernel weights nearby pixels more than distant ones:
\[G(x,y) = \frac{1}{2\pi\sigma_g^2} \exp\!\left(-\frac{x^2+y^2}{2\sigma_g^2}\right),\]
with \(\sigma_g\) controlling the width. It is the standard pre-segmentation denoiser for approximately Gaussian noise, and it has three properties no other kernel combines: it is separable (a 2D Gaussian is the product of two 1D Gaussians, reducing cost from \(O(k^2)\) to \(O(k)\) per pixel), it is the unique kernel that creates no new extrema as \(\sigma_g\) grows (the basis of scale-space theory), and its Fourier transform is again a Gaussian, so it introduces no ringing.
\[\text{FWHM}_{\text{after}}^2 = \text{FWHM}_{\text{before}}^2 + \text{FWHM}_{\text{kernel}}^2 .\]
The matched-filter theorem, and how much to smooth. To detect a known signal \(s\) in additive white noise, the SNR-optimal linear filter is a copy of the signal itself. For a roughly Gaussian lesion of width \(\sigma_{\text{les}}\), this means smoothing with \(\sigma_g \approx \sigma_{\text{les}}\), smoothing at the scale of the object you are looking for. Smoothing much less leaves noise; smoothing much more destroys the object. This is the quantitative justification for the 6–8 mm smoothing kernels conventional in fMRI and PET, and for not smoothing when the goal is boundary delineation.
par(mfrow = c(1, 4), mar = c(1, 1, 2.5, 1))
show_img(img_bc, "Input", zlim = c(0, 1))
show_img(gauss_blur(img_bc, 0.8), expression(sigma[g] == 0.8), zlim = c(0, 1))
show_img(gauss_blur(img_bc, 2.0), expression(sigma[g] == 2.0), zlim = c(0, 1))
show_img(gauss_blur(img_bc, 4.0), expression(sigma[g] == 4.0), zlim = c(0, 1))
## Detectability of the SMALL lesion versus smoothing width.
## Signal = peak value at the lesion center minus the parenchymal level.
## Noise = SD of the smoothed noise field, predicted exactly as sigma*||w||_2.
gk_2d <- function(s) { h <- ceiling(3*s); g <- exp(-(-h:h)^2/(2*s^2))
g <- g/sum(g); outer(g, g) }
par_core <- bin_erode(par_mask, struct_elem(3))
cen <- round(colMeans(which(lesion_mask, arr.ind = TRUE)))
sig_grid <- seq(0, 4.5, by = 0.25)
det_snr <- sapply(sig_grid, function(s) {
z <- if (s == 0) img_bc else gauss_blur(img_bc, s)
nz <- if (s == 0) sigma_noise else sigma_noise * sqrt(sum(gk_2d(s)^2))
(z[cen[1], cen[2]] - mean(z[par_core])) / nz
})
R_les <- sqrt(sum(lesion_mask) / pi) # equivalent radius, px
sigma_pred <- sqrt((R_les / 2)^2 + psf_sigma^2) # object scale after the PSF
plot(sig_grid, det_snr, type = "b", pch = 19, col = "darkgreen",
xlab = expression(smoothing~sigma[g]~(pixels)),
ylab = "detection SNR at the lesion center",
main = "Matched filtering: smooth at the object's own scale")
abline(v = sigma_pred, col = "red", lty = 2)
abline(v = sig_grid[which.max(det_snr)], col = "steelblue", lty = 3)
legend("bottomright", bty = "n",
legend = c(sprintf("predicted object scale = %.1f px", sigma_pred),
sprintf("empirical optimum = %.2f px", sig_grid[which.max(det_snr)])),
col = c("red", "steelblue"), lty = c(2, 3))
Detectability rises steeply, peaks in a broad plateau near the object’s own scale, and then decays as the lesion is smoothed away. The optimum is flat — which is fortunate, because it means the exact kernel width matters much less than being in the right order of magnitude. Smoothing at 0.25 px or at 10 px would both be badly wrong.
Nonlinear filters
Median filtering
A median filter replaces each pixel with the median of its neighborhood. Because the median is a rank statistic, not a weighted sum, median filtering is nonlinear and is not a convolution:
\[I_{\text{filtered}}[i,j] = \operatorname{median}\{\,I[u,v] : (u,v) \in \mathcal{N}_{i,j}\,\}.\]
For a neighborhood \(\{2,3,3,4,5,6,7,100,105\}\) the median is 5; the outliers 100 and 105 have almost no effect. This makes median filtering ideal for impulse/salt-and-pepper noise while preserving step edges, a step edge is exactly preserved by a median filter of any size, whereas any linear smoother blurs it.
## For the 3x3 case the median of 9 values is obtained by a sorting network:
## sort each of three triples, then take the median of (max of minima,
## median of medians, min of maxima). This is exact and fully vectorised.
.sort3 <- function(a, b, c) {
t <- pmin(a, b); b <- pmax(a, b); a <- t
t <- pmin(b, c); c <- pmax(b, c); b <- t
t <- pmin(a, b); b <- pmax(a, b); a <- t
list(a, b, c)
}
median_filter <- function(img, size = 3, mode = "replicate") {
p <- (size - 1) %/% 2
P <- pad_image(img, p, p, mode)
nr <- nrow(img); nc <- ncol(img)
sh <- vector("list", size * size); m <- 0L
for (u in 1:size) for (v in 1:size) {
m <- m + 1L; sh[[m]] <- P[u:(u + nr - 1), v:(v + nc - 1)]
}
if (size == 3) {
A <- .sort3(sh[[1]], sh[[2]], sh[[3]])
B <- .sort3(sh[[4]], sh[[5]], sh[[6]])
C <- .sort3(sh[[7]], sh[[8]], sh[[9]])
return(.sort3(pmax(A[[1]], B[[1]], C[[1]]),
.sort3(A[[2]], B[[2]], C[[2]])[[2]],
pmin(A[[3]], B[[3]], C[[3]]))[[2]])
}
S <- matrix(unlist(lapply(sh, as.vector)), ncol = size * size)
matrix(apply(S, 1, median), nr, nc)
}
## salt-and-pepper corruption, then mean vs median
set.seed(5)
sp <- img_bc
hit <- sample(length(sp), round(0.03 * length(sp)))
sp[hit] <- ifelse(runif(length(hit)) > 0.5, 1, 0)
c(mean_filter_RMSE = round(sqrt(mean((conv2(sp, matrix(1,3,3)/9) - img_bc)^2)), 4),
median_filter_RMSE = round(sqrt(mean((median_filter(sp, 3) - img_bc)^2)), 4))## mean_filter_RMSE median_filter_RMSE
## 0.0528 0.0414par(mfrow = c(1, 4), mar = c(1, 1, 2.5, 1))
show_img(sp, "3% salt & pepper", zlim = c(0, 1))
show_img(conv2(sp, matrix(1, 3, 3) / 9), "Mean 3x3 (smears)", zlim = c(0, 1))
show_img(median_filter(sp, 3), "Median 3x3", zlim = c(0, 1))
show_img(median_filter(img_bc, 3), "Median on clean input", zlim = c(0, 1))
Bilateral filtering: edge-preserving smoothing
A bilateral filter weights neighbors by both spatial distance and intensity similarity,
\[I'(p) = \frac{1}{Z_p}\sum_{q\in\mathcal N_p} \underbrace{e^{-\lVert p-q\rVert^2/2\sigma_s^2}}_{\text{spatial}}\; \underbrace{e^{-(I(p)-I(q))^2/2\sigma_r^2}}_{\text{range}}\;I(q), \qquad Z_p=\textstyle\sum_q(\cdots),\]
so voxels across a strong boundary contribute almost nothing. It smooths within tissue while leaving boundaries essentially untouched, the same idea underlies non-local means and, in a different guise, anisotropic diffusion and total-variation denoising.
Note that in the bilateral filter, the normalization factor, i.e., partition function, is \[{Z_p}=\sum_{q\in\mathcal N_p} \underbrace{e^{-\lVert p-q\rVert^2/2\sigma_s^2}}_{\text{spatial}}\; \underbrace{e^{-(I(p)-I(q))^2/2\sigma_r^2}}_{\text{range}}.\]
bilateral_filter <- function(img, sigma_s = 2, sigma_r = 0.05,
size = 2 * ceiling(2 * sigma_s) + 1, mode = "replicate") {
p <- (size - 1) %/% 2
P <- pad_image(img, p, p, mode)
nr <- nrow(img); nc <- ncol(img)
num <- matrix(0, nr, nc); den <- matrix(0, nr, nc)
for (u in -p:p) for (v in -p:p) {
shifted <- P[(u + p + 1):(u + p + nr), (v + p + 1):(v + p + nc)]
w <- exp(-(u^2 + v^2) / (2 * sigma_s^2)) *
exp(-(shifted - img)^2 / (2 * sigma_r^2))
num <- num + w * shifted; den <- den + w
}
num / den
}
img_bl <- bilateral_filter(img_bc, sigma_s = 2, sigma_r = 0.06)img_gs <- gauss_blur(img_bc, 1.5)
img_md <- median_filter(img_bc, 3)
par(mfrow = c(2, 2), mar = c(1, 1, 2.5, 1))
show_img(img_bc, "Input (bias-corrected)", zlim = c(0, 1))
show_img(img_gs, "Gaussian (sigma = 1.5)", zlim = c(0, 1))
show_img(img_md, "Median 3x3", zlim = c(0, 1))
show_img(img_bl, "Bilateral", zlim = c(0, 1))
par(mfrow = c(1, 1))
## quantitative comparison: noise suppression vs edge preservation
edge_band <- bin_dilate(tumor_mask, struct_elem(2)) & !bin_erode(tumor_mask, struct_elem(2))
score <- function(z, name) data.frame(
filter = name,
noise_SD_in_parenchyma = round(sd(z[bin_erode(par_mask, struct_elem(3))]), 4),
edge_gradient = round(mean(abs(conv2(z, Ktest))[edge_band]), 4),
lesion_contrast = round(mean(z[bin_erode(lesion_mask, struct_elem(1))]) -
mean(z[bin_erode(par_mask, struct_elem(3))]), 4))
do.call(rbind, list(score(img_bc, "none"), score(img_gs, "Gaussian 1.5"),
score(img_md, "median 3x3"), score(img_bl, "bilateral")))## filter noise_SD_in_parenchyma edge_gradient lesion_contrast
## 1 none 0.0421 0.4096 0.1409
## 2 Gaussian 1.5 0.0132 0.3005 0.1220
## 3 median 3x3 0.0199 0.3752 0.1370
## 4 bilateral 0.0194 0.4037 0.1299Reading the comparison table. The ideal denoiser has a low noise SD, a high edge gradient, and an unchanged lesion contrast. Gaussian smoothing wins on noise and loses on edges; the median filter is a compromise; the bilateral filter suppresses noise while retaining the boundary, at the cost of two extra parameters (\(\sigma_s,\sigma_r\)) that must be reported and held fixed across a study, because they silently change every downstream measurement.
| Feature | Gaussian smoothing | Median filtering | Bilateral filtering |
|---|---|---|---|
| Filter type | Linear | Nonlinear (rank) | Nonlinear (data-adaptive) |
| Convolution? | Yes | No | No |
| Best for | Gaussian noise | Impulse / outlier noise | Noise with edges to keep |
| Edge behavior | Blurs | Preserves step edges | Preserves strong edges |
| Noise SD prediction | Exact (\(\sigma\lVert w\rVert_2\)) | Approximate | Not available in closed form |
| Main risk | Erases fine detail | Deletes small isolated objects | Cartoon/staircase artifacts |
We carry img_dn <- median_filter(img_bc, 3) forward as the denoised image
for the segmentation and measurement sections: it removes pixel-scale noise while
keeping object boundaries sharp, and, unlike the bilateral filter, it has no
tuning parameters beyond the window size.
Edge enhancement and sharpening
Edges are rapid changes in intensity. The Sobel kernels approximate the horizontal and vertical intensity gradients,
\[K_x = \begin{bmatrix}-1&0&1\\-2&0&2\\-1&0&1\end{bmatrix}, \qquad K_y = \begin{bmatrix}-1&-2&-1\\0&0&0\\1&2&1\end{bmatrix},\]
with gradient magnitude \(|\nabla I| = \sqrt{(I*K_x)^2 + (I*K_y)^2}\) and orientation \(\theta = \operatorname{atan2}(I*K_y,\, I*K_x)\).
Calibrating a gradient. The Sobel kernel is a scaled derivative: it equals \(8\Delta\) times the true directional derivative for a unit-spacing grid (2 from the central difference, 4 from the smoothing column). To report a gradient in physical units (intensity per millimeter), divide by \(8\Delta x\). Reporting an uncalibrated “edge strength” is one of the commonest ways image-derived features become non-comparable across studies.
The Laplacian \(\nabla^2 I\) responds to second derivatives and highlights boundaries in all directions. Note the sign convention: the discrete Laplacian is
\[\nabla^2 \approx \begin{bmatrix}0&1&0\\1&-4&1\\0&1&0\end{bmatrix},\]
so the classical unsharp mask subtracts it, \(I_{\text{sharp}} = I - \alpha\,\nabla^2 I\), equivalently \(I_{\text{sharp}} = I + \alpha\,(I - G_\sigma * I)\) for a Gaussian \(G_\sigma\).
Kx <- matrix(c(-1, 0, 1, -2, 0, 2, -1, 0, 1), 3, 3, byrow = TRUE)
Ky <- matrix(c(-1, -2, -1, 0, 0, 0, 1, 2, 1), 3, 3, byrow = TRUE)
gx <- conv2(img_dn, Kx) / (8 * dx) # calibrated: intensity per mm
gy <- conv2(img_dn, Ky) / (8 * dy)
grad_mag <- sqrt(gx^2 + gy^2)
laplacian <- matrix(c(0, 1, 0, 1, -4, 1, 0, 1, 0), 3, 3, byrow = TRUE)
img_sharp <- img_dn - 0.7 * conv2(img_dn, laplacian)
## Laplacian-of-Gaussian: smooth first, so the second derivative is not noise
LoG <- conv2(gauss_blur(img_dn, 1.5), laplacian)
par(mfrow = c(1, 4), mar = c(1, 1, 2.5, 1))
show_img(img_dn, "Denoised", zlim = c(0, 1))
show_img(normalize01(grad_mag), "Sobel gradient magnitude", zlim = c(0, 1))
show_img(normalize01(img_sharp), "Unsharp mask (I - 0.7*Lap I)",zlim = c(0, 1))
show_img(normalize01(LoG), "Laplacian of Gaussian", zlim = c(0, 1))
par(mfrow = c(1, 1))
c(max_gradient_per_mm = round(max(grad_mag), 3),
noise_SD_before = round(sd(img_dn[bin_erode(par_mask, struct_elem(3))]), 4),
noise_SD_after_sharpening = round(sd(img_sharp[bin_erode(par_mask, struct_elem(3))]), 4))## max_gradient_per_mm noise_SD_before noise_SD_after_sharpening
## 0.1370 0.0199 0.0428Sharpening amplified the parenchymal noise substantially, which is the price of boosting high frequencies. That is acceptable for visual boundary delineation and unacceptable before an intensity measurement.
Clinical use. Edge enhancement helps visualize organ boundaries, vessel walls, lesion margins, or cell boundaries before segmentation. Use it cautiously when quantitative intensity measurements are required, because sharpening also amplifies noise and can create artificial edges (overshoot at every boundary — the spatial-domain face of Gibbs ringing).
Morphological operations
Binary images contain only two classes, foreground (1) and background (0). They support nonlinear morphological operations defined by a small mask, the structuring element \(B\):
- erosion \(A\ominus B\) shrinks foreground (a pixel survives only if the whole structuring element fits inside the object), removes thin protrusions and small specks;
- dilation \(A\oplus B\) grows foreground, fills gaps and connects nearby pieces;
- opening \(A\circ B = (A\ominus B)\oplus B\) removes small objects while approximately preserving the size of large ones;
- closing \(A\bullet B = (A\oplus B)\ominus B\) fills small holes and gaps.
Opening and closing are idempotent (\(A\circ B\circ B = A\circ B\)) and size-selective: opening with a disc of radius \(r\) deletes everything that cannot contain a disc of radius \(r\). That makes them a scale filter, not merely cosmetic cleanup.
noisy_mask <- img_dn > 0.5
par(mfrow = c(1, 4), mar = c(1, 1, 2.5, 1))
show_img(noisy_mask, "Raw mask (I > 0.5)", col = c("black","white"))
show_img(bin_erode(noisy_mask, struct_elem(2)), "Eroded (r = 2)", col = c("black","white"))
show_img(bin_open(noisy_mask, struct_elem(2)), "Opened (r = 2)", col = c("black","white"))
show_img(bin_close(noisy_mask, struct_elem(3)), "Closed (r = 3)", col = c("black","white"))
par(mfrow = c(1, 1))
data.frame(operation = c("raw", "erode r2", "open r2", "close r3"),
pixels = c(sum(noisy_mask), sum(bin_erode(noisy_mask, struct_elem(2))),
sum(bin_open(noisy_mask, struct_elem(2))),
sum(bin_close(noisy_mask, struct_elem(3)))))## operation pixels
## 1 raw 863
## 2 erode r2 458
## 3 open r2 642
## 4 close r3 879Morphology changes measurements. Erosion with a disc of radius \(r\) removes a boundary shell of thickness \(r\): for a convex object of area \(A\) and perimeter \(P\), \(\Delta A \approx -rP\). Opening a 10 mm lesion with a 2 mm element can cost several percent of its area. Cleaning operations belong in the reported pipeline and in the error budget, not in an unlogged “tidy-up” step.
library(EBImage)
shapes <- readImage(system.file("images", "shapes.png", package = "EBImage"))
logo <- shapes[110:512, 1:130]
kern <- makeBrush(5, shape = "diamond")
EBImage::display(combine(erode(logo, kern), dilate(logo, kern)),
method = "raster", all = TRUE)
Boundary effects and padding
When a kernel is centered near an image edge, part of its neighborhood falls outside the image and the software must decide how to fill those missing values. Common strategies are zero padding, mirror/reflection padding, replicate-edge padding, circular/periodic padding, and cropping the output. Boundary handling biases measurements near the field-of-view edge, and, in Fourier-domain filtering, circular padding is implicit, which wraps the opposite edge of the image into the calculation.
K21 <- matrix(1, 21, 21) / 441
par(mfrow = c(1, 3), mar = c(1, 1, 2.5, 1))
show_img(conv2(img_dn, K21, "zero"), "Zero padding (dark rim)", zlim = c(0, 1))
show_img(conv2(img_dn, K21, "replicate"), "Replicate padding", zlim = c(0, 1))
show_img(conv2(img_dn, K21, "reflect"), "Reflect padding", zlim = c(0, 1))
par(mfrow = c(1, 1))
edge_strip <- cbind(1:n, rep(3, n))
data.frame(padding = c("zero", "replicate", "reflect"),
mean_in_left_edge_column = round(c(
mean(conv2(img_dn, K21, "zero")[edge_strip]),
mean(conv2(img_dn, K21, "replicate")[edge_strip]),
mean(conv2(img_dn, K21, "reflect")[edge_strip])), 4),
true_mean = round(mean(img_dn[edge_strip]), 4))## padding mean_in_left_edge_column true_mean
## 1 zero 0.0340 0.0582
## 2 replicate 0.0569 0.0582
## 3 reflect 0.0569 0.0582Deconvolution: partially undoing the PSF
If blur is convolution with a known PSF \(h\), why not divide it out? In the frequency domain the naive inverse filter is \(\hat F = G/H\). It fails catastrophically because \(H(f)\to 0\) at high frequencies, so the ratio amplifies noise without bound. The Wiener filter regularizes it by trading bias against variance:
\[\hat F(f) = \frac{H^*(f)}{|H(f)|^2 + K}\,G(f), \qquad K \approx \frac{S_{\eta}(f)}{S_{f}(f)} \;=\; \text{noise-to-signal power ratio}.\]
For \(K\to0\) this is the inverse filter; for large \(K\) it becomes a smoother.
## build the transfer function H of the known PSF on the image grid
psf_on_grid <- function(sigma, n) {
h <- ceiling(3 * sigma); ax <- -h:h
g <- exp(-ax^2 / (2 * sigma^2)); g <- g / sum(g)
P <- matrix(0, n, n)
ctr<- n %/% 2 + 1
P[(ctr - h):(ctr + h), (ctr - h):(ctr + h)] <- outer(g, g)
P
}
## centring helpers for the DFT (used again throughout the Fourier section).
## fftshift moves the zero frequency to the center; ifftshift is its inverse
## (they differ for odd dimensions, so both are defined explicitly).
fftshift <- function(x) {
nr <- nrow(x); nc <- ncol(x)
x[c((floor(nr/2) + 1):nr, 1:floor(nr/2)),
c((floor(nc/2) + 1):nc, 1:floor(nc/2))]
}
ifftshift <- function(x) {
nr <- nrow(x); nc <- ncol(x)
x[c((ceiling(nr/2) + 1):nr, 1:ceiling(nr/2)),
c((ceiling(nc/2) + 1):nc, 1:ceiling(nc/2))]
}
Hf <- fft(ifftshift(psf_on_grid(psf_sigma, n)))
wiener <- function(g, H, K) Re(fft(fft(g) * Conj(H) / (Mod(H)^2 + K),
inverse = TRUE)) / length(g)
rmse <- function(a, b) sqrt(mean((a - b)^2))
## a single illustration at moderate noise
g_obs <- img_bc
rec_inv <- wiener(g_obs, Hf, 1e-6) # essentially the inverse filter
rec_w <- wiener(g_obs, Hf, 0.03)
par(mfrow = c(1, 4), mar = c(1, 1, 2.5, 1))
show_img(img_true, "Truth f", zlim = c(0, 1))
show_img(g_obs, "Observed (blur + noise)", zlim = c(0, 1))
show_img(rec_inv, "Inverse filter: unusable",zlim = c(0, 1))
show_img(rec_w, "Wiener (K = 0.03)", zlim = c(0, 1))
Whether deconvolution helps depends entirely on the SNR of the input. The table below repeats the experiment at four noise levels, reporting the RMSE against the known truth.
Ks <- c(1e-6, 1e-3, 1e-2, 3e-2, 1e-1, 3e-1)
study <- t(sapply(c(0.005, 0.010, 0.020, 0.040), function(sg) {
set.seed(9)
g <- gauss_blur(img_true, psf_sigma) + matrix(rnorm(n * n, 0, sg), n, n)
c(noise_SD = sg, observed = rmse(g, img_true),
setNames(sapply(Ks, function(k) rmse(wiener(g, Hf, k), img_true)),
paste0("K=", format(Ks, scientific = TRUE))))
}))
round(as.data.frame(study), 4)## noise_SD observed K=1e-06 K=1e-03 K=1e-02 K=3e-02 K=1e-01 K=3e-01
## 1 0.005 0.0310 1.097 0.0414 0.0251 0.0269 0.0412 0.0834
## 2 0.010 0.0322 2.194 0.0764 0.0324 0.0294 0.0417 0.0835
## 3 0.020 0.0365 4.387 0.1495 0.0520 0.0375 0.0435 0.0837
## 4 0.040 0.0503 8.775 0.2974 0.0966 0.0599 0.0501 0.0848cbind(as.data.frame(study)[, 1:2],
best_K = Ks[apply(study[, -(1:2)], 1, which.min)],
best_RMSE = round(apply(study[, -(1:2)], 1, min), 4),
improvement_pct = round(100 * (1 - apply(study[, -(1:2)], 1, min) / study[, 2]), 1))## noise_SD observed best_K best_RMSE improvement_pct
## 1 0.005 0.03097 0.01 0.0251 18.8
## 2 0.010 0.03216 0.03 0.0294 8.6
## 3 0.020 0.03652 0.03 0.0375 -2.8
## 4 0.040 0.05029 0.10 0.0501 0.5At \(\sigma = 0.005\) Wiener deconvolution reduces the error by roughly 20 %; by \(\sigma = 0.04\) there is nothing left to gain, because the noise the filter amplifies exceeds the resolution it restores. Deconvolution buys resolution with SNR, and only works when there is SNR to spend.
Deconvolution is regularized inversion, not magic. The best achievable \(K\) recovers some resolution and cannot recover frequencies where \(H\approx0\), those were destroyed at acquisition. Every deconvolution therefore assumes something (smoothness, positivity, sparsity) to fill the gap, and the assumption, not the data, determines what appears there. Report the PSF model and the regularization strength, and never quantify texture on a deconvolved image.
Choosing a filter
| Goal | Good choice | Reason |
|---|---|---|
| Reduce Gaussian noise | Gaussian smoothing at the object scale | matched filter; predictable \(\sigma_{\text{out}}\) |
| Reduce Poisson noise | Anscombe transform, then Gaussian | stabilizes variance first |
| Remove isolated bright/dark pixels | Median filter | rank statistic ignores outliers |
| Smooth without blurring boundaries | Bilateral / NLM / TV | range weighting stops at edges |
| Highlight boundaries | Sobel, Laplacian of Gaussian | first/second derivatives |
| Emphasize objects of one size | LoG or band-pass at that scale | scale-selective response |
| Clean a binary mask | Opening / closing with a sized element | size-selective, idempotent |
| Recover some resolution | Wiener / Richardson–Lucy | regularized inverse of a known PSF |
| Preserve quantitative intensities | Do as little as possible | most filters change the intensity distribution |
Checkpoint 5. An investigator applies a \(9\times9\) mean filter and reports that noise fell by a factor of 81. What did they get wrong, and what is the true factor? Answer: they used the window area; noise scales with \(\lVert w\rVert_2 = 1/k\), so the true reduction is a factor of 9, at the cost of an added blur of \(\text{FWHM}\approx k\sqrt{(k^2-1)/12}\cdot 2.355/k \approx 6\) pixels.
Fourier analysis and spatial-frequency filtering
The Fourier transform decomposes an image into spatial-frequency components. Low spatial frequencies represent broad, slowly varying intensity (smooth background, large anatomy); high spatial frequencies represent sharp edges, fine detail, and noise. This dual view is powerful because many image-processing tasks, and much of image formation, especially in MRI and CT, are most naturally understood as operations on spatial-frequency content.
The discrete Fourier transform
For a 1D signal \(y_l\), \(l=0,\dots,N-1\), using Euler’s formula \(e^{i\theta}=\cos\theta+i\sin\theta\), the forward and inverse DFTs are
\[Y_k = \sum_{l=0}^{N-1} y_l\, e^{-i 2\pi k l / N}, \qquad y_l = \frac{1}{N}\sum_{k=0}^{N-1} Y_k\, e^{\,i 2\pi k l / N}.\]
For a 2D image \(f[j,k]\) of size \(N\times M\),
\[F[u,v] = \sum_{j=0}^{N-1}\sum_{k=0}^{M-1} f[j,k]\, e^{-i 2\pi\left(\frac{uj}{N} + \frac{vk}{M}\right)},\]
with the inverse carrying a \(1/(NM)\) factor. (Different packages place the
normalization on the forward transform, the inverse, or split it symmetrically as
\(1/\sqrt{NM}\); R’s fft() puts it entirely on the inverse, which is why every
inverse transform in this chapter is divided by length(x).) The transform is
generally complex-valued, summarized by magnitude and phase:
\[|F[u,v]| = \sqrt{\operatorname{Re}(F)^2 + \operatorname{Im}(F)^2}, \qquad \phi[u,v] = \operatorname{atan2}\!\big(\operatorname{Im}(F),\,\operatorname{Re}(F)\big).\]
Physical frequency units
The DFT index \(u\) is not a frequency. With voxel spacing \(\Delta y\) along rows and \(N\) rows, the physical frequency of index \(u\) (after centring) is
\[f_u = \frac{u - \lfloor N/2\rfloor}{N\,\Delta y}\ \ [\text{cycles mm}^{-1}], \qquad |f|_{\max} = f_{\text{Ny}} = \frac{1}{2\Delta y},\]
and the frequency spacing \(\Delta f = 1/(N\Delta y) = 1/\text{FOV}\). In words:
- the field of view sets the frequency resolution (how finely we resolve frequencies);
- the voxel size sets the maximum frequency (how much detail is representable).
These two statements are the whole of sampling theory as it applies to imaging, and in MRI they are literally the k-space acquisition parameters.
freq_grid <- function(nr, nc, dy = 1, dx = 1) {
fu <- ((0:(nr - 1)) - floor(nr / 2)) / (nr * dy) # rows, cycles/mm
fv <- ((0:(nc - 1)) - floor(nc / 2)) / (nc * dx) # columns, cycles/mm
list(fu = fu, fv = fv, R = sqrt(outer(fu^2, fv^2, "+")))
}
FG <- freq_grid(n, n, dy, dx)
data.frame(FOV_mm = FOV,
voxel_mm = dx,
delta_f_per_mm = round(1 / FOV, 5),
Nyquist_per_mm = 1 / (2 * dx),
max_grid_freq = max(abs(FG$fv)),
DC_index = which(FG$fv == 0))## FOV_mm voxel_mm delta_f_per_mm Nyquist_per_mm max_grid_freq DC_index
## 1 160 1 0.00625 0.5 0.5 81F2 <- fft(img_dn)
mag2 <- Mod(F2)
phase2 <- Arg(F2)
par(mfrow = c(1, 3), mar = c(1, 1, 2.5, 1))
show_img(img_dn, "Image")
show_img(normalize01(log1p(fftshift(mag2))), "log |F| (centerd)")
show_img(fftshift(phase2), "phase (radians)", zlim = c(-pi, pi))
Inverse transform and reconstruction. Recombining magnitude and phase and inverting recovers the image exactly, confirming that the transform is information-preserving:
recon <- Re(fft(mag2 * exp(1i * phase2), inverse = TRUE)) / length(mag2)
c(max_abs_reconstruction_error = max(abs(recon - img_dn)))## max_abs_reconstruction_error
## 9.992e-16Magnitude tells you how much; phase tells you where
A classic experiment settles which half of the transform carries the structure: build a hybrid image from the magnitude of one image and the phase of another.
## a genuinely different second image with a comparably broadband spectrum
make_blobs <- function(seed = 3, nblob = 14) {
Xi <- col(matrix(0, n, n)); Yi <- row(matrix(0, n, n))
b <- matrix(0.10, n, n); set.seed(seed)
for (k in 1:nblob) {
cx <- sample(20:140, 1); cy <- sample(20:140, 1); r <- sample(6:16, 1)
b[(Xi - cx)^2 + (Yi - cy)^2 <= r^2] <- runif(1, 0.40, 0.95)
}
b[abs(Xi - Yi) < 3] <- 0.75
gauss_blur(b, 1.2)
}
imgB <- make_blobs()
FA_ <- fft(img_dn); FB_ <- fft(imgB)
hybrid_magA_phaseB <- Re(fft(Mod(FA_) * exp(1i * Arg(FB_)), inverse = TRUE)) / (n*n)
hybrid_magB_phaseA <- Re(fft(Mod(FB_) * exp(1i * Arg(FA_)), inverse = TRUE)) / (n*n)
par(mfrow = c(1, 4), mar = c(1, 1, 2.5, 1))
show_img(img_dn, "A: phantom")
show_img(imgB, "B: blobs")
show_img(hybrid_magA_phaseB, "|A| with phase of B")
show_img(hybrid_magB_phaseA, "|B| with phase of A")
par(mfrow = c(1, 1))
data.frame(
hybrid = c("|A| + phase(B)", "|B| + phase(A)"),
cor_with_A = round(c(cor(as.vector(hybrid_magA_phaseB), as.vector(img_dn)),
cor(as.vector(hybrid_magB_phaseA), as.vector(img_dn))), 3),
cor_with_B = round(c(cor(as.vector(hybrid_magA_phaseB), as.vector(imgB)),
cor(as.vector(hybrid_magB_phaseA), as.vector(imgB))), 3))## hybrid cor_with_A cor_with_B
## 1 |A| + phase(B) 0.354 0.740
## 2 |B| + phase(A) 0.740 0.355Each hybrid resembles its phase donor, not its magnitude donor. This is why phase errors (motion, \(B_0\) inhomogeneity, eddy currents) are so destructive in MRI, and why phase-based similarity measures are used in registration.
Spatial frequencies in images
| Spatial-frequency component | Image interpretation |
|---|---|
| Low frequency | smooth background, broad anatomy, slow intensity variation, bias field |
| Intermediate | organ boundaries, larger structures, gradual texture |
| High frequency | sharp edges, fine detail, small objects, noise |
Clinical interpretation. A smooth organ region contributes mainly low frequencies; a sharp lesion boundary contributes higher frequencies; random noise is spread across all frequencies (white noise has a flat spectrum). That is exactly why low-pass filtering reduces noise but blurs edges: it discards a band in which signal and noise coexist, and the signal in that band is the edge.
The convolution theorem
The convolution theorem is the bridge between spatial filtering and frequency analysis: convolution in the spatial domain equals multiplication in the frequency domain. If \(g=f*h\) then
\[G(f_x,f_y) = F(f_x,f_y)\,H(f_x,f_y),\]
and, conversely, multiplication in space equals convolution in frequency. Let us verify it numerically rather than take it on faith.
fft_conv <- function(img, kernel) {
nr <- nrow(img); nc <- ncol(img); kr <- nrow(kernel); kc <- ncol(kernel)
NR <- nr + kr - 1; NC <- nc + kc - 1 # zero-pad to avoid wrap-around
A <- matrix(0, NR, NC); A[1:nr, 1:nc] <- img
B <- matrix(0, NR, NC); B[1:kr, 1:kc] <- kernel
full <- Re(fft(fft(A) * fft(B), inverse = TRUE)) / (NR * NC)
pr <- (kr - 1) %/% 2; pc <- (kc - 1) %/% 2
full[(pr + 1):(pr + nr), (pc + 1):(pc + nc)]
}
Kbig <- gk_2d(2.5) # a 15 x 15 Gaussian
c(spatial_vs_fourier_max_difference =
max(abs(conv2(img_dn, Kbig, mode = "zero") - fft_conv(img_dn, Kbig))))## spatial_vs_fourier_max_difference
## 1.332e-15The two agree to machine precision. Beyond correctness, this is the reason FFT convolution is used for large kernels: the cost falls from \(O(N^2k^2)\) to \(O(N^2\log N)\), independent of kernel size.
Key idea. Filtering an image with a kernel is equivalent to multiplying its Fourier transform by the kernel’s Fourier transform. A smoothing kernel multiplies the spectrum by a function that attenuates high frequencies; a sharpening kernel boosts them. “What does this filter do?” is therefore always answerable by plotting \(|H(f)|\).
Point spread function and modulation transfer function
By the convolution theorem, blurring by the PSF becomes multiplication by its transform. The normalized magnitude
\[\text{MTF}(f) = \frac{|H(f)|}{|H(0)|}\]
is the modulation transfer function: the fraction of contrast that survives imaging at each spatial frequency. It is 1 at DC by construction and falls toward 0 at high frequencies. The frequency at which it falls to 10 % (\(f_{10}\)) or 50 % (\(f_{50}\)) is a standard resolution metric. A narrow PSF (sharp system) gives a broad MTF; a broad PSF (blurry system) gives a narrow MTF.
For a Gaussian PSF of width \(\sigma_h\) the MTF is available in closed form,
\[\text{MTF}(f) = \exp\!\big(-2\pi^2\sigma_h^2 f^2\big),\]
so we can check our numerics against theory exactly.
Hgrid <- fft(ifftshift(psf_on_grid(psf_sigma, n)))
MTF <- fftshift(Mod(Hgrid) / Mod(Hgrid)[1, 1])
prof <- MTF[n / 2 + 1, ]
fx <- FG$fv # cycles / mm
theory<- exp(-2 * pi^2 * (psf_sigma * dx)^2 * fx^2)
par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
show_img(MTF, "2D MTF (centerd)", zlim = c(0, 1))
plot(fx[fx >= 0], prof[fx >= 0], type = "l", lwd = 2, col = "steelblue",
xlab = "spatial frequency (cycles / mm)", ylab = "MTF",
main = "MTF: numerical vs analytic")
lines(fx[fx >= 0], theory[fx >= 0], col = "firebrick", lty = 2, lwd = 2)
abline(h = c(0.5, 0.1), lty = 3)
legend("topright", c("numerical", expression(exp(-2*pi^2*sigma^2*f^2))),
col = c("steelblue", "firebrick"), lty = c(1, 2), lwd = 2, bty = "n")
par(mfrow = c(1, 1))
f_at <- function(level) approx(prof[fx >= 0], fx[fx >= 0], xout = level)$y
data.frame(max_abs_MTF_error = signif(max(abs(prof - theory)), 3),
f50_cycles_per_mm = round(f_at(0.5), 3),
f10_cycles_per_mm = round(f_at(0.1), 3),
PSF_FWHM_mm = round(fwhm_of_sigma(psf_sigma) * dx, 2),
Nyquist_per_mm = 1 / (2 * dx))## max_abs_MTF_error f50_cycles_per_mm f10_cycles_per_mm PSF_FWHM_mm
## 1 0.000931 0.156 0.285 2.83
## Nyquist_per_mm
## 1 0.5Reading an MTF. \(f_{10}\approx 0.285\) cycles mm\(^{-1}\) corresponds to a line pair every \(1/(2\times 0.285)\approx 1.8\) mm, and \(f_{50}\approx 0.156\) to one every 3.2 mm, bracketing the 2.8 mm FWHM of the PSF, as they should. Note also that \(f_{10}\) sits far below the Nyquist frequency of 0.5 cycles mm\(^{-1}\): this system is blur-limited, not sampling-limited, so buying a finer matrix would add noise without adding information. When \(f_{10}\gtrsim f_{\text{Ny}}\) the reverse holds and the image is sampling-limited, i.e. it will alias.
Nyquist sampling and aliasing
A discrete image can faithfully represent only frequencies below \(f_{\text{Ny}}=1/(2\Delta x)\): the signal must be sampled at least twice per period. Content above \(f_{\text{Ny}}\) is not lost but folded: a frequency \(f_{\text{Ny}}+\delta\) appears at \(f_{\text{Ny}}-\delta\). That folding is aliasing, demonstrated visually in the inspection section, and it is irreversible, once folded, the true and alias frequencies are indistinguishable.
The engineering consequence is that an anti-alias (low-pass) filter must be applied before sampling or decimation, in hardware where possible (detector aperture, receiver bandwidth) and in software otherwise.
detail_period_px <- 3 # a detail repeating every 3 px
data.frame(period_px = detail_period_px,
frequency_cyc_per_px = round(1 / detail_period_px, 4),
frequency_cyc_per_mm = round(1 / (detail_period_px * dx), 4),
Nyquist_cyc_per_px = 0.5,
adequately_sampled = (1 / detail_period_px) < 0.5,
max_pixel_size_mm = round(detail_period_px * dx / 2, 3))## period_px frequency_cyc_per_px frequency_cyc_per_mm Nyquist_cyc_per_px
## 1 3 0.3333 0.3333 0.5
## adequately_sampled max_pixel_size_mm
## 1 TRUE 1.5Spatial-frequency filtering
Fourier-domain filtering modifies an image by altering selected frequencies. A low-pass filter keeps low frequencies (smooths, denoises); a high-pass filter keeps high frequencies (enhances edges, amplifies noise); a band-pass filter keeps a chosen range (emphasizes a particular size scale); and a notch filter removes a narrow band (useful for periodic stripe/grid artifacts).
Because the filter is a multiplication in frequency, its shape matters as much as its cutoff. A sharp ideal (brick-wall) cutoff corresponds to a sinc kernel in space and therefore rings, the same Gibbs phenomenon seen earlier. Smooth transitions (Gaussian, Butterworth) avoid it.
R <- FG$R # radial frequency, cycles/mm
fc <- 0.08 # cutoff
F0 <- fftshift(fft(img_dn))
apply_mask <- function(Fs, M) Re(fft(ifftshift(Fs * M), inverse = TRUE)) / (n * n)
butter <- function(R, fc, order = 2) 1 / (1 + (R / fc)^(2 * order))
lp_ideal <- apply_mask(F0, R <= fc)
lp_gauss <- apply_mask(F0, exp(-(R / fc)^2 / 2))
lp_butter <- apply_mask(F0, butter(R, fc, 2))
hp_ideal <- apply_mask(F0, R > fc)
bp <- apply_mask(F0, (R > 0.05) & (R < 0.15))
par(mfrow = c(2, 3), mar = c(1, 1, 2.5, 1))
show_img(img_dn, "Original", zlim = c(0, 1))
show_img(lp_ideal, "Ideal low-pass (rings)", zlim = c(0, 1))
show_img(lp_gauss, "Gaussian low-pass", zlim = c(0, 1))
show_img(lp_butter, "Butterworth low-pass n=2", zlim = c(0, 1))
show_img(hp_ideal, "High-pass", zlim = c(-0.2, 0.2))
show_img(bp, "Band-pass 0.05-0.15 /mm", zlim = c(-0.2, 0.2))
r0 <- 40
plot(1:n, img_dn[r0, ], type = "l", lwd = 2, col = "gray50", ylim = c(-0.15, 0.9),
xlab = "column", ylab = "intensity",
main = sprintf("Row %d profile: ideal cutoff rings, smooth cutoffs do not", r0))
lines(lp_ideal[r0, ], col = "firebrick", lwd = 2)
lines(lp_gauss[r0, ], col = "steelblue", lwd = 2)
lines(lp_butter[r0, ], col = "darkgreen", lwd = 2, lty = 2)
legend("topright", c("original", "ideal", "Gaussian", "Butterworth"),
col = c("gray50", "firebrick", "steelblue", "darkgreen"),
lty = c(1, 1, 1, 2), lwd = 2, bty = "n", cex = 0.8)
Ringing is easiest to quantify on a clean synthetic edge, where the true values are exactly 0 and 1 and any excursion outside \([0,1]\) is pure filter artifact.
disk_edge <- ((row(matrix(0, n, n)) - 80)^2 + (col(matrix(0, n, n)) - 80)^2) <= 30^2
Fd <- fftshift(fft(disk_edge * 1))
am <- function(M) Re(fft(ifftshift(Fd * M), inverse = TRUE)) / (n * n)
masks <- list(ideal = (R <= fc),
Gaussian = exp(-(R / fc)^2 / 2),
`Butterworth n=2`= butter(R, fc, 2),
`Butterworth n=8`= butter(R, fc, 8))
t(sapply(masks, function(M) {
z <- am(M)
c(max = round(max(z), 4), overshoot_pct = round(100 * (max(z) - 1), 1),
min = round(min(z), 4), undershoot_pct = round(-100 * min(z), 1))
}))## max overshoot_pct min undershoot_pct
## ideal 1.101 10.1 -0.0953 9.5
## Gaussian 1.000 0.0 0.0000 0.0
## Butterworth n=2 1.041 4.1 -0.0293 2.9
## Butterworth n=8 1.098 9.8 -0.0786 7.9The ideal (brick-wall) filter overshoots by about 10 %, the Gibbs constant again, while the Gaussian produces no overshoot at all. Butterworth filters interpolate between the two as their order increases: sharper frequency selectivity is always paid for with spatial ringing. This is a general fact about linear filters, not a numerical accident.
Notch filtering a periodic artifact
Periodic interference, RF spikes in MRI, grid lines in radiography, scan-line noise in microscopy, appears as isolated bright points in the spectrum and can be removed surgically, without touching the rest of the image.
stripe_art <- 0.10 * sin(2 * pi * 0.18 * col(matrix(0, n, n)))
img_striped <- img_dn + stripe_art
Fs <- fftshift(fft(img_striped))
## build a notch at the offending frequency pair (+/- 0.18 cycles/px = 0.18/dx per mm)
f_art <- 0.18 / dx
notch <- !(abs(abs(outer(rep(0, n), FG$fv, "+")) - f_art) < 0.012 &
abs(outer(FG$fu, rep(0, n), "+")) < 0.012)
img_notched <- apply_mask(Fs, notch)
par(mfrow = c(1, 4), mar = c(1, 1, 2.5, 1))
show_img(img_striped, "With periodic artifact", zlim = c(0, 1))
show_img(normalize01(log1p(Mod(Fs))), "Spectrum: two spikes")
show_img(img_notched, "After notch filter", zlim = c(0, 1))
show_img(img_striped - img_notched, "Removed component", zlim = c(-0.15, 0.15))
par(mfrow = c(1, 1))
c(RMSE_before = round(sqrt(mean((img_striped - img_dn)^2)), 4),
RMSE_after = round(sqrt(mean((img_notched - img_dn)^2)), 4))## RMSE_before RMSE_after
## 0.0709 0.0134Important caution. Frequency filtering can improve visualization but can also remove clinically relevant information or introduce ringing. Notch filtering in particular removes everything at that frequency, including genuine periodic anatomy (trabecular bone, muscle striations). Always compare a filtered image with the original, and look at the removed component, before any quantitative interpretation.
Worked one-dimensional examples: blurring and edge effects
These two classic exercises (after Hobbie & Roth) show, with a single spike, why removing high frequencies blurs an image while removing low frequencies produces edge/ringing effects. Let \(y_0=1\) and \(y_j=0\) for \(j=1,\dots,7\) on a period of \(N=8\). Since the DFT of a unit impulse is flat, \(Y_k=1\) for all \(k\), and the inverse transform can be written as a cosine series with \(a_k=1/8\) for every \(k\) and \(b_k=0\).
Low-pass (blurring). Keep only \(k\in\{0,1,7\}\) (note \(k=7\equiv k=-1\)):
\[y_j = \tfrac{1}{8}\big[\,1 + 2\cos(\pi j/4)\,\big].\]
High-pass (edge effects). Keep only the highest frequencies \(k\in\{3,4,5\}\); using \(\cos(5\pi j/4)=\cos(3\pi j/4)\) for integer \(j\),
\[y_j = \tfrac{1}{8}\big[2\cos(3\pi j/4) + \cos(\pi j)\big].\]
j <- 0:7
y_lp <- (1/8) * (1 + 2 * cos(j * pi / 4))
y_hp <- (1/8) * (2 * cos(3 * pi * j / 4) + cos(pi * j))
## verify against a direct DFT band-selection
Yk <- fft(c(1, rep(0, 7)))
keep <- function(ks) { Z <- rep(0 + 0i, 8); Z[ks + 1] <- Yk[ks + 1]
Re(fft(Z, inverse = TRUE)) / 8 }
c(low_pass_matches = max(abs(keep(c(0, 1, 7)) - y_lp)),
high_pass_matches = max(abs(keep(c(3, 4, 5)) - y_hp)))## low_pass_matches high_pass_matches
## 5.551e-17 2.220e-16par(mfrow = c(1, 2))
plot(j, y_lp, type = "b", pch = 19, col = "steelblue", ylim = c(-0.2, 0.45),
main = "Low-pass: the spike is smeared (blur)", xlab = "j", ylab = expression(y[j]))
abline(h = 0, col = "gray70")
plot(j, y_hp, type = "b", pch = 19, col = "firebrick", ylim = c(-0.2, 0.45),
main = "High-pass: overshoot beside the spike (edges)", xlab = "j", ylab = expression(y[j]))
abline(h = 0, col = "gray70")
The low-pass result spreads the original impulse over its neighbors (blur); the high-pass result is large and negative on either side of the impulse, exaggerating the change, the 1D analogue of edge enhancement and ringing.
Convolution and Fourier filtering: two views of one operation
| Image-space view | Frequency-space view |
|---|---|
| Apply a smoothing kernel | Multiply the spectrum by a low-pass envelope |
| Apply a sharpening kernel | Amplify high frequencies |
| Apply an edge detector | Multiply by a function vanishing at DC |
| Use a broad kernel | Narrow multiplier: stronger high-frequency suppression |
| Use a narrow kernel | Broad multiplier: effect localized in space |
| Cost \(O(N^2k^2)\) | Cost \(O(N^2\log N)\), independent of \(k\) |
Compare. Convolution is local and intuitive in image space; Fourier filtering is global and reveals how a filter reshapes the spatial-frequency content. Use the spatial view to reason about where, the frequency view to reason about what scale. One caveat: FFT filtering assumes the image is periodic, so the left edge wraps into the right. Pad before transforming when this matters.
Reconstruction from projections: the central-slice theorem
Fourier analysis does more than filter, it underlies tomographic reconstruction. CT, SPECT, and PET do not measure \(f(x,y)\) directly; they measure line integrals (projections)
\[p_\theta(s) = \int f\big(s\cos\theta - t\sin\theta,\; s\sin\theta + t\cos\theta\big)\,dt,\]
the Radon transform of \(f\). The central-slice (projection) theorem states that the 1D Fourier transform of \(p_\theta\) equals the radial slice, at angle \(\theta\), of the 2D Fourier transform of \(f\):
\[\mathcal{F}_{1}\{p_\theta\}(\rho) \;=\; F(\rho\cos\theta,\;\rho\sin\theta).\]
Collecting projections over many angles therefore fills the 2D frequency plane, and a 2D inverse transform recovers the image. Because those radial samples are denser near the origin (their density falls as \(1/\rho\)), the correct inversion weights each projection by \(|\rho|\) before back-projecting, the ramp filter of filtered back-projection (FBP):
\[f(x,y) = \int_0^{\pi} \Big(\mathcal{F}_1^{-1}\big\{|\rho|\, \mathcal{F}_1\{p_\theta\}\big\}\Big)(x\cos\theta + y\sin\theta)\;d\theta .\]
In CT the reconstructed quantity \(f\) is the attenuation coefficient \(\mu\); in PET/SPECT it is the radiotracer concentration. Let us build a miniature CT scanner and verify all of this on our own phantom.
## ---- forward projection: rotate, then integrate down the columns ---------
m_ct <- 96
sub <- img_true[round(seq(1, n, length.out = m_ct)),
round(seq(1, n, length.out = m_ct))]
thetas <- seq(0, 179, by = 1.5)
sinogram <- sapply(thetas, function(th) colSums(rotate_image(sub, -th)))
## ---- central-slice theorem, verified ------------------------------------
fftshift1 <- function(v) { N <- length(v); v[c((floor(N/2) + 1):N, 1:floor(N/2))] }
F2d <- fftshift(fft(sub)) # 2D FT of the object
ctr <- m_ct %/% 2 + 1
csl_check <- function(th_deg, slice) {
k <- which.min(abs(thetas - th_deg))
P1 <- fftshift1(fft(sinogram[, k])) # 1D FT of that projection
c(theta = thetas[k],
correlation = cor(Mod(P1), Mod(slice)),
max_relative_difference = max(abs(Mod(P1) - Mod(slice))) / max(Mod(slice)))
}
round(rbind(`theta = 0` = csl_check(0, F2d[ctr, ]),
`theta = 90` = csl_check(90, F2d[, ctr])), 6)## theta correlation max_relative_difference
## theta = 0 0 1 0
## theta = 90 90 1 0## ---- ramp-filtered back-projection ---------------------------------------
ramp_filter <- function(p, window = "hamming") {
N <- length(p); Np <- 2^ceiling(log2(2 * N))
pp <- c(p, rep(0, Np - N))
f <- c(0:(Np / 2), -(Np / 2 - 1):-1) / Np # cycles per sample
H <- abs(f)
if (window == "hamming") H <- H * (0.54 + 0.46 * cos(pi * f / max(abs(f))))
Re(fft(fft(pp) * H, inverse = TRUE))[1:N] / Np
}
back_project <- function(sino, thetas, filtered = TRUE) {
m <- nrow(sino); rec <- matrix(0, m, m)
for (k in seq_along(thetas)) {
q <- if (filtered) ramp_filter(sino[, k]) else sino[, k]
rec <- rec + rotate_image(matrix(rep(q, each = m), m, m), thetas[k])
}
rec * pi / length(thetas)
}
rec_fbp <- back_project(sinogram, thetas, filtered = TRUE)
rec_bp <- back_project(sinogram, thetas, filtered = FALSE)
rec_few <- back_project(sinogram[, seq(1, length(thetas), by = 10)],
thetas[seq(1, length(thetas), by = 10)], TRUE)
par(mfrow = c(2, 3), mar = c(1, 1, 2.5, 1))
show_img(sub, "Object f (down-sampled phantom)", zlim = c(0, 1))
image(x = thetas, y = 1:m_ct, z = t(sinogram), col = gray.colors(256),
xlab = "projection angle (deg)", ylab = "detector bin",
main = "Sinogram (Radon transform)")
show_img(rec_bp, "Unfiltered back-projection (1/r blur)")
show_img(rec_fbp, "Filtered back-projection", zlim = c(0, 1))
show_img(rec_few, "FBP, 12 angles: streaks", zlim = c(0, 1))
show_img(rec_fbp - sub, "FBP error", zlim = c(-0.2, 0.2))
par(mfrow = c(1, 1))
## quantitative agreement inside a central disc (away from the FOV boundary)
cc <- (m_ct + 1) / 2
disc <- (row(matrix(0, m_ct, m_ct)) - cc)^2 + (col(matrix(0, m_ct, m_ct)) - cc)^2 <
(0.42 * m_ct)^2
fitl <- lm(sub[disc] ~ rec_fbp[disc])
data.frame(row.names = NULL,
n_angles = length(thetas),
correlation_with_truth = round(cor(rec_fbp[disc], sub[disc]), 4),
recovered_slope = round(coef(fitl)[2], 3),
correlation_12_angles = round(cor(rec_few[disc], sub[disc]), 4),
correlation_unfiltered = round(cor(rec_bp[disc], sub[disc]), 4))## n_angles correlation_with_truth recovered_slope correlation_12_angles
## 1 120 0.9351 1.092 0.8557
## correlation_unfiltered
## 1 0.4331Three lessons are visible at once. Back-projection without the ramp filter produces the familiar \(1/r\) blur, because it over-counts low frequencies. Filtered back-projection recovers the object with a slope near 1, i.e. it is quantitatively calibrated. And with too few angles the frequency plane is sampled sparsely at large \(\rho\), producing the streak artifacts familiar from low-dose and sparse-view CT, a direct visual manifestation of angular undersampling.
How many views are enough? To sample the outermost frequency ring at the same density as the detector samples each projection, one needs \(N_\theta \gtrsim \frac{\pi}{2}N_{\text{det}}\) views. With 96 detector bins that is about 150 angles; our 120 is close, and 12 is hopeless. The ramp filter also explains why CT reconstruction amplifies noise (\(|\rho|\) grows with frequency), which is why clinical CT applies an apodization window, the Hamming taper used above, trading resolution for noise, exactly as in the Wiener filter.
Review questions. (1) What is convolution in image processing? (2) What is an image kernel? (3) Why does a mean filter reduce noise by \(1/k\) rather than \(1/k^2\)? (4) Why is median filtering not a convolution? (5) What distinguishes low from high spatial frequency, and what sets each limit in a real scan? (6) What does the 2D Fourier transform represent? (7) Why are Fourier transforms complex-valued, and which part carries the spatial structure? (8) State the convolution theorem and give one computational and one conceptual consequence. (9) Why does an ideal low-pass filter ring, and what would you use instead? (10) State the central-slice theorem and explain why FBP needs a ramp filter.
Checkpoint 6. An MRI protocol uses a 256 mm FOV with a \(128\times128\) matrix. What are the voxel size, the Nyquist frequency, and the frequency spacing? If the object contains structure at 0.3 cycles mm\(^{-1}\), what happens? Answers: 2 mm; \(f_{\text{Ny}} = 0.25\) cycles mm\(^{-1}\); \(\Delta f = 1/256\) cycles mm\(^{-1}\); the 0.3 cycles mm\(^{-1}\) content is above Nyquist and folds back to \(0.25-(0.3-0.25)=0.2\) cycles mm\(^{-1}\), appearing as a coarser, false pattern.
Registration and spatial alignment
Registration finds the spatial transformation \(T\) that brings a moving image into correspondence with a fixed (reference) image. It is essential for longitudinal studies, multimodal imaging, atlas-based analysis, motion correction, and radiation-treatment planning.
Every registration algorithm is the same optimization problem:
\[\hat T = \arg\min_{T}\; \underbrace{\mathcal{D}\big(I_{\text{fix}},\, I_{\text{mov}}\circ T\big)}_{\text{dissimilarity}} \;+\;\lambda\,\underbrace{\mathcal{R}(T)}_{\text{regularization}} ,\]
so specifying a registration means specifying four things: a transformation model (what \(T\) may be), a similarity metric (\(\mathcal D\)), an interpolator (how \(I_{\text{mov}}\circ T\) is evaluated off-grid), and an optimizer (how the search is conducted).
| Model | Free parameters (2D / 3D) | Preserves | Example use |
|---|---|---|---|
| Translation | 2 / 3 | everything but position | motion correction of a rigid organ |
| Rigid | 3 / 6 | distances and angles | repeat scans of the same head |
| Similarity | 4 / 7 | angles and shape | cross-scanner scaling |
| Affine | 6 / 12 | parallel lines | global size/shape normalization |
| B-spline / deformable | \(10^3\)–\(10^6\) | local topology (if regularized) | breathing, tumour growth, inter-subject |
Rigid transformations preserve distances and shapes. Affine adds global scaling and shear. Deformable registration allows local shape change and is needed when anatomy deforms or differs between subjects. Multimodal registration aligns images whose intensities have different physical meanings, so it relies on metrics that do not assume matching intensities.
Similarity metrics
The optimizer needs a number that says how well two images agree. Three are fundamental, and choosing the wrong one is the commonest cause of registration failure.
Sum of squared differences (SSD). Minimal when intensities match voxel-for-voxel. Optimal for same-modality, same-contrast images with additive Gaussian noise, and invalid whenever a global intensity change is present:
\[\text{SSD} = \frac{1}{|\Omega|}\sum_{x\in\Omega} \big(I_{\text{fix}}(x) - I_{\text{mov}}(T(x))\big)^2 .\]
Normalized cross-correlation (NCC). Invariant to affine intensity changes \(I\mapsto aI+b\), so it tolerates brightness and gain differences within a modality:
\[\text{NCC} = \frac{\sum_x (I_{\text{fix}}-\bar I_{\text{fix}}) (I_{\text{mov}}\!\circ\! T-\bar I_{\text{mov}})} {\sqrt{\sum_x (I_{\text{fix}}-\bar I_{\text{fix}})^2}\; \sqrt{\sum_x (I_{\text{mov}}\!\circ\! T-\bar I_{\text{mov}})^2}} \in[-1,1].\]
Mutual information (MI). Measures statistical dependence between the two intensity distributions and assumes only that the relationship is consistent, not that it is linear or even monotone. This makes it the standard metric for multimodal registration (CT–MRI, PET–MRI):
\[\text{MI}(A,B) = \sum_{a,b} p(a,b)\,\log_2\frac{p(a,b)}{p(a)\,p(b)} = H(A) + H(B) - H(A,B),\]
with \(H\) the Shannon entropy. Because MI grows when the overlap region shrinks (fewer voxels, lower joint entropy), practical implementations use normalized mutual information, \(\text{NMI}=\big(H(A)+H(B)\big)/H(A,B)\), which is far more robust to changing overlap.
ssd <- function(a, b) mean((a - b)^2)
ncc <- function(a, b) cor(as.vector(a), as.vector(b))
joint_hist <- function(a, b, bins = 32, range_a = range(a), range_b = range(b)) {
ia <- pmin(bins, pmax(1, floor((as.vector(a) - range_a[1]) /
(diff(range_a) + 1e-12) * bins) + 1))
ib <- pmin(bins, pmax(1, floor((as.vector(b) - range_b[1]) /
(diff(range_b) + 1e-12) * bins) + 1))
matrix(tabulate((ib - 1) * bins + ia, nbins = bins * bins), bins, bins)
}
entropy_bits <- function(p) { p <- p[p > 0]; -sum(p * log2(p)) }
mutual_information <- function(a, b, bins = 32, normalized = FALSE) {
J <- joint_hist(a, b, bins); P <- J / sum(J)
pa <- rowSums(P); pb <- colSums(P)
Ha <- entropy_bits(pa); Hb <- entropy_bits(pb); Hab <- entropy_bits(P)
if (normalized) (Ha + Hb) / Hab else Ha + Hb - Hab
}Why mutual information is required for multimodal data
We simulate a second modality whose intensity mapping is non-monotone with respect to the first: the surrounding medium is bright, the parenchyma dark, and the tumour brightest of all. This is exactly the situation between CT and a \(T_2\)-weighted or fluid-suppressed MRI, or between a transmission micrograph and a fluorescence channel, no single monotone curve maps one intensity scale onto the other.
## "modality B": a non-monotone remap of the same anatomy
modB <- matrix(0.60, n, n) # bright surrounding medium
modB[brain_mask] <- 0.25 # dark parenchyma
modB[vessel_mask] <- 0.75
modB[lesion_mask] <- 0.45
modB[tumor_mask] <- 0.85 # tumour still brightest
modB <- gauss_blur(modB, psf_sigma) + matrix(rnorm(n * n, 0, sigma_noise), n, n)
offsets <- -20:20
curves <- sapply(offsets, function(o) {
cand <- translate_image(modB, dx = o, dy = 0)
c(SSD = ssd(img_dn, cand), NCC = ncc(img_dn, cand),
MI = mutual_information(img_dn, cand),
NMI = mutual_information(img_dn, cand, normalized = TRUE))
})
par(mfrow = c(2, 3), mar = c(4, 4, 3, 1))
show_img(img_dn, "Modality A (fixed)", zlim = c(0, 1))
show_img(modB, "Modality B (moving)", zlim = c(0, 1))
JH <- joint_hist(img_dn, modB, bins = 48)
image(x = seq(0, 1, length.out = 48), y = seq(0, 1, length.out = 48),
z = log1p(JH), col = hcl.colors(64, "YlOrRd", rev = TRUE),
xlab = "A intensity (scaled)", ylab = "B intensity (scaled)",
main = "Joint histogram at alignment")
plot(offsets, curves["SSD", ], type = "b", pch = 19, col = "firebrick",
xlab = "horizontal offset (px)", ylab = "SSD", main = "SSD (minimize)")
abline(v = 0, lty = 2)
plot(offsets, curves["NCC", ], type = "b", pch = 19, col = "steelblue",
xlab = "horizontal offset (px)", ylab = "NCC", main = "NCC (maximize)")
abline(v = 0, lty = 2)
plot(offsets, curves["NMI", ], type = "b", pch = 19, col = "darkgreen",
xlab = "horizontal offset (px)", ylab = "NMI", main = "NMI (maximize)")
abline(v = 0, lty = 2)
par(mfrow = c(1, 1))
data.frame(metric = c("SSD (min)", "NCC (max)", "MI (max)", "NMI (max)"),
optimum_offset = c(offsets[which.min(curves["SSD", ])],
offsets[which.max(curves["NCC", ])],
offsets[which.max(curves["MI", ])],
offsets[which.max(curves["NMI", ])]),
correct = c(offsets[which.min(curves["SSD", ])],
offsets[which.max(curves["NCC", ])],
offsets[which.max(curves["MI", ])],
offsets[which.max(curves["NMI", ])]) == 0)## metric optimum_offset correct
## 1 SSD (min) 20 FALSE
## 2 NCC (max) 20 FALSE
## 3 MI (max) 0 TRUE
## 4 NMI (max) 0 TRUESSD and NCC are misled: both are optimized at the edge of the search range rather than at zero offset, because at true alignment bright tissue in A faces dark tissue in B and the correlation is strongly negative (NCC \(\approx -0.65\)). Mutual information and NMI locate the true optimum, because they require only that the mapping between intensities be statistically consistent, not linear, not even monotone. Note also that taking \(|\)NCC\(|\) would rescue this particular example but fails as soon as the relationship is non-monotone within one image, which is the general multimodal case.
Rigid registration in practice
We now register a shifted and rotated copy of our own image, so the answer is known exactly.
fixed <- img_dn
true_dx <- 9; true_dy <- -6
moving <- translate_image(fixed, dx = true_dx, dy = true_dy)
search <- expand.grid(dx = -15:15, dy = -15:15)
score <- apply(search, 1, function(p)
ncc(fixed, translate_image(moving, dx = -p[1], dy = -p[2])))
best <- search[which.max(score), ]
registered <- translate_image(moving, dx = -best$dx, dy = -best$dy)
par(mfrow = c(1, 4), mar = c(1, 1, 2.5, 1))
show_img(fixed, "Fixed", zlim = c(0, 1))
show_img(moving, sprintf("Moving (dx=%d, dy=%d)", true_dx, true_dy), zlim = c(0, 1))
show_img(registered, sprintf("Registered (recovered %d, %d)", best$dx, best$dy), zlim = c(0, 1))
show_img(fixed - registered, "Residual difference", zlim = c(-0.3, 0.3))
true_angle <- 17
moving_rot <- rotate_image(fixed, true_angle)
angles <- seq(-30, 30, by = 2)
rscore <- sapply(angles, function(a) ncc(fixed, rotate_image(moving_rot, -a)))
## sub-sample-accuracy refinement: fit a parabola through the peak and its neighbors
k <- which.max(rscore)
trio <- rscore[(k - 1):(k + 1)]
delta<- 0.5 * (trio[1] - trio[3]) / (trio[1] - 2 * trio[2] + trio[3])
angle_refined <- angles[k] + delta * diff(angles)[1]
par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
plot(angles, rscore, type = "b", pch = 19, col = "darkgreen",
xlab = "candidate correction angle (degrees)", ylab = "NCC with fixed image",
main = "Similarity vs. rotation")
abline(v = true_angle, col = "red", lty = 2)
points(angle_refined, max(rscore), col = "blue", pch = 8, cex = 1.5)
show_img(rotate_image(moving_rot, -angle_refined),
sprintf("Registered at %.2f deg", angle_refined), zlim = c(0, 1))
par(mfrow = c(1, 1))
data.frame(true_angle = true_angle,
grid_optimum = angles[k],
grid_step = diff(angles)[1],
parabolic_refinement = round(angle_refined, 3),
residual_error_deg = round(angle_refined - true_angle, 3))## true_angle grid_optimum grid_step parabolic_refinement residual_error_deg
## 1 17 16 2 16.95 -0.052Sub-voxel accuracy comes from interpolating the metric, not the image. A coarse parameter grid plus a quadratic fit through the peak recovers the transformation to a small fraction of the grid step, here to \(\sim0.05^{\circ}\) from a \(2^{\circ}\) search. Real optimizers (Powell, gradient descent, LBFGS) do the same thing continuously.
The optimization landscape: local minima and multiresolution
Exhaustive search is only possible in one or two dimensions. Real registrations have 6–12 (rigid/affine) or thousands (deformable) of parameters and must use iterative optimization, which can be trapped by local optima, especially in images with repeating structure or a large initial misalignment.
The standard remedy is a multiresolution (pyramid) strategy: register heavily smoothed, down-sampled versions first, where the landscape is smooth and the capture range is wide, then refine at successively finer scales.
To make the difficulty visible we add a periodic texture, ribs, trabecular bone, a scanner grid, muscle striations, to the phantom. Periodic structure creates spurious similarity peaks at multiples of its period, which is precisely when local optimizers fail.
texture <- 0.25 * sin(2 * pi * col(matrix(0, n, n)) / 11) * brain_mask
A_tex <- img_dn + texture
B_tex <- translate_image(A_tex, dx = 26, dy = 0) # true shift = 26 px
offs <- -40:40
scales<- c(0, 1, 3, 6)
land <- sapply(scales, function(s) {
A <- if (s == 0) A_tex else gauss_blur(A_tex, s)
B <- if (s == 0) B_tex else gauss_blur(B_tex, s)
sapply(offs, function(o) ncc(A, translate_image(B, dx = -o, dy = 0)))
})
n_local_max <- function(v) sum(diff(sign(diff(v))) == -2)
matplot(offs, land, type = "l", lwd = 2, lty = 1,
col = c("firebrick", "darkorange", "seagreen", "steelblue"),
xlab = "candidate horizontal correction (px)", ylab = "NCC",
main = "Smoothing removes spurious optima and widens the basin")
abline(v = 26, lty = 2)
legend("topleft", sprintf("sigma = %g", scales), bty = "n", lwd = 2,
col = c("firebrick", "darkorange", "seagreen", "steelblue"))
data.frame(smoothing_sigma = scales,
local_maxima = apply(land, 2, n_local_max),
argmax_offset = offs[apply(land, 2, which.max)])## smoothing_sigma local_maxima argmax_offset
## 1 0 7 26
## 2 1 7 26
## 3 3 2 25
## 4 6 1 23A local optimizer makes the point sharply. A simple hill-climb started 12 pixels away stalls on a spurious peak at full resolution, but a coarse-to-fine schedule walks all the way to the truth:
hill_climb <- function(A, B, start) {
x <- start
repeat {
cand <- c(x - 1, x, x + 1)
sc <- sapply(cand, function(o) ncc(A, translate_image(B, dx = -o, dy = 0)))
if (which.max(sc) == 2) return(x)
x <- cand[which.max(sc)]
}
}
naive <- hill_climb(A_tex, B_tex, 12)
step1 <- hill_climb(gauss_blur(A_tex, 6), gauss_blur(B_tex, 6), 12)
step2 <- hill_climb(gauss_blur(A_tex, 2), gauss_blur(B_tex, 2), step1)
step3 <- hill_climb(A_tex, B_tex, step2)
data.frame(strategy = c("single-scale hill climb", "coarse-to-fine pyramid"),
start = c(12, 12),
result = c(naive, step3),
truth = 26,
path = c(paste(12, naive, sep = " -> "),
paste(12, step1, step2, step3, sep = " -> ")))## strategy start result truth path
## 1 single-scale hill climb 12 15 26 12 -> 15
## 2 coarse-to-fine pyramid 12 26 26 12 -> 23 -> 26 -> 26Note also that the coarse estimate is slightly biased (smoothing shifts the optimum by a pixel or two). A pyramid must therefore always end at full resolution; the coarse levels supply a basin of attraction, not the answer.
Interpolation, resampling, and masked registration
Because the moving image rarely lands exactly on the fixed grid, registration interpolates: linear or spline for intensities, nearest-neighbor for label maps. To keep measurements faithful, compose all transforms and resample once — we measured earlier what repeated resampling costs.
When only part of the image should drive alignment, a surgical cavity, the scanner table, air, a large tumour that has changed, a moved fiducial marker — masked registration restricts the similarity metric to a chosen region.
Two caveats are worth stating plainly, because masking is often applied reflexively. First, for a low-dimensional, high-contrast problem such as the rigid alignment above, the dominant cue is the object outline and the optimum is usually recovered with or without a mask; masking earns its keep in deformable registration, where a contaminated region locally drags the warp, and in low-contrast or low-overlap problems. Second, masking shrinks the sample: a mask that removes most of the informative voxels can make the metric noisier and the optimum less well determined than leaving the contamination in. Masking is a bias–variance trade-off, and like every other choice in the pipeline it should be reported.
A related trap is overlap dependence. Plain mutual information tends to increase as the overlap region shrinks, because fewer voxels means lower joint entropy, so an optimizer maximizing MI can be rewarded for sliding the images apart. This is exactly what NMI is designed to prevent, and it is why NMI, not MI, is the default in production registration packages.
Deformable registration and the Jacobian
A deformable transform is a displacement field \(u(x)\) with \(T(x)=x+u(x)\). Its local behavior is described by the Jacobian matrix \(J = \partial T/\partial x = I + \nabla u\), whose determinant has a direct physical reading:
- \(\det J > 1\) : local expansion (tissue growth, inflation);
- \(\det J = 1\) : volume preserved;
- \(0<\det J<1\) : local contraction (atrophy, deflation);
- \(\det J \le 0\) : the mapping has folded, anatomically impossible and a sign that regularization was too weak.
This is the mathematical foundation of tensor-based morphometry, where \(\log\det J\) is used as a voxel-wise measure of local volume change (e.g. regional brain atrophy).
## a smooth, known deformation: local expansion around the tumour center
amp <- 6; wid <- 22
c_r <- round(mean(which(tumor_mask, arr.ind = TRUE)[, 1]))
c_c <- round(mean(which(tumor_mask, arr.ind = TRUE)[, 2]))
gaussb <- exp(-((row(matrix(0, n, n)) - c_r)^2 +
(col(matrix(0, n, n)) - c_c)^2) / (2 * wid^2))
ux <- amp * gaussb * (col(matrix(0, n, n)) - c_c) / wid # displacement, columns
uy <- amp * gaussb * (row(matrix(0, n, n)) - c_r) / wid # displacement, rows
warped <- warp_bilinear(fixed, .row_idx(n, n) + uy, .col_idx(n, n) + ux)
## finite-difference Jacobian determinant of T(x) = x + u(x)
d <- function(M, axis) if (axis == 1)
(rbind(M[-1, ], M[n, ]) - rbind(M[1, ], M[-n, ])) / 2 else
(cbind(M[, -1], M[, n]) - cbind(M[, 1], M[, -n])) / 2
J11 <- 1 + d(uy, 1); J12 <- d(uy, 2)
J21 <- d(ux, 1); J22 <- 1 + d(ux, 2)
detJ <- J11 * J22 - J12 * J21
par(mfrow = c(1, 4), mar = c(1, 1, 2.5, 1))
show_img(fixed, "Before deformation", zlim = c(0, 1))
show_img(warped, "After deformation", zlim = c(0, 1))
show_img(sqrt(ux^2+uy^2),"Displacement magnitude (px)")
show_img(detJ, "det J (>1 expand, <1 shrink)", zlim = c(0.5, 1.6),
col = hcl.colors(256, "Blue-Red"))
par(mfrow = c(1, 1))
data.frame(min_detJ = round(min(detJ), 3), max_detJ = round(max(detJ), 3),
folded_voxels = sum(detJ <= 0),
implied_volume_change_pct = round(100 * (sum(detJ) / (n * n) - 1), 3))## min_detJ max_detJ folded_voxels implied_volume_change_pct
## 1 0.922 1.619 0 0.373Deformable registration can create the finding. A sufficiently flexible warp can match any two images, including matching a tumour onto healthy tissue. Always report the transformation model, the regularization (\(\mathcal{R}(T)\), typically bending energy \(\int\lVert\nabla^2 u\rVert^2\)), the number of degrees of freedom, and the fraction of voxels with \(\det J\le 0\). A registration with no negative Jacobians and a physically plausible \(\det J\) range is evidence; a Dice score alone is not.
Registration quality checks
Numerical convergence is not alignment. Report a geometric error where possible: the fiducial registration error (FRE) is the residual at the landmarks used to fit the transform, while the target registration error (TRE) is the error at independent points of clinical interest, and TRE is the one that matters, because FRE can be driven to zero by overfitting. Crucially, for a rotational residual \(\delta\theta\) the TRE at a point \(r\) from the center of rotation is \(r\,\delta\theta\): a registration that is excellent at the center can be poor at the periphery, so TRE must be reported where the measurement is made, not where the optimizer looked.
checkerboard <- function(A, B, k = 20) {
tile <- ((row(A) - 1) %/% k + (col(A) - 1) %/% k) %% 2 == 0
ifelse(tile, A, B)
}
misaligned <- translate_image(fixed, dx = 4, dy = 3)
par(mfrow = c(1, 3), mar = c(1, 1, 2.5, 1))
show_img(checkerboard(fixed, registered), "Checkerboard: well registered", zlim = c(0, 1))
show_img(checkerboard(fixed, misaligned), "Checkerboard: 5 px error", zlim = c(0, 1))
show_img(fixed - misaligned, "Difference image", zlim = c(-0.3, 0.3))
par(mfrow = c(1, 1))
## TRE implied by the residual rotation error recovered earlier.
## For a rotation error d(theta), the displacement at a point r mm from the
## center of rotation is r * d(theta) -- so TRE grows linearly with distance.
theta_err <- (angle_refined - true_angle) * pi / 180
centroid <- function(m) colMeans(which(m, arr.ind = TRUE))
cen <- c((n + 1) / 2, (n + 1) / 2)
lms <- rbind(tumour = centroid(tumor_mask),
lesion = centroid(lesion_mask),
`FOV corner` = c(10, 10))
radius_mm <- sqrt(((lms[, 1] - cen[1]) * dy)^2 + ((lms[, 2] - cen[2]) * dx)^2)
data.frame(landmark = rownames(lms),
radius_mm = round(radius_mm, 1),
TRE_mm = round(radius_mm * abs(theta_err), 3),
row.names = NULL)## landmark radius_mm TRE_mm
## 1 tumour 38.2 0.034
## 2 lesion 42.9 0.039
## 3 FOV corner 99.7 0.090## angle_error_deg
## -0.0516Check alignment, don’t assume it. Overlay edges of the registered moving image on the fixed image; inspect a checkerboard and a difference image; verify that independent landmarks coincide and report the TRE in millimeters; state the final metric value and the metric used; and confirm the result is anatomically plausible, not merely numerically optimal. For deformable registration, additionally report the Jacobian statistics.
Multimodal image fusion
Fusion overlays complementary information after registration, CT–MRI, PET–CT, PET–MRI, or functional maps on anatomical MRI. A typical workflow: preprocess each image, register the moving image to the fixed image, resample the moving image into the fixed space once, overlay or extract matched-region measurements, and validate anatomical plausibility.
Fusion is powerful (CT supplies high-resolution anatomy and attenuation correction, MRI supplies soft-tissue contrast, PET supplies molecular information) but can mislead: a misregistration of one PET voxel is several millimeters of anatomy, and resolution mismatches mean the fused display invites the eye to attribute PET signal to CT structures far smaller than the PET PSF.
## a low-resolution "functional" map fused onto the high-resolution anatomy
func <- gauss_blur(1.0 * tumor_mask + 0.35 * lesion_mask, 4) +
matrix(rnorm(n * n, 0, 0.01), n, n)
fuse <- function(anat, funct, alpha = 0.55, thresh = 0.12) {
a <- normalize01(anat); f <- normalize01(funct)
m <- f > thresh
r <- a; g <- a; b <- a
r[m] <- (1 - alpha) * a[m] + alpha * f[m]
g[m] <- (1 - alpha) * a[m] + alpha * (1 - f[m]) * 0.6
b[m] <- (1 - alpha) * a[m] * 0.4
arr <- array(c(r, g, b), dim = c(nrow(anat), ncol(anat), 3))
plot(as.raster(arr))
}
par(mfrow = c(1, 3), mar = c(1, 1, 2.5, 1))
show_img(img_dn, "Anatomy (high resolution)", zlim = c(0, 1))
show_img(func, "Function (low resolution)")
fuse(img_dn, func); title("Fused overlay")
Checkpoint 7. You must align a baseline and a 12-month follow-up brain MRI of the same subject to measure hippocampal atrophy. Which transformation model, which metric, and which single design decision most threatens the validity of the result? Answer: rigid (6 DOF), the skull does not deform; NCC or SSD (same modality); and the threat is asymmetric processing, resampling only the follow-up image biases it by one interpolation relative to the baseline. Register both to a subject-specific midpoint space, or use an inverse-consistent method.
Segmentation: from intensities to labeled structures
Segmentation separates structures of interest from the rest of the image — organs, tumours, lesions, vessels, tissue compartments, or single cells. It defines what will be measured, so it is the single most consequential step for quantitative accuracy.
Key idea. Segmentation defines what gets measured. A poor segmentation produces poor measurements even when the measurement formula is exactly correct — and, unlike noise, segmentation error does not average away over repeats.
Common methods span thresholding, region growing, clustering, edge detection, active contours and level sets, watershed, atlas-based segmentation, and deep-learning approaches such as U-Net. By degree of automation they range from manual through semi-automatic to fully automatic and AI-assisted.
Thresholding
Thresholding classifies pixels by whether their value is above, below, or within a range:
\[M(x,y) = \begin{cases} 1, & I(x,y) \ge T\\ 0, & I(x,y) < T. \end{cases}\]
The operation is general; the meaning of the threshold is modality-specific (HU in CT, sequence-dependent intensity in MRI, SUV in PET, echo intensity in ultrasound, fluorescence intensity in microscopy).
Global thresholding and its failure mode
A single threshold for the whole image works when the histogram is cleanly separated and the intensity scale is spatially uniform, which is exactly what a bias field destroys.
T_global <- 0.62
par(mfrow = c(1, 4), mar = c(4, 4, 2.5, 1))
show_img(img_dn, "Denoised, bias-corrected", zlim = c(0, 1))
hist(img_dn[head_mask], breaks = 60, col = "gray80",
main = "Histogram inside the head", xlab = "intensity")
abline(v = T_global, col = "red", lwd = 2)
show_img(img_dn > T_global, "Threshold on corrected image", col = c("black", "white"))
show_img(img_noisy > T_global, "Same threshold, UNcorrected", col = c("black", "white"))
Adaptive (local) thresholding
Adaptive thresholding compares each pixel to a local background estimate, so the effective threshold varies across the image. It rescues a global method under shading, at the cost of introducing a window-size parameter and of responding to local contrast rather than absolute value, which makes it unsuited to calibrated modalities where the absolute value is the point.
local_bg <- box_blur(img_noisy, 31)
mask_adaptive <- img_noisy > (local_bg + 0.08)
par(mfrow = c(1, 3), mar = c(1, 1, 2.5, 1))
show_img(img_noisy, "Uncorrected input", zlim = c(0, 1))
show_img(local_bg, "Local background (31 px)", zlim = c(0, 1))
show_img(mask_adaptive, "Adaptive threshold", col = c("black", "white"))
Otsu’s method
When the histogram is bimodal, Otsu’s method chooses the threshold that minimizes the within-class variance. Writing \(\omega_0(T),\omega_1(T)\) for the class probabilities and \(\mu_0,\mu_1\) for the class means, the total variance decomposes as
\[\sigma_{\text{total}}^2 = \underbrace{\omega_0\sigma_0^2+\omega_1\sigma_1^2}_{\text{within}} + \underbrace{\omega_0\omega_1(\mu_0-\mu_1)^2}_{\text{between}},\]
and since \(\sigma^2_{\text{total}}\) does not depend on \(T\), minimizing the within-class variance is equivalent to maximizing the between-class variance
\[\sigma_B^2(T) = \omega_0(T)\,\omega_1(T)\,\big(\mu_0(T)-\mu_1(T)\big)^2 = \frac{\big(\mu_T\,\omega_0(T) - \mu(T)\big)^2}{\omega_0(T)\big(1-\omega_0(T)\big)},\]
which can be evaluated for every candidate \(T\) in one pass over the cumulative
histogram. Our otsu_threshold() (defined earlier) does exactly this.
otsu_curve <- function(x, nbins = 128) {
v <- as.vector(x); br <- seq(min(v), max(v), length.out = nbins + 1)
h <- hist(v, breaks = br, plot = FALSE); p <- h$counts / sum(h$counts)
w0 <- cumsum(p); mu <- cumsum(p * h$mids); muT <- mu[nbins]
sb <- (muT * w0 - mu)^2 / (w0 * (1 - w0)); sb[!is.finite(sb)] <- NA
data.frame(threshold = h$mids, between_class_variance = sb)
}
oc <- otsu_curve(img_noisy)
T_o <- otsu_threshold(img_noisy)
par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
plot(oc, type = "l", lwd = 2, col = "steelblue",
main = "Otsu: maximize between-class variance", xlab = "threshold")
abline(v = T_o, col = "red", lty = 2)
show_img(img_noisy > T_o, sprintf("Otsu mask (T = %.3f)", T_o), col = c("black", "white"))
Multi-level Otsu
Real tissue is rarely two classes. Extending the same criterion to \(k\) classes — maximize \(\sum_c \omega_c(\mu_c-\mu_T)^2\) over all \(k-1\) cut points, separates parenchyma, intermediate structures, and the tumour in one step.
multi_otsu <- function(x, nbins = 128) {
v <- as.vector(x); v <- v[is.finite(v)]
br <- seq(min(v), max(v), length.out = nbins + 1)
h <- hist(v, breaks = br, plot = FALSE)
p <- h$counts / sum(h$counts); mids <- h$mids
w <- cumsum(p); mu <- cumsum(p * mids); muT <- mu[nbins]
best <- -Inf; cuts <- c(NA, NA)
for (a in 2:(nbins - 2)) for (b in (a + 1):(nbins - 1)) {
w1 <- w[a]; w2 <- w[b] - w[a]; w3 <- 1 - w[b]
if (w1 <= 0 || w2 <= 0 || w3 <= 0) next
m1 <- mu[a] / w1; m2 <- (mu[b] - mu[a]) / w2; m3 <- (muT - mu[b]) / w3
s <- w1 * (m1 - muT)^2 + w2 * (m2 - muT)^2 + w3 * (m3 - muT)^2
if (s > best) { best <- s; cuts <- c(mids[a], mids[b]) }
}
cuts
}
cuts <- multi_otsu(img_dn[head_mask])
data.frame(cut_1 = round(cuts[1], 3), cut_2 = round(cuts[2], 3),
true_parenchyma = 0.40, true_vessel = 0.62, true_tumour = 0.80)## cut_1 cut_2 true_parenchyma true_vessel true_tumour
## 1 0.412 0.59 0.4 0.62 0.8seg <- matrix(0L, n, n)
seg[head_mask] <- 1L
seg[head_mask & img_dn > cuts[1]] <- 2L
seg[head_mask & img_dn > cuts[2]] <- 3L
par(mfrow = c(1, 2), mar = c(1, 1, 2.5, 1))
show_img(img_dn, "Denoised input", zlim = c(0, 1))
show_img(seg, "Three-class Otsu labels",
col = c("black", "steelblue", "orange", "firebrick"), zlim = c(0, 3))
par(mfrow = c(1, 1))
table(label = c("background", "parenchyma", "intermediate", "bright")[seg + 1L])## label
## background bright intermediate parenchyma
## 10253 541 4541 10265The ROC view: what is the best possible threshold?
Because we have ground truth, we can ask what any intensity threshold could achieve. The ROC curve traces sensitivity against 1 − specificity over all thresholds; the area under the curve (AUC) summarizes separability, and Youden’s index (\(J=\text{sens}+\text{spec}-1\)) identifies the best operating point. This also quantifies exactly how much bias correction was worth.
roc_curve <- function(score, truth) {
o <- order(score, decreasing = TRUE)
tt <- truth[o]
data.frame(threshold = score[o],
tpr = cumsum(tt) / sum(tt),
fpr = cumsum(!tt) / sum(!tt))
}
auc_of <- function(r)
sum(diff(c(0, r$fpr)) * (r$tpr + c(0, head(r$tpr, -1))) / 2)
abnormal <- tumor_mask | vessel_mask | lesion_mask
r_raw <- roc_curve(img_noisy[head_mask], abnormal[head_mask])
r_bc <- roc_curve(img_dn[head_mask], abnormal[head_mask])
plot(r_raw$fpr, r_raw$tpr, type = "l", lwd = 2, col = "firebrick",
xlab = "1 - specificity", ylab = "sensitivity",
main = "Detecting abnormal tissue by intensity alone")
lines(r_bc$fpr, r_bc$tpr, lwd = 2, col = "steelblue")
abline(0, 1, lty = 3)
legend("bottomright", bty = "n", lwd = 2, col = c("firebrick", "steelblue"),
legend = c(sprintf("uncorrected AUC = %.3f", auc_of(r_raw)),
sprintf("bias-corrected AUC = %.3f", auc_of(r_bc))))
youden <- function(r) { j <- r$tpr - r$fpr; k <- which.max(j)
c(threshold = r$threshold[k], sensitivity = r$tpr[k], specificity = 1 - r$fpr[k]) }
round(rbind(uncorrected = youden(r_raw), corrected = youden(r_bc)), 3)## threshold sensitivity specificity
## uncorrected 0.475 0.886 0.908
## corrected 0.454 0.993 0.974Bias correction raises the AUC substantially: the same threshold operation becomes a better classifier once the nuisance field is removed. This is the cleanest possible demonstration that preprocessing is not cosmetic.
Binary masks and regions of interest
A mask is a (usually binary) image marking the pixels of a region of interest (ROI):
\[M(x,y,z) = \begin{cases} 1, & \text{voxel belongs to the region}\\ 0, & \text{otherwise.}\end{cases}\]
Masks matter because they define exactly which voxels enter later measurements. Masks should be saved and reported alongside the original image, its geometry, and the preprocessing steps, because together they record what was measured.
Cleaning binary masks
Raw threshold masks contain specks (false positives) and holes (false negatives). Morphological opening removes specks; closing and hole-filling repair interiors. Cleaning a mask before measurement materially changes the result, so it belongs in the reported pipeline.
Connected components and object labeling
A binary mask marks foreground pixels but does not distinguish one object from another. Connected-component labeling groups neighboring foreground pixels and assigns each group a unique integer, the basis for counting objects, removing small artifacts, and measuring objects individually. Connectivity must be stated: in 2D, 4-connectivity (edge neighbors) or 8-connectivity (including diagonals); in 3D, 6-, 18-, or 26-connectivity. The choice changes object counts, especially for thin structures.
## Teaching implementation: explicit flood fill, 4-connectivity.
label_components <- function(mask) {
nr <- nrow(mask); nc <- ncol(mask)
lab <- matrix(0L, nr, nc); cur <- 0L
for (s in which(mask)) {
if (lab[s] != 0L) next
cur <- cur + 1L; stack <- s; lab[s] <- cur
while (length(stack)) {
p <- stack[length(stack)]; stack <- stack[-length(stack)]
i <- ((p - 1L) %% nr) + 1L
j <- ((p - 1L) %/% nr) + 1L
for (d in list(c(-1L, 0L), c(1L, 0L), c(0L, -1L), c(0L, 1L))) {
ni <- i + d[1]; nj <- j + d[2]
if (ni >= 1 && ni <= nr && nj >= 1 && nj <= nc &&
mask[ni, nj] && lab[ni, nj] == 0L) {
lab[ni, nj] <- cur; stack <- c(stack, (nj - 1L) * nr + ni)
}
}
}
}
lab
}
## Production implementation: vectorised label propagation (same partition,
## no per-pixel R loop). Each sweep is four shifted pmax operations.
label_cc <- function(mask) {
nr <- nrow(mask); nc <- ncol(mask)
lab <- matrix(0L, nr, nc); lab[mask] <- which(mask)
repeat {
P <- pad_image(lab, 1, 1, "zero")
nb <- pmax(P[1:nr, 2:(nc + 1)], P[3:(nr + 2), 2:(nc + 1)],
P[2:(nr + 1), 1:nc], P[2:(nr + 1), 3:(nc + 2)], lab)
nb[!mask] <- 0L
if (identical(as.vector(nb), as.vector(lab))) break
lab <- nb
}
u <- sort(unique(lab[lab > 0]))
if (length(u)) lab[lab > 0] <- match(lab[lab > 0], u)
matrix(as.integer(lab), nr, nc)
}
fill_holes <- function(mask) {
nr <- nrow(mask); nc <- ncol(mask)
bg <- !mask
lab <- label_cc(bg)
edge <- unique(c(lab[1, ], lab[nr, ], lab[, 1], lab[, nc]))
edge <- edge[edge > 0]
mask | (bg & !matrix(lab %in% edge, nr, nc))
}
## agreement between the two implementations, on a small test mask
demo <- matrix(FALSE, 50, 50)
for (ctr in list(c(12, 12), c(35, 15), c(25, 38)))
demo[(row(demo) - ctr[1])^2 + (col(demo) - ctr[2])^2 <= 36] <- TRUE
c(objects_teaching = max(label_components(demo)),
objects_fast = max(label_cc(demo)),
same_partition = identical(as.vector(table(label_components(demo))),
as.vector(table(label_cc(demo)))))## objects_teaching objects_fast same_partition
## 3 3 1bright <- img_dn > cuts[2] & head_mask
lab_b <- label_cc(bright)
sizes <- tabulate(lab_b[lab_b > 0])
par(mfrow = c(1, 3), mar = c(1, 1, 2.5, 1))
show_img(demo, "Test mask (3 disks)", col = c("black", "white"))
show_img(label_cc(demo), "Labelled components",
col = c("black", "steelblue", "orange", "seagreen"), zlim = c(0, 3))
show_img(lab_b > 0, "Bright class in the phantom", col = c("black", "white"))
par(mfrow = c(1, 1))
data.frame(n_components = max(lab_b), largest_px = max(sizes),
components_over_20px = sum(sizes >= 20))## n_components largest_px components_over_20px
## 1 1 541 1Selecting the target object
“Largest component” is a common but fragile rule. A more defensible rule uses the property that actually defines the target. Below we compare selecting the largest component against selecting the one with the highest mean intensity, which encodes the clinical statement “the lesion is the brightest sizeable structure.” At a permissive threshold the two rules disagree completely: the largest component is a sprawling parenchymal blob, while the brightest is the tumour.
dice_coef <- function(a, b) 2 * sum(a & b) / (sum(a) + sum(b))
select_component <- function(labels, img, min_px = 20,
rule = c("brightest", "largest")) {
rule <- match.arg(rule)
sz <- tabulate(labels[labels > 0])
cand <- which(sz >= min_px)
if (!length(cand)) return(matrix(FALSE, nrow(labels), ncol(labels)))
k <- if (rule == "largest") cand[which.max(sz[cand])]
else cand[which.max(sapply(cand, function(k) mean(img[labels == k])))]
labels == k
}
compare_rules <- function(T) {
lb <- label_cc(img_dn > T & head_mask)
L <- select_component(lb, img_dn, rule = "largest")
B <- select_component(lb, img_dn, rule = "brightest")
data.frame(threshold = T, n_components = max(lb),
largest_px = sum(L), largest_Dice = dice_coef(L, tumor_mask),
brightest_px = sum(B), brightest_Dice = dice_coef(B, tumor_mask))
}
round(rbind(compare_rules(cuts[1]), compare_rules(cuts[2])), 3)## threshold n_components largest_px largest_Dice brightest_px brightest_Dice
## 1 0.412 226 1243 0.000 811 0.811
## 2 0.590 1 541 0.987 541 0.987At the permissive cut the “largest component” rule fails outright; at the strict cut both rules find the tumour. In practice this role is played by a user-placed seed, a bounding box, or an atlas prior, and whichever it is should be reported, because it is part of the segmentation definition.
lab_strict <- label_cc(img_dn > cuts[2] & head_mask)
mask_tumour_raw <- select_component(lab_strict, img_dn, rule = "brightest")
mask_tumour <- fill_holes(bin_open(mask_tumour_raw, struct_elem(1)))
data.frame(raw_px = sum(mask_tumour_raw), cleaned_px = sum(mask_tumour),
true_px = sum(tumor_mask),
Dice = round(dice_coef(mask_tumour, tumor_mask), 3))## raw_px cleaned_px true_px Dice
## 1 541 541 553 0.987Region-based segmentation
Region growing
Region growing starts from one or more seed voxels and iteratively annexes neighboring voxels that are similar enough. Two variants behave very differently: comparing each candidate to the fixed seed value is stable but fails on gradients, while comparing to the running region mean adapts but is prone to leakage across weak boundaries. Both are sensitive to the seed and to the tolerance.
region_grow <- function(img, seed, tol = 0.08, adaptive = FALSE) {
nr <- nrow(img); nc <- ncol(img)
visited <- matrix(FALSE, nr, nc); region <- matrix(FALSE, nr, nc)
ref <- img[seed[1], seed[2]]; s_sum <- ref; s_n <- 1
stack <- matrix(seed, ncol = 2)
while (nrow(stack) > 0) {
p <- stack[nrow(stack), ]; stack <- stack[-nrow(stack), , drop = FALSE]
i <- p[1]; j <- p[2]
if (i < 1 || i > nr || j < 1 || j > nc || visited[i, j]) next
visited[i, j] <- TRUE
target <- if (adaptive) s_sum / s_n else ref
if (abs(img[i, j] - target) <= tol) {
region[i, j] <- TRUE; s_sum <- s_sum + img[i, j]; s_n <- s_n + 1
stack <- rbind(stack, c(i + 1, j), c(i - 1, j), c(i, j + 1), c(i, j - 1))
}
}
region
}
## Seed placed at the tumour centroid. NOTE the index order: the first index is
## the ROW (vertical), the second the COLUMN (horizontal). Swapping them is the
## classic region-growing bug -- the seed silently lands outside the target.
seed <- round(colMeans(which(tumor_mask, arr.ind = TRUE)))
grown_ok <- region_grow(img_dn, seed, tol = 0.18, adaptive = FALSE)
grown_leak <- region_grow(img_dn, seed, tol = 0.30, adaptive = TRUE)
par(mfrow = c(1, 3), mar = c(1, 1, 2.5, 1))
show_img(img_dn, "Input (seed marked)", zlim = c(0, 1))
points(seed[2], n - seed[1] + 1, col = "red", pch = 3, cex = 2, lwd = 2)
show_img(grown_ok, "Well-chosen tolerance (0.18)", col = c("black", "white"))
show_img(grown_leak, "Running mean, tol 0.30: LEAKED", col = c("black", "white"))
Region growing is exquisitely sensitive to its tolerance, and the two variants fail differently. Sweeping the tolerance makes the cliff visible: the running-mean variant tracks a slowly drifting intensity for a while and then, at a critical tolerance, escapes through a weak boundary and floods the entire organ.
tols <- seq(0.06, 0.32, by = 0.02)
sweep_rg <- t(sapply(tols, function(t) c(
fixed = dice_coef(region_grow(img_dn, seed, t, FALSE), tumor_mask),
adaptive = dice_coef(region_grow(img_dn, seed, t, TRUE), tumor_mask))))
matplot(tols, sweep_rg, type = "b", pch = 19, lty = 1, lwd = 2,
col = c("steelblue", "firebrick"), ylim = c(0, 1),
xlab = "intensity tolerance", ylab = "Dice vs. ground truth",
main = "Region growing: a narrow window of good tolerances")
legend("bottomleft", c("fixed seed reference", "running region mean"),
col = c("steelblue", "firebrick"), lwd = 2, pch = 19, bty = "n")
## tolerance fixed adaptive
## 1 0.06 0.900 0.897
## 2 0.08 0.915 0.915
## 3 0.10 0.931 0.936
## 4 0.12 0.954 0.961
## 5 0.14 0.970 0.973
## 6 0.16 0.977 0.990
## 7 0.18 0.991 0.993
## 8 0.20 0.993 0.986
## 9 0.22 0.986 0.971
## 10 0.24 0.973 0.954
## 11 0.26 0.963 0.934
## 12 0.28 0.949 0.068
## 13 0.30 0.933 0.041
## 14 0.32 0.912 0.015A parameter with a cliff is a parameter that must be validated. The running-mean variant is better than the fixed variant over part of the range and catastrophically worse just beyond it. Any pipeline containing such a parameter needs the sensitivity analysis of the next section, not a single default value.
Clustering: k-means and Gaussian mixtures
K-means groups voxels into \(k\) clusters by minimizing within-cluster sum of squares. It is unsupervised and fast, but \(k\) must be chosen, the result depends on initialization, and, importantly, it assigns hard labels and ignores spatial context, so it cannot represent partial-volume voxels. A Gaussian mixture model fitted by EM gives soft memberships \(p(\text{class}\mid I)\), which is the natural representation of a partial-volume voxel, and is the basis of standard tissue-classification tools; adding a Markov random field prior introduces the missing spatial regularity.
vals <- img_dn[head_mask]
set.seed(1)
km <- kmeans(vals, centers = 3, nstart = 20)
ord <- rank(tapply(vals, km$cluster, mean)) # relabel by mean intensity
seg_km <- matrix(0L, n, n); seg_km[head_mask] <- ord[km$cluster]
par(mfrow = c(1, 2), mar = c(1, 1, 2.5, 1))
show_img(seg_km, "K-means (k = 3) inside the head",
col = c("black", "steelblue", "orange", "firebrick"), zlim = c(0, 3))
show_img(seg, "Multi-level Otsu, for comparison",
col = c("black", "steelblue", "orange", "firebrick"), zlim = c(0, 3))
par(mfrow = c(1, 1))
data.frame(cluster = 1:3, center = round(sort(km$centers[, 1]), 3),
n_voxels = as.vector(table(ord[km$cluster])))## cluster center n_voxels
## 3 1 0.388 10924
## 2 2 0.442 3882
## 1 3 0.746 541Edge- and contour-based segmentation
Rather than grouping similar interiors, edge-based methods locate boundaries where the intensity gradient is large, then link edge pixels into closed contours. Active contours (“snakes”) and level sets evolve a curve \(C\) to minimize an energy of the form
\[E(C) = \underbrace{\alpha\!\int\!\lVert C'(s)\rVert^2 ds + \beta\!\int\!\lVert C''(s)\rVert^2 ds}_{\text{internal: smoothness}} \;-\; \underbrace{\gamma\!\int\!\lVert\nabla I(C(s))\rVert^2 ds}_{\text{external: attraction to edges}} ,\]
trading boundary smoothness against edge attraction. Region-based variants (Chan–Vese) replace the gradient term with a piecewise-constant fit, which makes them far more robust on noisy or weakly-edged images, the usual situation in medicine.
Edge-based segmentation excels on high-contrast boundaries but struggles in noisy images, where the gradient is dominated by noise and contours fragment or leak. Edges should therefore be computed after denoising, and preferably with a scale-selective operator such as the Laplacian of Gaussian.
edge_strength <- normalize01(sqrt(conv2(img_dn, Kx)^2 + conv2(img_dn, Ky)^2))
noisy_edges <- normalize01(sqrt(conv2(img_noisy, Kx)^2 + conv2(img_noisy, Ky)^2))
par(mfrow = c(1, 4), mar = c(1, 1, 2.5, 1))
show_img(noisy_edges, "Gradient of the RAW image", zlim = c(0, 1))
show_img(edge_strength, "Gradient after denoising", zlim = c(0, 1))
show_img(edge_strength > 0.25, "Thresholded edges", col = c("black", "white"))
show_img(bin_close(edge_strength > 0.25, struct_elem(2)),
"After closing (linking)", col = c("black", "white"))
Deep learning changes the estimator, not the problem. A U-Net learns the mapping from image patches to labels directly, and typically optimizes a soft Dice loss, \(\mathcal{L}=1-\frac{2\sum_i p_i g_i}{\sum_i p_i+\sum_i g_i}\), or a combination of Dice and cross-entropy. It removes the need to hand-design features, but it does not remove the need for calibration, validation on unseen scanners, uncertainty quantification, or reporting, and it adds a dependence on the training distribution that every other method in this section lacks.
Cell segmentation with EBImage (optional)
Thresholding and connected components cannot separate objects that touch. The watershed transform treats the image (or the distance map of a binary mask) as a topographic surface and floods basins from their minima; the ridges where floods meet become object boundaries. The distance map assigns each foreground pixel its distance to the nearest background pixel, so basins center on object cores and ridges fall at the necks between touching objects. The practical pitfall is over-segmentation, noise creates spurious minima, and the remedy is marker-controlled watershed, flooding only from trusted markers.
Given a set of seeds (nuclei, cell centers, landmarks), a Voronoi
tessellation assigns every pixel to its nearest seed. EBImage’s propagate()
generalizes this to a Voronoi diagram on a Riemann manifold, where the distance
between neighboring points \((x,y,z)\) and \((x+dx,y+dy,z+dz)\) is
\[ds = \sqrt{\tfrac{2}{\lambda+1}\big[\lambda(dx^2+dy^2) + dz^2\big]},\]
with elevation \(z\) given by image intensity and \(\lambda\) trading lateral distance against intensity change. Large \(\lambda\) ignores intensity and gives an ordinary spatial Voronoi diagram; small \(\lambda\) makes boundaries hug intensity edges — ideal for growing cytoplasm masks outward from nuclear seeds.
library(EBImage)
th <- otsu(nuc)
nuc_th <- combine(mapply(function(frame, t) frame > t,
getFrames(nuc), th, SIMPLIFY = FALSE))
EBImage::display(nuc_th, method = "raster", all = TRUE)
library(EBImage)
nmask_bin <- thresh(nuc, w = 10, h = 10, offset = 0.05)
nmask_bin <- opening(nmask_bin, makeBrush(5, shape = "disc"))
nmask_bin <- fillHull(nmask_bin)
wlabels <- watershed(distmap(nmask_bin), tolerance = 1, ext = 1)
EBImage::display(colorLabels(wlabels), method = "raster", all = TRUE)
library(EBImage)
nuc2 <- readImage(system.file("images", "nuclei.tif", package = "EBImage"))
cel <- readImage(system.file("images", "cells.tif", package = "EBImage"))
cells <- rgbImage(green = 1.5 * cel, blue = nuc2)
nmask <- thresh(nuc2, w = 10, h = 10, offset = 0.05)
nmask <- opening(nmask, makeBrush(5, shape = "disc"))
nmask <- fillHull(nmask)
nmask <- bwlabel(nmask) # nuclei = seeds
ctmask <- opening(cel > 0.1, makeBrush(5, shape = "disc"))
cmask <- propagate(cel, seeds = nmask, mask = ctmask) # cytoplasm
segmented <- paintObjects(cmask, cells, col = "#ff00ff")
segmented <- paintObjects(nmask, segmented, col = "#ffff00")
EBImage::display(segmented, method = "raster", all = TRUE)
Applications and characteristic failure modes
| Target | Typical approach | Characteristic difficulty |
|---|---|---|
| Cells and nuclei (microscopy) | threshold + watershed / seeded propagation | many touching objects |
| Organs (CT/MRI) | atlas- and model-based, deep learning | inter-subject variability |
| Tumours and lesions | region/level-set, deep learning, manual | heterogeneous interior, weak margins |
| Vessels and tubular structures | vesselness (Frangi) filters, tracking | thin, branching, partial volume |
| Tissue compartments | GMM + spatial prior | partial volume, bias field |
Watch for these failure modes. Under-segmentation merges objects that should be separate. Over-segmentation splits one object into many. Leakage through weak boundaries pulls the region into adjacent tissue. Inclusion of background or artifact adds non-target voxels. And sensitivity to preprocessing means the same image yields different masks under different smoothing or normalization, so preprocessing must be fixed, documented, and identical across every image in a comparison.
Validating segmentation results
A mask can look plausible yet bias measurements, so segmentation should be validated whenever it feeds quantitative analysis. Validation rests on a reference (ground truth, often expert annotation) and on metrics of two distinct kinds.
Overlap metrics measure how much of the volume agrees:
\[\text{Dice}(A,B)=\frac{2|A\cap B|}{|A|+|B|},\qquad \text{Jaccard}(A,B)=\frac{|A\cap B|}{|A\cup B|},\qquad \text{Jaccard}=\frac{\text{Dice}}{2-\text{Dice}} .\]
Distance metrics measure how far the boundaries are apart, which overlap metrics cannot see:
\[\text{HD}(A,B)=\max\Big\{\sup_{a\in\partial A}\inf_{b\in\partial B}d(a,b),\; \sup_{b\in\partial B}\inf_{a\in\partial A}d(a,b)\Big\},\]
with the 95th-percentile variant HD95 used in practice because the maximum is dominated by a single outlier voxel, and the average symmetric surface distance (ASSD) giving the typical boundary error in millimeters.
jaccard_coef <- function(a, b) sum(a & b) / sum(a | b)
boundary_pixels <- function(mask) {
nr <- nrow(mask); nc <- ncol(mask)
P <- pad_image(mask * 1, 1, 1, "zero")
nb <- P[1:nr, 2:(nc+1)] + P[3:(nr+2), 2:(nc+1)] +
P[2:(nr+1), 1:nc] + P[2:(nr+1), 3:(nc+2)]
which(mask & nb < 4, arr.ind = TRUE)
}
surface_distances <- function(A, B, dy = 1, dx = 1) {
a <- boundary_pixels(A); b <- boundary_pixels(B)
if (!nrow(a) || !nrow(b)) return(c(HD = NA, HD95 = NA, ASSD = NA))
D2 <- outer(a[, 1] * dy, b[, 1] * dy, "-")^2 +
outer(a[, 2] * dx, b[, 2] * dx, "-")^2
dab <- sqrt(apply(D2, 1, min)); dba <- sqrt(apply(D2, 2, min))
c(HD = max(c(dab, dba)),
HD95 = max(quantile(dab, 0.95), quantile(dba, 0.95)),
ASSD = (sum(dab) + sum(dba)) / (length(dab) + length(dba)))
}
confusion <- function(est, gt) {
TP <- sum(est & gt); FP <- sum(est & !gt)
FN <- sum(!est & gt); TN <- sum(!est & !gt)
c(Dice = dice_coef(est, gt), Jaccard = jaccard_coef(est, gt),
Sensitivity = TP / (TP + FN), Specificity = TN / (TN + FP),
Precision = TP / (TP + FP))
}
round(rbind(
`threshold + CC (corrected)` = c(confusion(mask_tumour, tumor_mask),
surface_distances(mask_tumour, tumor_mask, dy, dx)),
`region growing (tol 0.18)` = c(confusion(grown_ok, tumor_mask),
surface_distances(grown_ok, tumor_mask, dy, dx))
), 3)## Dice Jaccard Sensitivity Specificity Precision HD
## threshold + CC (corrected) 0.987 0.975 0.976 1 0.998 1
## region growing (tol 0.18) 0.991 0.982 0.984 1 0.998 1
## HD95 ASSD
## threshold + CC (corrected) 1 0.173
## region growing (tol 0.18) 1 0.126Dice is size-dependent, do not compare it across object sizes. For an object of radius \(R\) whose boundary is misplaced by a uniform \(\varepsilon\), the Dice deficit scales as \(1-\text{Dice}\approx \varepsilon/R\) in 2D (and \(\tfrac{3}{2}\varepsilon/R\) in 3D). A one-pixel boundary error costs almost nothing on a large organ and is catastrophic on a small lesion. The table below makes the point concretely.
dice_for_radius <- function(R, eps = 1) {
S <- 4 * R + 21; cc <- (S + 1) / 2
d <- (row(matrix(0, S, S)) - cc)^2 + (col(matrix(0, S, S)) - cc)^2
A <- d <= R^2
B <- d <= (R + eps)^2 # boundary displaced by eps px
c(radius_px = R, true_px = sum(A), Dice = dice_coef(A, B),
volume_error_pct = 100 * (sum(B) / sum(A) - 1))
}
round(as.data.frame(t(sapply(c(3, 5, 10, 20, 40), dice_for_radius))), 3)## radius_px true_px Dice volume_error_pct
## 1 3 29 0.744 68.966
## 2 5 81 0.835 39.506
## 3 10 317 0.914 18.927
## 4 20 1257 0.956 9.228
## 5 40 5025 0.977 4.697What a good validation reports. A Dice score alone is not a validation. Report (i) an overlap metric, (ii) a boundary metric in millimeters, (iii) the bias in the derived measurement (volume, mean intensity), and (iv) the human inter-observer and intra-observer variability on the same data, which sets the practical ceiling: an algorithm that agrees with experts as well as experts agree with each other is doing as well as the reference permits.
Checkpoint 8. Two methods segment the same 8 mm lesion. Method A: Dice 0.88, HD95 2 mm, volume bias +1 %. Method B: Dice 0.86, HD95 9 mm, volume bias −12 %. Which would you use for a volumetric response study, and why is Dice the wrong tiebreaker? Answer: A, its boundary error and volume bias are both far smaller. Dice differs by 0.02 while the volume bias, which is what a response criterion actually measures, differs by 13 percentage points.
From segmentation to measurement: quantification and uncertainty
A validated mask becomes a measurement only when combined with metadata. This section turns masks into numbers, area, volume, intensity, shape, texture, and, in keeping with this edition’s theme, treats the uncertainty of those numbers as a first-class quantity rather than an afterthought.
The deliverable is a triple. A quantitative imaging result is not a number. It is (value, uncertainty, validation): a calibrated estimate, a defensible interval, and evidence that the pipeline that produced it measures what it claims to measure.
Area, volume, and intensity
For a 2D mask with pixel spacing \(\Delta x,\Delta y\) and \(N\) foreground pixels, the area is \(A = N\,\Delta x\,\Delta y\). For a 3D mask, \(V = N\,\Delta x\,\Delta y\,\Delta z\) (or \(N\,|\det(RS)|\) for oblique geometry). The mean intensity inside the mask, \(\bar I = \frac1N\sum_{(x,y)\in M} I(x,y)\), summarizes the ROI signal.
n_pixels <- sum(mask_tumour)
area_mm2 <- n_pixels * dx * dy
mean_I <- mean(img_dn[mask_tumour])
data.frame(segmented_pixels = n_pixels,
pixel_spacing_mm = dx,
area_mm2 = round(area_mm2, 1),
true_area_mm2 = sum(tumor_mask) * dx * dy,
area_bias_pct = round(100 * (area_mm2 / (sum(tumor_mask) * dx * dy) - 1), 2),
mean_intensity = round(mean_I, 4),
equiv_diameter_mm= round(2 * sqrt(area_mm2 / pi), 2))## segmented_pixels pixel_spacing_mm area_mm2 true_area_mm2 area_bias_pct
## 1 541 1 541 553 -2.17
## mean_intensity equiv_diameter_mm
## 1 0.7461 26.25Clinical example. The same 2,000-pixel lesion is \(125\text{ mm}^2\) at \(0.25\text{ mm}\) spacing and \(500\text{ mm}^2\) at \(0.5\text{ mm}\) spacing. A volume reported without its voxel geometry is not a measurement.
Partial-volume-aware (sub-voxel) volumetry
Binary counting forces every boundary voxel to be entirely in or entirely out, which discretizes the volume in steps of one voxel and biases small objects. A fractional estimate instead uses the intensity itself: if a voxel’s value is a linear mixture of object and background, \(I = f\,I_{\text{obj}} + (1-f)\,I_{\text{bg}}\), then
\[\hat f = \operatorname{clip}\!\left(\frac{I - I_{\text{bg}}} {I_{\text{obj}} - I_{\text{bg}}},\,0,\,1\right), \qquad \hat A = \Delta x\,\Delta y\sum_{\text{voxels}} \hat f .\]
This is the two-compartment model behind sub-voxel tissue quantification, and it usually reduces both the bias and the quantization noise of volumetry.
I_obj <- mean(img_dn[bin_erode(tumor_mask, struct_elem(3))]) # pure-tumour value
I_bg <- mean(img_dn[bin_erode(par_mask, struct_elem(3))]) # pure-parenchyma
roi_band <- bin_dilate(mask_tumour, struct_elem(4)) # object + its rim
frac <- pmin(pmax((img_dn - I_bg) / (I_obj - I_bg), 0), 1)
area_frac<- sum(frac[roi_band]) * dx * dy
truth <- sum(tumor_mask) * dx * dy
data.frame(method = c("binary voxel count", "fractional (partial-volume)"),
area_mm2 = round(c(area_mm2, area_frac), 1),
truth_mm2 = truth,
bias_pct = round(100 * (c(area_mm2, area_frac) / truth - 1), 2))## method area_mm2 truth_mm2 bias_pct
## 1 binary voxel count 541.0 553 -2.17
## 2 fractional (partial-volume) 539.6 553 -2.42Shape features and the digitization bias of the perimeter
Shape descriptors quantify geometry: area/volume, perimeter/surface area, circularity \(4\pi A/P^2\) (1 for a disk, smaller for irregular shapes), elongation, and sphericity in 3D. They help distinguish smooth benign lesions from spiculated malignant ones.
Perimeter, however, is the most frequently mis-measured quantity in image analysis. Counting foreground–background transitions gives the crack boundary length, which for a smooth curve overestimates the true perimeter by a factor \(4/\pi \approx 1.27\), regardless of resolution.
Why \(4/\pi\), and how to fix it. The Cauchy–Crofton formula states that a curve of length \(L\) crossed by a grid of parallel lines with perpendicular spacing \(\lambda\) produces, on average over orientation,
\[L = \frac{\pi}{2}\,\big\langle N_\theta\,\lambda_\theta \big\rangle_\theta ,\]
where \(N_\theta\) counts intersections. Using only the horizontal and vertical pixel grids (\(\lambda = 1\)) gives \(\hat L = \frac{\pi}{4}(N_h+N_v) = \frac{\pi}{4}P_{\text{crack}}\), so the familiar transition count must be multiplied by \(\pi/4\). Adding the two diagonal directions (\(\lambda = 1/\sqrt2\)) gives the more isotropic four-direction estimator
\[\hat L = \frac{\pi}{8}\Big(N_h + N_v + \tfrac{1}{\sqrt2}(N_{\nearrow}+N_{\searrow})\Big).\]
crossings <- function(mask) {
m <- mask * 1; nr <- nrow(m); nc <- ncol(m)
P <- pad_image(m, 1, 1, "zero")
c(Nh = sum(abs(P[2:(nr+1), 1:(nc+1)] - P[2:(nr+1), 2:(nc+2)])), # horizontal
Nv = sum(abs(P[1:(nr+1), 2:(nc+1)] - P[2:(nr+2), 2:(nc+1)])), # vertical
Nd1 = sum(abs(P[1:(nr+1), 1:(nc+1)] - P[2:(nr+2), 2:(nc+2)])), # "\\"
Nd2 = sum(abs(P[2:(nr+2), 1:(nc+1)] - P[1:(nr+1), 2:(nc+2)]))) # "/"
}
perimeter_crack <- function(mask, dy = 1, dx = 1) {
N <- unname(crossings(mask)); N[1] * dy + N[2] * dx }
perimeter_crofton <- function(mask, dy = 1, dx = 1) {
N <- unname(crossings(mask))
(pi / 8) * (N[1] * dy + N[2] * dx + (N[3] + N[4]) * sqrt(dx * dy) / sqrt(2)) }
## convergence on digital disks of known perimeter
conv_tab <- t(sapply(c(5, 10, 20, 40, 80), function(R) {
S <- 2*R + 21; cc <- (S + 1)/2
disk <- (row(matrix(0,S,S)) - cc)^2 + (col(matrix(0,S,S)) - cc)^2 <= R^2
c(R = R, true = 2*pi*R,
crack = perimeter_crack(disk),
crack_x_pi_4 = (pi/4) * perimeter_crack(disk),
crofton4 = perimeter_crofton(disk))
}))
conv_tab <- as.data.frame(conv_tab)
conv_tab$err_crack_pct <- round(100 * (conv_tab$crack / conv_tab$true - 1), 1)
conv_tab$err_crofton_pct <- round(100 * (conv_tab$crofton4 / conv_tab$true - 1), 1)
round(conv_tab, 1)## R true crack crack_x_pi_4 crofton4 err_crack_pct err_crofton_pct
## 1 5 31.4 44 34.6 33.9 40.1 8.0
## 2 10 62.8 84 66.0 65.2 33.7 3.8
## 3 20 125.7 164 128.8 127.7 30.5 1.6
## 4 40 251.3 324 254.5 252.7 28.9 0.6
## 5 80 502.7 644 505.8 505.0 28.1 0.5The naive crack perimeter is 28–40 % too long and does not improve with resolution; the Crofton estimator converges. The consequence for circularity is severe:
circularity <- function(A, P) 4 * pi * A / P^2
A_px <- sum(mask_tumour)
data.frame(
estimator = c("crack transitions", "Crofton (4 directions)"),
perimeter_mm= round(c(perimeter_crack(mask_tumour, dy, dx),
perimeter_crofton(mask_tumour, dy, dx)), 1),
circularity = round(c(circularity(A_px, perimeter_crack(mask_tumour)),
circularity(A_px, perimeter_crofton(mask_tumour))), 3),
note = c("biased low by ~(pi/4)^2 = 0.62", "unbiased for smooth shapes"))## estimator perimeter_mm circularity
## 1 crack transitions 108.0 0.583
## 2 Crofton (4 directions) 84.6 0.949
## note
## 1 biased low by ~(pi/4)^2 = 0.62
## 2 unbiased for smooth shapesCheck what your software computes. A “circularity” of 0.62 for a perfect circle is the signature of an uncorrected crack perimeter. Shape features are routinely compared across studies and software packages that use different perimeter conventions, a difference of 27 % in \(P\) becomes 61 % in \(4\pi A/P^2\).
Texture features
Texture captures spatial intensity patterns that the mean misses (homogeneous vs. heterogeneous tumours). A classic family derives from the gray-level co-occurrence matrix (GLCM) \(P(i,j)\): the normalized frequency with which intensity level \(i\) occurs a fixed offset away from level \(j\). Standard summaries include
\[\text{contrast}=\sum_{i,j}P(i,j)(i-j)^2,\quad \text{homogeneity}=\sum_{i,j}\frac{P(i,j)}{1+(i-j)^2},\] \[\text{energy}=\sum_{i,j}P(i,j)^2,\quad \text{entropy}=-\sum_{i,j}P(i,j)\log_2 P(i,j).\]
Three implementation choices change the numbers substantially and must always be reported: the number of gray levels, whether the matrix is symmetrized (\(P + P^{\top}\)), and whether it is averaged over multiple directions (which makes it rotation-robust).
glcm_features <- function(img, mask = NULL, levels = 32,
offsets = list(c(0,1), c(-1,1), c(-1,0), c(-1,-1))) {
v <- if (is.null(mask)) img else { z <- img; z[!mask] <- NA; z }
rng <- range(v, na.rm = TRUE)
q <- floor((v - rng[1]) / (diff(rng) + 1e-12) * (levels - 1e-9)) + 1
q[q > levels] <- levels
nr <- nrow(q); nc <- ncol(q); P <- matrix(0, levels, levels)
for (o in offsets) {
di <- o[1]; dj <- o[2]
i1 <- max(1, 1 - di):min(nr, nr - di)
j1 <- max(1, 1 - dj):min(nc, nc - dj)
a <- q[i1, j1]; b <- q[i1 + di, j1 + dj]
ok <- !is.na(a) & !is.na(b)
if (any(ok)) {
tt <- table(factor(a[ok], 1:levels), factor(b[ok], 1:levels))
P <- P + tt + t(tt) # symmetrize
}
}
P <- P / sum(P); ii <- row(P); jj <- col(P)
mu <- sum(P * ii); sdv <- sqrt(sum(P * (ii - mu)^2))
c(contrast = sum(P * (ii - jj)^2),
homogeneity = sum(P / (1 + (ii - jj)^2)),
energy = sum(P^2),
entropy_bits= -sum(P[P > 0] * log2(P[P > 0])),
correlation = sum(P * (ii - mu) * (jj - mu)) / sdv^2)
}
## sensitivity of texture features to two arbitrary implementation choices
rbind(
`denoised, 32 levels` = glcm_features(img_dn, mask_tumour, levels = 32),
`denoised, 8 levels` = glcm_features(img_dn, mask_tumour, levels = 8),
`raw, 32 levels` = glcm_features(img_noisy, mask_tumour, levels = 32),
`smoothed, 32 levels` = glcm_features(gauss_blur(img_dn, 2), mask_tumour, levels = 32)
)## contrast homogeneity energy entropy_bits correlation
## denoised, 32 levels 23.379 0.3938 0.015286 7.419 0.6728
## denoised, 8 levels 1.475 0.6972 0.134522 3.998 0.6780
## raw, 32 levels 30.892 0.1993 0.005044 8.194 0.3498
## smoothed, 32 levels 17.727 0.4525 0.030642 7.053 0.8747Texture features are not intrinsic properties of the tissue. Changing the number of gray levels or adding a mild smoothing step alters contrast, entropy, and correlation by tens of percent, far more than most reported biological effects. This is why the Image Biomarker Standardisation Initiative (IBSI) exists, and why an unstandardized radiomic signature rarely replicates.
Uncertainty in image-derived measurements
Every number above is an estimate. Reporting it without an uncertainty overstates what the image can tell us. We follow the standard metrological framework (the Guide to the Expression of Uncertainty in Measurement, GUM), which classifies components by how they are evaluated, not by their cause:
- Type A, evaluated by statistical analysis of repeated observations (noise-driven variability, test–retest scatter);
- Type B, evaluated by other means: calibration certificates, manufacturer specifications, physical reasoning, or expert judgement (voxel-size accuracy, scanner calibration, partial-volume bias).
For a measurand \(y = f(x_1,\dots,x_m)\) the combined standard uncertainty is
\[u_c^2(y) = \sum_{i=1}^{m}\left(\frac{\partial f}{\partial x_i}\right)^{\!2}u^2(x_i) \;+\;2\sum_{i<j}\frac{\partial f}{\partial x_i}\frac{\partial f}{\partial x_j}\, u(x_i,x_j),\]
and an expanded uncertainty \(U = k\,u_c\) with coverage factor \(k = 2\) gives an interval of roughly 95 % confidence. The cross-terms matter: pipeline steps are generally not independent.
Worked example: uncertainty of a volume
Volume is \(V = N v\), a product, so relative uncertainties add in quadrature:
\[\left(\frac{u_V}{V}\right)^2 = \left(\frac{u_N}{N}\right)^2 + \left(\frac{u_v}{v}\right)^2 .\]
The voxel-volume term is Type B (from the geometric calibration of the scanner); the count term is Type A and is dominated by the boundary. For an object of volume \(V\) and surface area \(S\) whose boundary position is uncertain by \(\varepsilon\),
\[u_V \approx \varepsilon\,S \quad\Longrightarrow\quad \frac{u_V}{V} \approx \varepsilon\,\frac{S}{V} = \frac{3\varepsilon}{R} \ \ \text{for a sphere of radius } R .\]
This single formula explains why small and convoluted objects are intrinsically harder to measure than large compact ones, an effect of geometry, not of algorithm quality.
u_rel_volume <- function(R_mm, eps_mm, u_calib_rel = 0.01) {
u_N_rel <- 3 * eps_mm / R_mm # boundary term (sphere)
sqrt(u_N_rel^2 + u_calib_rel^2)
}
tab <- expand.grid(R_mm = c(2, 5, 10, 20, 40), eps_mm = c(0.5, 1.0))
tab$u_rel_pct <- round(100 * u_rel_volume(tab$R_mm, tab$eps_mm), 1)
tab$V_mm3 <- round(4/3 * pi * tab$R_mm^3, 0)
tab$U95_mm3 <- round(2 * tab$V_mm3 * u_rel_volume(tab$R_mm, tab$eps_mm), 0)
tab## R_mm eps_mm u_rel_pct V_mm3 U95_mm3
## 1 2 0.5 75.0 34 51
## 2 5 0.5 30.0 524 315
## 3 10 0.5 15.0 4189 1259
## 4 20 0.5 7.6 33510 5071
## 5 40 0.5 3.9 268083 20809
## 6 2 1.0 150.0 34 102
## 7 5 1.0 60.0 524 629
## 8 10 1.0 30.0 4189 2515
## 9 20 1.0 15.0 33510 10075
## 10 40 1.0 7.6 268083 40568The trap that catches everyone: correlated noise and the effective sample size
The standard error of a mean over \(N\) voxels is \(\sigma/\sqrt{N}\) only if the voxels are independent. Reconstruction, interpolation, and any smoothing induce spatial correlation, so the true standard error is
\[\text{SE}(\bar I) = \frac{\sigma}{\sqrt{N_{\text{eff}}}}, \qquad N_{\text{eff}} = \frac{N}{1 + 2\sum_{k\ge1}\rho_k} \;\approx\; \frac{N}{4\pi\sigma_g^2}\]
for a field smoothed with a 2D Gaussian of width \(\sigma_g\) pixels. The naive formula therefore understates the standard error by a factor \(2\sigma_g\sqrt{\pi}\), a factor of 7 for \(\sigma_g = 2\) px. Let us verify this by brute force.
N2 <- 128
roi <- matrix(FALSE, N2, N2); roi[35:94, 35:94] <- TRUE # 60 x 60 ROI
set.seed(77)
check_neff <- function(sigma_g, reps = 120) {
means <- replicate(reps, {
z <- matrix(rnorm(N2 * N2, 0, 1), N2, N2)
z <- if (sigma_g > 0) gauss_blur(z, sigma_g) else z
mean(z[roi])
})
one <- matrix(rnorm(N2 * N2, 0, 1), N2, N2)
one <- if (sigma_g > 0) gauss_blur(one, sigma_g) else one
naive <- sd(one[roi]) / sqrt(sum(roi))
c(sigma_g = sigma_g, naive_SE = naive, empirical_SE = sd(means),
ratio = sd(means) / naive, theory = if (sigma_g > 0) 2 * sigma_g * sqrt(pi) else 1)
}
neff <- as.data.frame(t(sapply(c(0, 1, 2, 3), check_neff)))
round(neff, 4)## sigma_g naive_SE empirical_SE ratio theory
## 1 0 0.0164 0.0164 0.997 1.000
## 2 1 0.0047 0.0173 3.683 3.545
## 3 2 0.0023 0.0164 6.989 7.090
## 4 3 0.0016 0.0168 10.501 10.635plot(neff$sigma_g, neff$ratio, type = "b", pch = 19, col = "firebrick",
xlab = expression(smoothing~sigma[g]~(pixels)),
ylab = "true SE / naive SE",
main = "Smoothing inflates the standard error of an ROI mean")
lines(neff$sigma_g, neff$theory, col = "steelblue", lty = 2, lwd = 2)
legend("topleft", c("empirical", expression(2*sigma[g]*sqrt(pi))),
col = c("firebrick", "steelblue"), lty = c(1, 2), pch = c(19, NA), bty = "n")
This is the single most common statistical error in quantitative imaging. Computing “mean \(\pm\) SD\(/\sqrt{N}\)” over a smoothed ROI can understate the standard error by an order of magnitude, turning noise into a significant finding. Either estimate \(N_{\text{eff}}\) from the autocorrelation, or, far safer, obtain the uncertainty from repeated measurements rather than from within-image voxel statistics.
Sensitivity analysis: how much does the answer depend on the choices?
The cheapest useful uncertainty estimate sweeps each free parameter across a plausible band and reports the spread of the resulting measurement. Here we sweep the segmentation threshold, re-running the complete object-selection step each time so that the measurement is the one we would actually report.
measure_at_threshold <- function(T) {
lb <- label_cc(img_dn > T & head_mask)
m <- select_component(lb, img_dn, rule = "brightest")
if (!sum(m)) return(c(area_mm2 = NA, Dice = NA))
m <- fill_holes(bin_open(m, struct_elem(1)))
c(area_mm2 = sum(m) * dx * dy, Dice = dice_coef(m, tumor_mask))
}
band <- seq(cuts[2] - 0.06, cuts[2] + 0.06, by = 0.01)
sweep <- as.data.frame(cbind(threshold = band, t(sapply(band, measure_at_threshold))))
round(sweep, 3)## threshold area_mm2 Dice
## 1 0.53 579 0.977
## 2 0.54 571 0.984
## 3 0.55 568 0.987
## 4 0.56 560 0.990
## 5 0.57 556 0.992
## 6 0.58 549 0.989
## 7 0.59 541 0.987
## 8 0.60 531 0.980
## 9 0.61 525 0.974
## 10 0.62 523 0.972
## 11 0.63 516 0.965
## 12 0.64 512 0.962
## 13 0.65 502 0.952plot(sweep$threshold, sweep$area_mm2, type = "b", pch = 19, col = "purple",
xlab = "segmentation threshold", ylab = expression(measured~area~(mm^2)),
main = "Measurement sensitivity to the threshold")
abline(v = cuts[2], col = "red", lty = 2)
abline(h = sum(tumor_mask) * dx * dy, col = "darkgreen", lty = 3)
legend("topright", c("Otsu cut", "true area"), col = c("red", "darkgreen"),
lty = c(2, 3), bty = "n")
u_thresh <- (max(sweep$area_mm2) - min(sweep$area_mm2)) / 2
c(area_at_Otsu_mm2 = round(measure_at_threshold(cuts[2])["area_mm2"], 1),
half_range_mm2 = round(u_thresh, 1),
half_range_pct = round(100 * u_thresh / mean(sweep$area_mm2), 1))## area_at_Otsu_mm2.area_mm2 half_range_mm2 half_range_pct
## 541.0 38.5 7.1Monte-Carlo propagation through the whole pipeline
A sensitivity sweep varies one knob. To capture the combined effect of noise propagating through bias correction, denoising, thresholding, component selection and mask cleaning, we re-run the entire pipeline over independent noise realizations of the same underlying object. This is Type A uncertainty evaluated honestly, and it also exposes the bias that no within-image statistic can see.
make_phantom <- function(tumor_a = 0.20, tumor_b = 0.14, seed = 1,
sigma = sigma_noise, psf = psf_sigma) {
brain <- ellipse(X, Y, 0.00, 0.00, 0.82, 0.95)
tumour <- ellipse(X, Y, 0.34, -0.34, tumor_a, tumor_b, angle = 0.5)
vessel <- abs(Y - 0.35 * sin(5 * X)) < 0.022 & abs(X) < 0.70 & brain
lesion <- ellipse(X, Y, -0.42, 0.34, 0.06, 0.06)
f <- matrix(0.05, n, n)
f[brain] <- 0.40; f[vessel] <- 0.62; f[lesion] <- 0.55; f[tumour] <- 0.80
set.seed(seed)
list(g = bias_field * gauss_blur(f, psf) + matrix(rnorm(n*n, 0, sigma), n, n),
truth = tumour)
}
## the complete measurement pipeline, exactly as assembled in this chapter
measure_lesion <- function(g) {
hm <- fill_holes(bin_open(g > otsu_threshold(g), struct_elem(2)))
d <- median_filter(g / estimate_bias(g, hm, degree = 2), 3)
ct <- multi_otsu(d[hm])
lb <- label_cc(d > ct[2] & hm)
m <- select_component(lb, d, rule = "brightest")
if (!sum(m)) return(list(area = NA, mask = NULL, mean_I = NA))
m <- fill_holes(bin_open(m, struct_elem(1)))
list(area = sum(m) * dx * dy, mask = m, mean_I = mean(d[m]))
}n_mc <- 30
mc <- t(sapply(1:n_mc, function(s) {
p <- make_phantom(seed = 1000 + s)
r <- measure_lesion(p$g)
c(area = r$area, mean_I = r$mean_I, Dice = dice_coef(r$mask, p$truth))
}))
truth_area <- sum(make_phantom(seed = 1)$truth) * dx * dy
par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
hist(mc[, "area"], breaks = 12, col = "gray80",
main = "Monte-Carlo distribution of the measured area",
xlab = expression(area~(mm^2)))
abline(v = truth_area, col = "darkgreen", lwd = 2)
abline(v = mean(mc[, "area"]), col = "firebrick", lwd = 2, lty = 2)
legend("topright", c("truth", "mean estimate"), col = c("darkgreen", "firebrick"),
lwd = 2, lty = c(1, 2), bty = "n", cex = 0.8)
plot(mc[, "area"], mc[, "Dice"], pch = 19, col = "steelblue",
xlab = expression(measured~area~(mm^2)), ylab = "Dice",
main = "Overlap vs. measured area")
par(mfrow = c(1, 1))
data.frame(
replicates = n_mc,
truth_mm2 = truth_area,
mean_mm2 = round(mean(mc[, "area"]), 1),
bias_mm2 = round(mean(mc[, "area"]) - truth_area, 1),
bias_pct = round(100 * (mean(mc[, "area"]) / truth_area - 1), 2),
SD_mm2 = round(sd(mc[, "area"]), 2),
CV_pct = round(100 * sd(mc[, "area"]) / mean(mc[, "area"]), 2),
mean_Dice = round(mean(mc[, "Dice"]), 3))## replicates truth_mm2 mean_mm2 bias_mm2 bias_pct SD_mm2 CV_pct mean_Dice
## 1 30 553 540.2 -12.8 -2.31 7.02 1.3 0.985Bias and precision are different quantities and need different evidence. The Monte-Carlo spread (a fraction of a percent here) measures precision only. The systematic offset from the truth, driven by the PSF and by where the threshold sits on the blurred edge, is bias, and no amount of repetition reveals it. Bias requires a phantom with known truth; precision requires repeats. A quantitative imaging biomarker needs both.
Repeatability: Bland–Altman, RC, wCV, and ICC
Precision claims are conventionally reported using test–retest data across several subjects. We simulate ten “patients” with different tumour sizes, each imaged twice with independent noise.
sizes <- seq(0.12, 0.26, length.out = 10)
tr <- t(sapply(seq_along(sizes), function(k) {
p1 <- make_phantom(sizes[k], sizes[k] * 0.7, seed = 100 + k)
p2 <- make_phantom(sizes[k], sizes[k] * 0.7, seed = 200 + k)
c(truth = sum(p1$truth) * dx * dy,
test = measure_lesion(p1$g)$area,
retest= measure_lesion(p2$g)$area)
}))
tr <- as.data.frame(tr)
d <- tr$test - tr$retest
mn <- (tr$test + tr$retest) / 2
bias <- mean(d); s_d <- sd(d)
LoA <- bias + c(-1.96, 1.96) * s_d
RC <- 1.96 * s_d # repeatability coefficient
wCV <- (s_d / sqrt(2)) / mean(mn) # within-subject coefficient of variation
icc_agreement <- function(a, b) {
M <- cbind(a, b); k <- 2; N <- nrow(M); gm <- mean(M)
MSB <- k * sum((rowMeans(M) - gm)^2) / (N - 1)
MSW <- sum((M - rowMeans(M))^2) / (N * (k - 1))
(MSB - MSW) / (MSB + (k - 1) * MSW)
}
par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
plot(mn, d, pch = 19, col = "steelblue", ylim = range(c(d, LoA)) * 1.2,
xlab = expression(mean~of~test~and~retest~(mm^2)),
ylab = expression(test - retest~(mm^2)),
main = "Bland-Altman")
abline(h = bias, col = "firebrick", lwd = 2)
abline(h = LoA, col = "firebrick", lty = 2)
plot(tr$truth, mn, pch = 19, col = "darkgreen",
xlab = expression(true~area~(mm^2)), ylab = expression(measured~area~(mm^2)),
main = "Accuracy: measured vs. true")
abline(0, 1, lty = 2)
par(mfrow = c(1, 1))
data.frame(
bias_mm2 = round(bias, 2),
SD_of_differences = round(s_d, 2),
LoA_lower = round(LoA[1], 1),
LoA_upper = round(LoA[2], 1),
repeatability_coef = round(RC, 1),
within_subject_CV_pct = round(100 * wCV, 2),
ICC_agreement = round(icc_agreement(tr$test, tr$retest), 4),
mean_accuracy_bias_pct = round(100 * mean(mn / tr$truth - 1), 2))## bias_mm2 SD_of_differences LoA_lower LoA_upper repeatability_coef
## 1 -2.9 12.29 -27 21.2 24.1
## within_subject_CV_pct ICC_agreement mean_accuracy_bias_pct
## 1 1.69 0.9988 -3.22What each statistic is for.
- Bland–Altman bias and limits of agreement, how much two measurements of the same subject may differ. The LoA, not the mean difference, is what a clinician needs.
- Repeatability coefficient \(\text{RC}=1.96\sqrt{2}\,\sigma_w = 2.77\,\sigma_w\) , the smallest change in an individual that exceeds measurement noise. A reported “20 % progression” is meaningless if it is inside the RC.
- Within-subject CV, the same quantity expressed relatively, so it can be compared across lesion sizes and sites.
- ICC, the fraction of total variance that is between subjects. It is high here because the simulated tumours differ greatly in size; a high ICC in a homogeneous population is much harder to achieve, which is why ICC must always be reported with the range of the measurand.
Precision without accuracy is the classic failure. Our pipeline is extremely repeatable (within-subject CV around 1 %) and systematically biased by a few percent because thresholding a blurred edge under-counts the object. A biomarker can be exquisitely reproducible and still wrong, which is fine for tracking change within a patient and fatal for comparing against an absolute cut-off.
Error amplification: why diameters are worse than volumes
Response criteria such as RECIST use the longest diameter, but the biology scales with volume, \(V \propto D^3\). Differentiating,
\[\frac{u_V}{V} = 3\,\frac{u_D}{D},\]
so a 10 % diameter error becomes a 30 % volume error, and, conversely, the familiar RECIST thresholds correspond to much larger volume changes than they appear to.
data.frame(
criterion = c("RECIST partial response (-30% diameter)",
"RECIST progressive disease (+20% diameter)",
"measurement error of +/-1 mm on a 10 mm lesion"),
diameter_change_pct = c(-30, 20, 10),
volume_change_pct = round(100 * ((1 + c(-0.30, 0.20, 0.10))^3 - 1), 1))## criterion diameter_change_pct
## 1 RECIST partial response (-30% diameter) -30
## 2 RECIST progressive disease (+20% diameter) 20
## 3 measurement error of +/-1 mm on a 10 mm lesion 10
## volume_change_pct
## 1 -65.7
## 2 72.8
## 3 33.1Take-home. A measurement without an uncertainty is incomplete. A reproducible pipeline plus a sensitivity analysis and a Monte-Carlo or test–retest study turns “the lesion is 545 mm²” into “the lesion is 545 ± 9 mm² (k = 2), with a partial-volume bias of about −2 %, Dice 0.99 against a known phantom”, which is what downstream modelling and clinical decisions actually require.
Reporting checklist for a quantitative imaging measurement
State all of it, or the number is not reproducible.
- Acquisition: modality, scanner, sequence/protocol, reconstruction kernel, voxel size, slice thickness.
- Calibration: intensity units and how they were obtained (rescale slope/intercept, SUV formula, phantom calibration).
- Preprocessing: resampling (interpolator, target grid, how many times), bias correction (method, degree), denoising (filter and parameters), normalization.
- Segmentation: method, all parameters, seed/prior source, mask cleaning operations, and who or what drew it.
- Measurement: the formula, the perimeter/surface convention, and any partial-volume correction.
- Uncertainty: Type A (repeats or Monte Carlo) and Type B (calibration) components, the coverage factor, and the repeatability coefficient.
- Validation: overlap and boundary metrics against a stated reference, plus inter-observer variability on the same data.
Checkpoint 9. A paper reports mean ADC in a tumour ROI as \(1.12 \pm 0.02\times10^{-3}\) mm²/s, with the uncertainty computed as SD\(/\sqrt{N}\) over 900 voxels of a Gaussian-smoothed map (\(\sigma_g = 2\) px). By roughly what factor is the reported uncertainty too small, and what should they have done? Answer: about \(2\sigma_g\sqrt\pi \approx 7\), so the honest figure is nearer \(\pm 0.14\). They should have estimated \(N_{\text{eff}}\), or better, reported a test–retest repeatability coefficient.
Feature extraction, radiomics, and modeling
The preceding sections converted a raw acquisition into a validated mask and a set of numbers with uncertainties. This section closes the loop from images to inference: how individual measurements become feature vectors, how those vectors are condensed through dimensionality reduction, and how they feed predictive models, together with the specific ways this pipeline produces results that do not replicate.
From a region to a feature vector
A radiomic analysis represents each segmented region not by a single number but by a high-dimensional vector:
- First-order (histogram) features, mean, variance, skewness, kurtosis, energy, entropy, and percentiles of the intensities inside the mask. These ignore spatial arrangement entirely.
- Shape features, area/volume, perimeter/surface area, compactness, sphericity, elongation. These ignore intensity entirely.
- Texture features, GLCM contrast, homogeneity, energy, entropy, correlation, plus run-length and size-zone descriptors. These capture spatial pattern.
- Higher-order / filtered features, the same statistics recomputed after wavelet or Laplacian-of-Gaussian filtering, which emphasize structure at chosen spatial scales.
The function below assembles a compact feature vector for our phantom lesion by reusing helpers built earlier in the chapter. Note that it uses the Crofton perimeter and the symmetrized, multi-direction GLCM, a real study should likewise use a curated, standardized library (for example one compliant with the Image Biomarker Standardisation Initiative).
first_order <- function(img, mask) {
v <- img[mask]; m <- mean(v); s <- sd(v)
c(mean = m, sd = s,
skewness = mean((v - m)^3) / s^3,
kurtosis = mean((v - m)^4) / s^4 - 3,
p10 = unname(quantile(v, 0.10)), p90 = unname(quantile(v, 0.90)),
IQR = unname(diff(quantile(v, c(0.25, 0.75)))))
}
shape_feats <- function(mask, dy = 1, dx = 1) {
A <- sum(mask) * dy * dx
P <- perimeter_crofton(mask, dy, dx)
ij <- which(mask, arr.ind = TRUE)
ij_mm <- cbind(ij[, 1] * dy, ij[, 2] * dx) # rows -> y (mm), cols -> x (mm)
ev <- sort(eigen(cov(ij_mm), symmetric = TRUE)$values, decreasing = TRUE)
c(area = A, perimeter = P, circularity = 4 * pi * A / P^2,
equiv_diameter = 2 * sqrt(A / pi),
elongation = sqrt(ev[2] / ev[1]))
}
feature_vector <- c(first_order(img_dn, mask_tumour),
shape_feats(mask_tumour, dy, dx),
glcm_features(img_dn, mask_tumour, levels = 32))
round(feature_vector, 4)## mean sd skewness kurtosis p10
## 0.7461 0.0485 -1.4655 1.4508 0.6585
## p90 IQR area perimeter circularity
## 0.7868 0.0416 541.0000 84.6189 0.9495
## equiv_diameter elongation contrast homogeneity energy
## 26.2454 0.7032 23.3791 0.3938 0.0153
## entropy_bits correlation
## 7.4190 0.6728Radiomics workflow. (1) acquire and reconstruct; (2) harmonize, resample to a common voxel size and normalize intensities; (3) segment the region of interest; (4) extract standardized features; (5) select / reduce features; (6) model; (7) validate on independent data. Steps 1–3 are the image-processing pipeline of this entire chapter; steps 4–7 are statistical learning. Errors in the early steps propagate into every downstream feature.
Feature stability: which features survive a repeat scan?
Before asking whether a feature is predictive, ask whether it is reproducible. A feature whose test–retest ICC is 0.4 cannot support a prediction, however strong its apparent association. We already have the machinery: re-run the pipeline over independent noise realizations and compute, for each feature, the fraction of variance that is between-subject rather than within-subject.
## 8 "subjects" of differing lesion size, each imaged twice
sizes8 <- seq(0.13, 0.25, length.out = 8)
extract_all <- function(g) {
r <- measure_lesion(g)
if (is.null(r$mask)) return(rep(NA, 17))
d <- median_filter(g / estimate_bias(g, fill_holes(bin_open(
g > otsu_threshold(g), struct_elem(2))), 2), 3)
c(first_order(d, r$mask), shape_feats(r$mask, dy, dx),
glcm_features(d, r$mask, levels = 32))
}
F1 <- t(sapply(seq_along(sizes8), function(k)
extract_all(make_phantom(sizes8[k], sizes8[k]*0.7, seed = 300 + k)$g)))
F2 <- t(sapply(seq_along(sizes8), function(k)
extract_all(make_phantom(sizes8[k], sizes8[k]*0.7, seed = 400 + k)$g)))
icc_features <- sapply(seq_len(ncol(F1)), function(j)
icc_agreement(F1[, j], F2[, j]))
names(icc_features) <- names(feature_vector)
op <- par(mar = c(4, 9, 3, 1))
barplot(sort(icc_features), horiz = TRUE, las = 1, xlim = c(0, 1),
col = ifelse(sort(icc_features) > 0.9, "seagreen",
ifelse(sort(icc_features) > 0.75, "goldenrod", "firebrick")),
xlab = "test-retest ICC", main = "Feature stability under repeat imaging",
cex.names = 0.7)
abline(v = c(0.75, 0.9), lty = 2)
## equiv_diameter area perimeter kurtosis p90
## 1.000 0.999 0.999 0.917 0.917
## mean IQR skewness homogeneity contrast
## 0.913 0.897 0.897 0.811 0.780
## sd correlation circularity p10 elongation
## 0.610 0.601 0.522 0.443 0.313
## entropy_bits energy
## 0.188 0.017Shape and first-order features are highly repeatable here; several texture features are not. A defensible radiomics study pre-filters on stability (typically ICC > 0.75 or 0.9 on test–retest or phantom data) before any association testing, otherwise unstable features consume degrees of freedom and generate false positives.
Batch effects and reproducibility. Radiomic features are notoriously sensitive to scanner, reconstruction kernel, slice thickness, and segmentation. Features that look predictive can simply be encoding the site rather than the biology. Fix and document voxel size, normalization, and segmentation; use phantom or test–retest data to discard unstable features; and prefer explicit harmonization (e.g. ComBat-style adjustment) before modeling.
Dimensionality reduction
A feature matrix typically has many more columns than independent directions. Dimensionality reduction compresses these into a few informative coordinates for visualization, denoising, and as input to a classifier.
| Method | Idea | Preserves | Typical use |
|---|---|---|---|
| PCA | orthogonal directions of maximum variance | global linear structure | decorrelation, denoising, first look |
| ICA | statistically independent (non-Gaussian) sources | independent signals | source separation (e.g. fMRI) |
| t-SNE | match local neighbor probabilities in low-D | local clusters | visualizing cluster structure |
| UMAP | manifold graph + low-D layout | local + some global | fast embedding of large data |
set.seed(7)
## NOTE: a local sample size is named n_obs, NOT n -- overwriting the image
## size n is a classic and silent source of downstream breakage.
n_obs <- 60
latent <- rep(1:3, each = n_obs / 3) # three latent "tissue classes"
F_a <- latent + rnorm(n_obs, 0, 0.30)
F_b <- 2 * latent + rnorm(n_obs, 0, 0.40) # correlated with F_a
F_c <- rnorm(n_obs, 0, 1.00) # pure noise feature
Xf <- scale(cbind(F_a, F_b, F_c))
pc <- prcomp(Xf)
par(mfrow = c(1, 2), mar = c(4, 4, 3, 1))
barplot(summary(pc)$importance[2, ], col = "steelblue", ylim = c(0, 1),
ylab = "proportion of variance", main = "PCA scree")
plot(pc$x[, 1], pc$x[, 2], col = latent, pch = 19,
xlab = "PC1", ylab = "PC2", main = "Correlated features collapse onto PC1")
par(mfrow = c(1, 1))
c(image_size_n_unchanged = n, PC1_variance = round(summary(pc)$importance[2, 1], 3))## image_size_n_unchanged PC1_variance
## 160.000 0.657Nonlinear embedding (optional, network-dependent)
When classes lie on a curved manifold, linear projections blur them together.
t-SNE instead preserves local neighborhoods. The canonical demonstration
embeds thousands of \(28\times28\) handwritten digit images into the plane; it is
adapted from
DSPA2 Chapter 4
and gated behind run_heavy because it downloads a 42,000-image dataset.
library(Rtsne); library(httr)
httr::set_config(config(ssl_verifypeer = 0L))
zipURL <- paste0("https://socr.umich.edu/people/dinov/2017/Spring/",
"DSPA_HS650/data/DigitRecognizer_TrainingData.zip")
zipFile <- file.path(tempdir(), "digits.zip")
invisible(GET(zipURL, write_disk(zipFile, overwrite = TRUE)))
train <- read.csv(unzip(zipFile, exdir = tempdir()))
labels <- factor(train$label)
Xd <- as.matrix(train[, -1]) / 255
set.seed(1); idx <- sample(nrow(Xd), 5000)
emb <- Rtsne(Xd[idx, ], dims = 2, perplexity = 30, max_iter = 500,
check_duplicates = FALSE)
plot(emb$Y, col = rainbow(10)[as.integer(labels[idx])], pch = 19, cex = 0.5,
xlab = "t-SNE 1", ylab = "t-SNE 2",
main = "t-SNE embedding of handwritten digits")
legend("topright", legend = 0:9, col = rainbow(10), pch = 19, cex = 0.6, ncol = 2)Embeddings are for insight, not measurement. t-SNE and UMAP distances and cluster sizes are not faithful to the original feature space; they are tuned to reveal local neighborhoods. Use them to explore and to sanity-check class separability, but never read quantitative distances off the plot, and always re-run with several seeds and perplexities before trusting a structure.
From features to prediction, and how it goes wrong
The feature vector is the input to a predictive model, logistic regression, penalized regression, random forests, support-vector machines, or a neural network, that maps imaging features to an outcome. The modeling machinery is developed in later chapters; here we demonstrate the failure mode that dominates image-derived modeling.
Selection leakage, demonstrated
With hundreds of features and tens of patients, some feature will correlate with the outcome by chance. If features are chosen using the whole dataset and only then cross-validated, the test folds have already influenced the model. The result is a confident, publishable, and entirely spurious accuracy.
set.seed(11)
n_pat <- 40; n_feat <- 500
Xr <- matrix(rnorm(n_pat * n_feat), n_pat, n_feat) # PURE NOISE
yr <- rep(0:1, each = n_pat / 2) # label independent of X
cv_accuracy <- function(select_inside_fold) {
folds <- sample(rep(1:5, length.out = n_pat))
if (!select_inside_fold) { # LEAKY: select once, on all data
r_all <- abs(apply(Xr, 2, function(z) cor(z, yr)))
keep_all <- order(r_all, decreasing = TRUE)[1:5]
}
acc <- sapply(1:5, function(k) {
tr <- folds != k; te <- !tr
keep <- if (select_inside_fold) {
r <- abs(apply(Xr[tr, , drop = FALSE], 2, function(z) cor(z, yr[tr])))
order(r, decreasing = TRUE)[1:5]
} else keep_all
fit <- glm(yr[tr] ~ ., data = data.frame(Xr[tr, keep, drop = FALSE]),
family = binomial)
pred <- predict(fit, newdata = data.frame(Xr[te, keep, drop = FALSE]),
type = "response") > 0.5
mean(pred == (yr[te] == 1))
})
mean(acc)
}
res <- replicate(20, c(leaky = cv_accuracy(FALSE), honest = cv_accuracy(TRUE)))
boxplot(t(res), col = c("firebrick", "seagreen"), ylim = c(0, 1),
ylab = "cross-validated accuracy",
main = "Feature selection outside vs. inside the CV loop")
abline(h = 0.5, lty = 2)
## leaky honest
## mean_accuracy 0.775 0.480
## sd 0.024 0.073The data contain no signal whatsoever, yet selecting features before cross-validating yields accuracies far above chance, while selecting inside each fold correctly returns roughly 0.5. Every preprocessing step that looks at the outcome or at the test data, feature selection, normalization constants, PCA rotations, imputation, harmonization, must be fitted inside the training fold.
Cautions specific to image-derived models.
- Leakage through preprocessing. Fit every step inside the training fold only, and use nested cross-validation when hyperparameters are tuned.
- Overfitting wide matrices. With hundreds of features and tens of patients, spurious predictors abound. Use penalization, stability pre-filtering, and report the events-per-variable ratio.
- Site as a confounder. If scanners differ between classes, the model may learn the scanner. Harmonize, and validate on data from an unseen site.
- Segmentation dependence. Perturb the mask and confirm predictions are stable; a model that flips with a one-voxel boundary change is not clinically usable.
- Optimistic reference standards. If the outcome label was itself derived from the image, the model may simply be re-deriving it.
Checkpoint 10. A radiomics paper reports AUC 0.91 in 60 patients using 12 features chosen from 1,200, with “5-fold cross-validation”. What single question determines whether this is credible? Answer: whether the 12 features were selected inside each training fold. If not, the AUC is uninterpretable, the demonstration above obtains far-above-chance accuracy from pure noise under exactly this protocol.
Volumetric visualization and 3D rendering
Visualization appears at both ends of the pipeline: as inspection during processing and as communication of the final result. Because most clinical data are volumetric, a few standard renderings recur across modalities, and each discards different information, so the choice of view is itself an analytical decision.
Overlays: showing a result in context
A segmentation is most interpretable when displayed on top of the image it came from, so the reader can judge boundary placement directly. Two conventions are standard: a translucent tint and a contour outline. The contour is usually preferable for quality control because it does not obscure the very voxels whose classification is in question.
Both helpers below draw on top of show_img(), which guarantees that the overlay
and the underlying image use the same coordinate convention. Building an
overlay with an independently constructed raster is a frequent source of silently
transposed or vertically flipped figures.
## Both functions assume show_img() has established x = 1..ncol, y = 1..nrow
## with matrix row 1 drawn at the TOP of the plot.
overlay_tint <- function(img, mask, color = "red", alpha = 0.40,
main = "", zlim = c(0, 1)) {
show_img(img, main = main, zlim = zlim)
rgba <- matrix(NA_character_, nrow(mask), ncol(mask))
rgba[mask] <- adjustcolor(color, alpha.f = alpha)
rasterImage(as.raster(rgba), 0.5, 0.5, ncol(mask) + 0.5, nrow(mask) + 0.5,
interpolate = FALSE)
}
overlay_contour <- function(img, mask, color = "red", lwd = 2,
main = "", zlim = c(0, 1), add = FALSE) {
if (!add) show_img(img, main = main, zlim = zlim)
contour(x = 1:ncol(mask), y = 1:nrow(mask),
z = t((mask * 1)[nrow(mask):1, , drop = FALSE]),
levels = 0.5, drawlabels = FALSE, col = color, lwd = lwd, add = TRUE)
}par(mfrow = c(1, 3), mar = c(1, 1, 2.5, 1))
overlay_tint(img_dn, mask_tumour, main = "Tinted overlay")
overlay_contour(img_dn, mask_tumour, main = "Contour overlay")
overlay_contour(img_dn, tumor_mask, color = "green", main = "Estimate vs. truth")
overlay_contour(img_dn, mask_tumour, color = "red", add = TRUE)
legend("bottomleft", c("truth", "estimate"), col = c("green", "red"),
lwd = 2, bty = "n", text.col = "white", cex = 0.9)
The third panel is the figure that a reader of a segmentation paper actually needs: both contours on the same anatomy, at the same window, so that the location of the disagreement, not merely its magnitude, is visible.
A 3D phantom, orthogonal views, and projections
Two display conventions dominate volumetric data.
- Orthogonal (tri-planar) views show the axial, coronal, and sagittal slices through a chosen voxel simultaneously. They preserve quantitative intensity and are the workhorse of radiological reading.
- Maximum-intensity projection (MIP) collapses a slab along the viewing direction by taking the brightest voxel along each ray, \(\text{MIP}(x,y)=\max_z I(x,y,z)\). MIP is ideal for bright, sparse structures such as contrast-filled vessels (CT/MR angiography) or focal uptake in PET, but it discards depth and all intensity information except the maximum, so it systematically over-represents noise spikes and cannot be measured on.
## a small 3D phantom: ellipsoidal "head", spherical lesion, helical vessel
nx <- 64; ny <- 64; nz <- 48
dz <- 2.0 # mm, anisotropic slices
gx <- seq(-1, 1, length.out = nx)
gy <- seq(-1, 1, length.out = ny)
gz <- seq(-1, 1, length.out = nz)
Xv <- array(rep(gx, times = ny * nz), c(nx, ny, nz))
Yv <- array(rep(rep(gy, each = nx), times = nz), c(nx, ny, nz))
Zv <- array(rep(gz, each = nx * ny), c(nx, ny, nz))
head3 <- (Xv / 0.85)^2 + (Yv / 0.85)^2 + (Zv / 0.90)^2 <= 1
lesion3 <- ((Xv - 0.35)^2 + (Yv + 0.20)^2 + (Zv - 0.10)^2) <= 0.16^2
vessel3 <- (abs(Xv - 0.45 * cos(3 * Zv)) < 0.05) &
(abs(Yv - 0.45 * sin(3 * Zv)) < 0.05) & head3
vol <- array(0.05, c(nx, ny, nz))
vol[head3] <- 0.40
vol[vessel3] <- 0.85
vol[lesion3] <- 0.75
set.seed(9)
vol <- vol + array(rnorm(nx * ny * nz, 0, 0.03), c(nx, ny, nz))
c(dim_x = nx, dim_y = ny, dim_z = nz,
voxel_mm = paste(dx, dy, dz, sep = " x "),
voxel_volume_mm3 = dx * dy * dz,
lesion_voxels = sum(lesion3),
lesion_volume_mm3 = round(sum(lesion3) * dx * dy * dz, 1))## dim_x dim_y dim_z voxel_mm
## "64" "64" "48" "1 x 1 x 2"
## voxel_volume_mm3 lesion_voxels lesion_volume_mm3
## "2" "403" "806"ctr <- c(round(nx * 0.67), round(ny * 0.40), round(nz * 0.55))
par(mfrow = c(1, 3), mar = c(1, 1, 2.5, 1))
show_img(t(vol[, , ctr[3]]), sprintf("Axial (z = %d)", ctr[3]), zlim = c(0, 1))
show_img(t(vol[, ctr[2], ]), sprintf("Coronal (y = %d)", ctr[2]), zlim = c(0, 1))
show_img(t(vol[ctr[1], , ]), sprintf("Sagittal (x = %d)",ctr[1]), zlim = c(0, 1))
Anisotropic voxels distort every non-axial view. Our slices are 2 mm thick but 1 mm in plane, so the coronal and sagittal panels above are stretched by a factor of two along the through-plane axis unless the display corrects for it. Any measurement made on such a reformat, a length, an angle, a boundary, must use the physical spacing, not the pixel count.
mip <- function(v, axis) apply(v, setdiff(1:3, axis), max)
par(mfrow = c(1, 3), mar = c(1, 1, 2.5, 1))
show_img(t(mip(vol, 3)), "MIP along z (axial view)", zlim = c(0, 1))
show_img(t(mip(vol, 2)), "MIP along y (coronal view)", zlim = c(0, 1))
show_img(t(mip(vol, 1)), "MIP along x (sagittal view)", zlim = c(0, 1))
par(mfrow = c(1, 1))
## what MIP costs: the projected value is a maximum of many noisy voxels
lesion_true <- 0.75
mip_z <- mip(vol, 3)
data.frame(
quantity = c("true lesion value", "mean over lesion voxels", "MIP over the lesion column"),
value = round(c(lesion_true,
mean(vol[lesion3]),
max(mip_z[round(nx*0.5):round(nx*0.85),
round(ny*0.25):round(ny*0.55)])), 3))## quantity value
## 1 true lesion value 0.750
## 2 mean over lesion voxels 0.752
## 3 MIP over the lesion column 0.937The MIP value exceeds the true lesion intensity because a maximum over \(n_z\) noisy voxels is a biased estimator of the underlying value, the expected maximum of \(m\) draws from \(\mathcal N(\mu,\sigma^2)\) grows roughly as \(\mu + \sigma\sqrt{2\ln m}\). MIP is a detection display, never a quantification one.
Surface and volume rendering
Two families turn a volume into a 3D picture.
Surface rendering first segments a structure, extracts an iso-surface (e.g. by marching cubes), and shades the resulting mesh. It is compact and interactive, but it is only as good as the segmentation, and the iso-value chosen becomes an invisible, unreported parameter that determines the apparent size of everything shown.
Volume rendering instead casts rays through the entire volume, accumulating color and opacity from a transfer function \(\;I \mapsto (r,g,b,\alpha)\) that maps intensity to appearance. Along a ray the standard front-to-back compositing recursion is
\[C_{\text{out}} = C_{\text{in}} + (1 - \alpha_{\text{in}})\,\alpha_i C_i, \qquad \alpha_{\text{out}} = \alpha_{\text{in}} + (1-\alpha_{\text{in}})\,\alpha_i .\]
It shows the full data without committing to a hard segmentation, but it is more computationally demanding and extremely sensitive to the transfer function — which is why two volume renderings of the same CT can suggest quite different lesion sizes.
## a simple ray-cast volume rendering with two different transfer functions
render <- function(v, center, width, gain = 1) {
alpha <- pmin(pmax((v - (center - width/2)) / width, 0), 1) * gain
out <- matrix(0, dim(v)[1], dim(v)[2]); acc <- matrix(0, dim(v)[1], dim(v)[2])
for (k in 1:dim(v)[3]) {
a <- alpha[, , k]
out <- out + (1 - acc) * a * v[, , k]
acc <- acc + (1 - acc) * a
}
out / pmax(acc, 1e-6)
}
par(mfrow = c(1, 2), mar = c(1, 1, 2.5, 1))
show_img(t(render(vol, center = 0.55, width = 0.60, gain = 0.25)),
"Transfer function A (broad, soft)")
show_img(t(render(vol, center = 0.80, width = 0.20, gain = 0.60)),
"Transfer function B (narrow, opaque)")
Both panels render identical data. The transfer function alone changes which structures dominate, the visual analogue of the windowing lesson from the inspection section, and an equally strong argument for reporting the rendering parameters alongside any 3D figure used as evidence.
## A real 3D PET volume from SOCR, rendered as tri-planar orthographic views.
library(brainR); library(httr)
httr::set_config(config(ssl_verifypeer = 0L))
url <- "https://socr.umich.edu/HTML5/BrainViewer/data/PET_FDG_3D_vol.nii.gz"
dest <- file.path(tempdir(), "PET_FDG_3D_vol.nii.gz")
invisible(GET(url, write_disk(dest, overwrite = TRUE)))
pet_vol <- readNIfTI(dest, reorient = FALSE)
mid <- as.integer(dim(pet_vol) / 2 + 0.5)
orthographic(pet_vol, xyz = mid, zlim = range(pet_vol) * 0.9)Single-particle and molecular reconstruction
At the microscopic end of the scale, the same “projections to volume” logic
powers cryo-electron microscopy: thousands of noisy 2D projections of
identical particles in random orientations are aligned, their orientations
estimated, and the images back-projected into a 3D density map from which an
atomic model is built. The mathematics is the central-slice theorem of the
Fourier section, applied with unknown projection angles, and the resulting
structures (distributed as PDB entries, and renderable interactively with
packages such as r3dmol) connect medical image processing directly to
structural biology. Both invert a set of projections to recover an unknown
density; both are limited by SNR, by angular sampling, and by the regularization
implicitly imposed when the data run out.
Pick the rendering to match the question. Use orthogonal slices for quantitative reading and exact localization; MIP for detecting bright, sparse structures such as vessels or focal uptake; surface rendering when a clean segmentation exists and shape or topology is the message; and volume rendering to preview a full dataset without committing to a segmentation. No single view is best, each trades completeness against interpretability, and none of the projective views should ever be measured on.
Checkpoint 11. A vascular surgeon measures a stenosis on a MIP of a CT angiogram and reports 70 % narrowing. Give two reasons the number may be wrong. Answer: MIP takes the maximum along each ray, so an adjacent bright structure (calcification, bone, another vessel) can be projected into the lumen and mask the narrowing; and MIP discards depth, so an oblique vessel is foreshortened and its apparent diameter depends on the projection direction. Measurement should be made on cross-sectional reformats with known spacing.
Capstone: an end-to-end quantitative pipeline
The individual operations in this chapter are most instructive when assembled into a single reproducible workflow. The capstone below takes the synthetic image from acquisition to a calibrated measurement with an uncertainty and a validation score, reusing only functions and ground-truth masks defined earlier. Every stage corresponds to a section above.
## (0) SETUP ----------------------------------------------------------------
stopifnot(exists("img_noisy"), exists("tumor_mask"), exists("measure_lesion"))
truth_area_mm2 <- sum(tumor_mask) * dx * dy
## (1) IMPORT + CALIBRATION -------------------------------------------------
## g is the acquired image; dx, dy come from the metadata (never assumed).
g <- img_noisy
voxel_area_mm2 <- dx * dy
## (2) INSPECT ---------------------------------------------------------------
sigma_hat <- mad(conv2(g, lap_unit)[par_mask]) # robust noise estimate
head_m <- fill_holes(bin_open(g > otsu_threshold(g), struct_elem(2)))
snr_hat <- mean(g[head_m]) / sigma_hat
## (3) PREPROCESS: log-domain polynomial bias correction ---------------------
b_hat <- estimate_bias(g, head_m, degree = 2)
g_bc <- g / b_hat
## (4) DENOISE: edge-preserving 3x3 median (NOT size 1, which is the identity)
g_dn <- median_filter(g_bc, 3)
## (5) SEGMENT: three-class Otsu inside the tissue mask ----------------------
cuts_c <- multi_otsu(g_dn[head_m])
labels <- label_cc(g_dn > cuts_c[2] & head_m)
## (6) SELECT + CLEAN: brightest sizeable component, opened and hole-filled --
mask_c <- fill_holes(bin_open(select_component(labels, g_dn, rule = "brightest"),
struct_elem(1)))
## (7) VALIDATE against the known truth --------------------------------------
val <- c(confusion(mask_c, tumor_mask),
surface_distances(mask_c, tumor_mask, dy, dx))
## (8) MEASURE with calibrated spacing ---------------------------------------
area_c <- sum(mask_c) * voxel_area_mm2
meanI_c <- mean(g_dn[mask_c])
diam_c <- 2 * sqrt(area_c / pi)
## (9) UNCERTAINTY -----------------------------------------------------------
## Type A(i): sensitivity of the reported value to the threshold choice
band_c <- seq(cuts_c[2] - 0.05, cuts_c[2] + 0.05, by = 0.01)
areas_c <- sapply(band_c, function(t) {
m <- select_component(label_cc(g_dn > t & head_m), g_dn, rule = "brightest")
if (!sum(m)) NA else sum(fill_holes(bin_open(m, struct_elem(1)))) * voxel_area_mm2
})
u_threshold <- sd(areas_c, na.rm = TRUE)
## Type A(ii): noise propagated through the whole pipeline (Monte Carlo)
mc_area <- sapply(1:12, function(s) measure_lesion(make_phantom(seed = 2000 + s)$g)$area)
u_noise <- sd(mc_area)
## Type B: geometric calibration, taken as 1% relative on the voxel area
u_calib <- 0.01 * area_c
## combined standard uncertainty and expanded (k = 2) uncertainty
u_c <- sqrt(u_threshold^2 + u_noise^2 + u_calib^2)
U95 <- 2 * u_creport <- data.frame(
quantity = c("SNR (estimated)", "Dice", "Jaccard", "HD95 (mm)", "ASSD (mm)",
"area (mm^2)", "equivalent diameter (mm)", "mean intensity",
"u(threshold) mm^2", "u(noise) mm^2", "u(calibration) mm^2",
"combined u_c (mm^2)", "expanded U, k=2 (mm^2)",
"true area (mm^2)", "bias vs truth (%)"),
value = round(c(snr_hat, val["Dice"], val["Jaccard"], val["HD95"], val["ASSD"],
area_c, diam_c, meanI_c,
u_threshold, u_noise, u_calib, u_c, U95,
truth_area_mm2, 100 * (area_c / truth_area_mm2 - 1)), 3),
row.names = NULL)
report## quantity value
## 1 SNR (estimated) 10.575
## 2 Dice 0.987
## 3 Jaccard 0.975
## 4 HD95 (mm) 1.000
## 5 ASSD (mm) 0.173
## 6 area (mm^2) 541.000
## 7 equivalent diameter (mm) 26.245
## 8 mean intensity 0.746
## 9 u(threshold) mm^2 21.035
## 10 u(noise) mm^2 5.797
## 11 u(calibration) mm^2 5.410
## 12 combined u_c (mm^2) 22.480
## 13 expanded U, k=2 (mm^2) 44.961
## 14 true area (mm^2) 553.000
## 15 bias vs truth (%) -2.170cat(sprintf(
"\nREPORTED RESULT: lesion area = %.0f +/- %.0f mm^2 (k = 2)\n%s\n",
area_c, U95,
sprintf(" equivalent diameter %.1f mm | Dice %.3f, HD95 %.1f mm vs. reference | known bias %+.1f%%",
diam_c, val["Dice"], val["HD95"], 100 * (area_c / truth_area_mm2 - 1))))##
## REPORTED RESULT: lesion area = 541 +/- 45 mm^2 (k = 2)
## equivalent diameter 26.2 mm | Dice 0.987, HD95 1.0 mm vs. reference | known bias -2.2%par(mfrow = c(1, 3), mar = c(1, 1, 2.5, 1))
show_img(g, "1. Acquired image g", zlim = c(0, 1))
show_img(g_dn, "4. Bias-corrected + denoised", zlim = c(0, 1))
overlay_contour(g_dn, tumor_mask, color = "green",
main = "7-9. Result vs. reference")
overlay_contour(g_dn, mask_c, color = "red", add = TRUE)
The whole chapter in one object. The final report is not a single number but a calibrated measurement with an uncertainty and a validation score, produced by a pipeline whose every step is explicit and reproducible. Note also that the expanded uncertainty does not cover the bias: the systematic partial-volume under-estimate is a separate line in the report, because it is corrected by better physics, not by more repeats. That triple, (value, uncertainty, validation), plus a stated bias, is the deliverable that downstream modelling and clinical decision-making actually require.
Instructor’s guide
Suggested two-week module.
| Session | Content | Active-learning artifact |
|---|---|---|
| 1 | Forward model, geometry, calibration (HU/SUV), PSF/FWHM, partial volume | Checkpoints 1–2; recompute the recovery-coefficient curve for a PET-like 5 mm FWHM |
| 2 | Inspection, noise statistics, SNR/CNR/Rose, artifacts | Checkpoint 3; the auto-scaling trap figure as a “spot the error” exercise |
| 3 | Preprocessing, interpolation, bias correction | Checkpoint 4; Activity 2 |
| 4 | Filtering, noise propagation, deconvolution | Checkpoint 5; Activity 1 |
| 5 | Fourier, MTF, Nyquist, Gibbs, FBP | Checkpoint 6; Activity 4 |
| 6 | Registration, metrics, Jacobian | Checkpoint 7; Activity 5 |
| 7 | Segmentation and validation | Checkpoint 8; Activity 3 |
| 8 | Measurement and uncertainty; capstone | Checkpoints 9–10; Activity 6 and the capstone |
Common student misconceptions, and the antidote in this chapter.
- “Smoothing improves the image.” → the noise-propagation table plus the matched-filter curve: smoothing improves detection at one scale and destroys resolution everywhere.
- “A 9×9 mean filter cuts noise 81-fold.” → \(\sigma\lVert w\rVert_2\); the factor is 9.
- “Deconvolution recovers the true image.” → the SNR study: it helps only when there is SNR to spend, and never restores frequencies where MTF ≈ 0.
- “Dice 0.9 means the volume is accurate to 10 %.” → the size-dependence table and the Method A/B comparison in Checkpoint 8.
- “SE = SD/√N.” → the effective-sample-size experiment; the factor is ~7 after modest smoothing.
- “High reproducibility means the measurement is right.” → the test–retest table: CV ≈ 1.7 % with a −3 % bias.
- “Cross-validation protects against overfitting.” → the leakage demonstration reaching 0.78 accuracy on pure noise.
Assessment rubric for a student mini-project (segment and measure a structure in any image of their choosing): (a) geometry and units correct and stated (20 %); (b) preprocessing justified and reported (15 %); (c) segmentation with stated parameters and a cleaning step (15 %); (d) validation with an overlap and a boundary metric (20 %); (e) an uncertainty budget with at least one Type A and one Type B component (20 %); (f) an honest limitations paragraph naming the dominant error source (10 %).
Running the code. The entire teaching thread is base R and knits in well
under two minutes on a laptop. Set run_heavy <- TRUE in the setup chunk to
enable the network-dependent t-SNE and NIfTI demonstrations. The EBImage chunks
self-disable if Bioconductor is unavailable.
Learning activities
These activities use only the synthetic phantom and helper functions built in this chapter, so they are fully reproducible and self-checking against ground truth. Each maps to a section above.
- Denoising trade-offs. Apply mean, Gaussian, median, and bilateral filters
of matching support to
img_bc. For each, report (i) the measured noise SD in the eroded parenchyma, (ii) the predicted SD from \(\sigma\lVert w\rVert_2\) where it applies, (iii) the mean edge gradient across the tumour boundary, and- the Dice of an Otsu segmentation against
tumor_mask. Which filter best preserves the small lesion, and does the ranking change if you evaluate on the tumour instead? Hint: the answer depends on object scale, relate it to the matched-filter curve.
- the Dice of an Otsu segmentation against
- Bias correction and thresholding. Segment the tumour with (a) a fixed
global threshold, (b) an adaptive local threshold, and (c) multi-level Otsu,
each applied to both
img_noisyandimg_bc. Tabulate six Dice values. Explain which combinations fail and why, using the bias field as the source of difficulty. Then repeat with the bias amplitude doubled. - Segmentation validation. Perturb
mask_tumourby dilating and eroding it by 1 and 2 pixels. Plot Dice, HD95, ASSD, and the area error against the perturbation. Which metric is most sensitive to a systematic boundary shift, and which is blind to it? Repeat for the small lesion and explain the difference using the \(1-\text{Dice}\approx\varepsilon/R\) relation. - Fourier filtering and sampling. Apply ideal, Gaussian, and Butterworth low-pass filters at cutoffs of 0.05, 0.10, and 0.20 cycles mm\(^{-1}\). For each, report the overshoot on a clean disk edge and convert the cutoff into a spatial scale in millimeters. Then down-sample the phantom by 2, 3, and 4 with and without an anti-alias pre-filter, and identify at which factor the vessel becomes aliased.
- Registration. Apply a known rigid transform (translation + rotation) to
img_dn, then recover it by maximizing NCC on a coarse grid with parabolic refinement. Report the residual TRE at three landmarks at different radii, and verify the linear growth of TRE with radius. Repeat withmodBas the moving image and confirm that NCC fails while NMI succeeds. - Uncertainty budget. Using
make_phantom()andmeasure_lesion(), build a complete error budget for the lesion area: Monte-Carlo noise term, threshold sensitivity term, and an assumed 1 % calibration term. Then double the noise SD and recompute. Which term dominates in each regime, and what does that imply about where to spend effort, better hardware, or a better algorithm? - Modality interpretation. In a short paragraph, explain why CT, MRI, and PET can share the same processing algorithms yet require different interpretation of intensity, and give one measurement that would be meaningful in one modality and misleading in another.
- Reproduce a published number. Take any quantitative-imaging claim from the literature and complete the reporting checklist for it. Which of the seven items are missing, and could you reproduce the number from what is reported?
Summary
General image processing is the computational bridge between image acquisition and biomedical interpretation. The operations developed here, import and calibration, inspection, preprocessing and normalization, filtering, Fourier analysis, registration, segmentation, measurement, feature extraction, and visualization, are shared across CT, MRI, PET, ultrasound, and microscopy. They are unified by a single object, the image as a calibrated tensor, and by a single equation, the forward model \(g=\mathcal S[b\cdot(h*f)]+\eta\).
Four themes recur.
First, the mathematics is general but the meaning is modality-specific: the same convolution that sharpens a micrograph must be interpreted differently when the intensities are Hounsfield units, relaxation-weighted signal, or tracer concentration, and the same normalization that is mandatory for conventional MRI would destroy the calibration of CT.
Second, every operation has a computable cost. Smoothing reduces noise by \(\lVert w\rVert_2\) and adds resolution in quadrature; sampling below Nyquist folds frequencies irreversibly; truncating k-space rings at 9 % regardless of how much is kept; deconvolution buys resolution with SNR; interpolation blurs, and repeated interpolation blurs repeatedly. These are not rules of thumb but theorems, and each was verified numerically above.
Third, every measurement carries an uncertainty, and accuracy is not precision. A calibrated value with a validation score, an expanded uncertainty, and a stated bias is the real deliverable. Precision comes from repeats; bias comes only from a known reference. Our own pipeline was reproducible to about 1.7 % and biased by about −3 %, and no amount of repetition would have revealed the second number.
Fourth, reproducibility is a requirement, not a virtue: because features propagate from every preprocessing choice into downstream models, the pipeline must be explicit, documented, and testable. Good image analysis therefore demands both computational skill and an understanding of the underlying imaging physics.
Suggested problems
Calibration, geometry, and measurement
- Field of view. An image is \(512\times512\) pixels with pixel spacing \(0.4\text{ mm}\times0.4\text{ mm}\). What is the field of view in millimeters?
- Area. A segmented region contains 2,000 pixels at \(0.25\text{ mm}\) spacing. What is its area in mm\(^2\)?
- Voxel volume. A CT scan has voxels of \(0.7\times0.7\times1.5\text{ mm}\). A lesion contains 8,000 voxels. Give its volume in mm\(^3\) and cm\(^3\).
- Metadata interpretation. In image A a structure is 80 pixels wide at \(0.5\text{ mm}\) spacing; in image B it is 100 pixels wide at \(0.4\text{ mm}\). Are the physical widths the same?
- Oblique geometry. An acquisition has \(R\) = a \(15^\circ\) rotation about \(z\) and \(S=\operatorname{diag}(0.8,0.8,2.5)\). What is the physical voxel volume, and how would you detect a left–right flip from the affine?
- Hounsfield units. Water has \(\mu = 0.019\text{ mm}^{-1}\) at the effective energy used. A voxel has \(\mu = 0.0209\text{ mm}^{-1}\). What is its HU value, and what tissue is that consistent with?
- Resolution and dose. A protocol gives SNR 20 at \(1\times1\times1\text{ mm}\). What SNR results at \(0.5\times0.5\times0.5\text{ mm}\) at the same dose, and by what factor must dose increase to restore SNR 20?
Filtering and convolution
- Kernel interpretation. \(K=\tfrac19\mathbf 1_{3\times3}\). What type of filter is this, and what is its effect on random noise versus sharp edges?
- Normalization. Why must a smoothing kernel’s entries sum to one?
- Noise propagation. By what factor does a \(5\times5\) mean filter reduce white noise? A Gaussian with \(\sigma_g=2\) px?
- Median filter. For the neighborhood \(\{4,5,5,6,7,8,9,100,120\}\), what value does a \(3\times3\) median filter assign to the center pixel, and why does this suppress salt-and-pepper noise?
- Filter choice. A micrograph contains isolated bright hot pixels. Which of mean, Gaussian, or median filtering would you choose, and why?
- Resolution budget. A scanner has a 3.0 mm FWHM PSF. You post-smooth with a 4.0 mm FWHM Gaussian. What is the final resolution?
- Clinical caution. A lesion is only a few pixels wide. Explain why aggressive smoothing could bias a measurement of its size and its mean intensity.
Fourier analysis, sampling, and registration
- Spatial frequency. An image has strong random noise. Would a low-pass or a high-pass filter better reduce it, and what side effect follows?
- Convolution theorem. Explain in words what it means that convolution in image space corresponds to multiplication in frequency space, and give one computational and one conceptual consequence.
- Nyquist. A detail repeats every 3 pixels of 0.5 mm. What is its spatial frequency in cycles mm\(^{-1}\), and what is the largest pixel size that avoids aliasing it?
- MTF. A Gaussian PSF has \(\sigma_h = 1.0\) mm. At what spatial frequency does the MTF fall to 10 %?
- Gibbs. You double the amount of k-space acquired. What happens to the width of the ringing, and what happens to its amplitude?
- Registration metric. You must align a \(T_1\) MRI to a PET scan of the same patient. Why is mutual information preferred over SSD, and why is normalized MI preferred over plain MI?
- TRE. A rigid registration has a residual rotation error of \(0.5^\circ\) about the image center. What is the target registration error at a structure 80 mm from the center?
Segmentation, uncertainty, and modeling
- Quantitative caution. A researcher sharpens an MRI and then measures mean intensity in a segmented region. Why is this problematic?
- Validation. Two segmentations of the same lesion give Dice 0.62 and 0.91 against an expert mask. Interpret these, and state one situation in which a high Dice still hides a clinically important error.
- Perimeter bias. A perfectly circular object is measured with a transition-counting perimeter. What circularity value results, and what is the correction?
- Uncertainty propagation. A volume is \(V=Nv\) with relative uncertainties \(u_N/N=4\%\) and \(u_v/v=3\%\). Estimate the relative uncertainty in \(V\) and the expanded uncertainty at \(k=2\).
- Boundary uncertainty. A spherical lesion of radius 5 mm has a boundary localized to \(\pm0.8\) mm. Estimate the relative volume uncertainty. Repeat for a 20 mm lesion.
- Effective sample size. An ROI contains 900 voxels of a map smoothed with \(\sigma_g=1.5\) px. By what factor does SD\(/\sqrt N\) understate the true standard error?
- Repeatability. Test–retest gives a within-subject SD of 40 mm\(^3\). What is the repeatability coefficient, and what is the smallest change in an individual you could call real?
- RECIST. A lesion’s longest diameter falls by 30 %. By what percentage has its volume fallen, assuming isotropic shrinkage?
- Radiomics reproducibility. A model trained on features from Scanner A fails on Scanner B. List two plausible causes from this chapter and one mitigation for each.
Answers to selected problems
- \(512\times0.4 = 204.8\) mm in each direction.
- \(2000\times0.25\times0.25 = 125\text{ mm}^2\).
- \(8000\times0.7\times0.7\times1.5 = 5880\text{ mm}^3 = 5.88\text{ cm}^3\).
- A: \(80\times0.5 = 40\) mm; B: \(100\times0.4 = 40\) mm, yes, identical.
- \(|\det(RS)| = 0.8\times0.8\times2.5 = 1.6\text{ mm}^3\) (a rotation has \(\det R = 1\), so the rotation does not change the volume). A flip is signalled by \(\det A < 0\).
- \(1000(0.0209-0.019)/0.019 = +100\) HU, consistent with dense soft tissue or contrast-enhanced tissue.
- Voxel volume falls by 8, so SNR falls by \(\sqrt8 = 2.83\) to about 7.1. To restore SNR 20 requires \((20/7.1)^2 \approx 8\times\) the dose.
- A mean (box) low-pass filter: it averages away zero-mean noise and blurs edges, since it cannot distinguish a real step from a fluctuation.
- So the filter has unit DC gain, a uniform region maps to itself, and the mean intensity (and hence any intensity-based measurement) is preserved.
- Mean \(5\times5\): \(\lVert w\rVert_2 = 1/5\), so noise falls by 5. Gaussian \(\sigma_g = 2\): \(\sigma_{\text{out}} \approx \sigma/(2\sigma_g\sqrt\pi) = \sigma/7.09\).
- The median is 7. The two outliers (100, 120) occupy only rank positions, not weight, so they cannot pull the estimate.
- Median, hot pixels are impulse noise; averaging spreads a single 4095-valued pixel over its whole neighborhood instead of removing it.
- \(\sqrt{3.0^2+4.0^2} = 5.0\) mm.
- Smoothing convolves the lesion with the kernel: the apparent size grows (the profile widens) while the peak and mean intensity fall (partial-volume dilution). Both biases worsen as the lesion approaches the kernel width.
- Low-pass, because white noise is spread over all frequencies while the signal is concentrated at low ones; the side effect is loss of edges and fine detail.
- Filtering multiplies the spectrum by the kernel’s transform. Computationally, large-kernel convolution costs \(O(N^2\log N)\) via the FFT instead of \(O(N^2k^2)\); conceptually, “what does this filter do?” is answered by plotting \(|H(f)|\).
- Period \(= 3\times0.5 = 1.5\) mm, so \(f = 0.667\) cycles mm\(^{-1}\). Nyquist requires \(\Delta x \le 1/(2f) = 0.75\) mm.
- \(\text{MTF}=e^{-2\pi^2\sigma_h^2f^2}=0.1 \Rightarrow f=\sqrt{\ln 10/(2\pi^2\sigma_h^2)} = 0.342\) cycles mm\(^{-1}\).
- The ringing becomes narrower (confined closer to the edge) but its fractional overshoot stays at about 8.95 %, the Gibbs constant is independent of the truncation width.
- CT/PET/MRI intensities are related non-monotonically, so SSD and correlation are optimized away from true alignment; MI assumes only statistical dependence. NMI is preferred because plain MI rises as the overlap region shrinks, rewarding the optimizer for sliding the images apart.
- \(\text{TRE} = r\,\delta\theta = 80\times(0.5\times\pi/180) = 0.70\) mm.
- Sharpening is a high-pass boost: it amplifies noise, creates overshoot at every boundary, and changes the intensity distribution inside the ROI, so the “mean intensity” is a property of the filter as much as of the tissue.
- 0.62 is poor (roughly a boundary error of a third of the radius); 0.91 is good. A high Dice can still hide a clinically important error when the disagreement is concentrated in a small but critical region (a margin involving a vessel), or when it is a systematic boundary shift that biases volume while barely moving Dice on a large object.
- Crack perimeter overestimates by \(4/\pi\), so \(4\pi A/P^2\) is underestimated by \((\pi/4)^2 = 0.617\), a perfect circle reads about 0.62. Multiply the transition count by \(\pi/4\) (Cauchy–Crofton), or use the four-direction estimator.
- \(\sqrt{4^2+3^2} = 5\%\); expanded \(U = 2\times5\% = 10\%\).
- \(u_V/V \approx 3\varepsilon/R = 3\times0.8/5 = 48\%\) for the 5 mm lesion, and \(3\times0.8/20 = 12\%\) for the 20 mm one, the same boundary precision gives a fourfold difference in relative volume uncertainty.
- \(2\sigma_g\sqrt\pi = 2\times1.5\times1.772 = 5.3\).
- \(\text{RC} = 2.77\,\sigma_w = 111\text{ mm}^3\). Only a change larger than that in an individual can be attributed to something other than measurement error.
- \((1-0.30)^3 - 1 = -65.7\%\).
- Causes: different reconstruction kernel or slice thickness changing texture features; different intensity scaling or coil bias changing first-order features. Mitigations: resample to a common voxel size and standardize the discretization (IBSI-compliant extraction); harmonize with ComBat-style adjustment and pre-filter to features with high test–retest ICC.
Glossary
| Term | Meaning |
|---|---|
| ASSD | average symmetric surface distance between two segmentation boundaries |
| Bias field | slowly varying multiplicative sensitivity variation |
| CNR | contrast-to-noise ratio, \(\lvert\mu_A-\mu_B\rvert/\sigma\) |
| CRC | contrast recovery coefficient, the partial-volume attenuation of contrast |
| Dice / Jaccard | overlap agreement indices; \(J = D/(2-D)\) |
| FBP | filtered back-projection, ramp-filtered inversion of the Radon transform |
| FWHM | full width at half maximum, \(2\sqrt{2\ln2}\,\sigma\) for a Gaussian |
| GLCM | gray-level co-occurrence matrix, the basis of classical texture features |
| GUM | Guide to the Expression of Uncertainty in Measurement (Type A / Type B) |
| HD95 | 95th-percentile Hausdorff distance |
| HU | Hounsfield unit, \(1000(\mu-\mu_w)/\mu_w\) |
| ICC | intraclass correlation coefficient, between-subject variance fraction |
| LSI | linear shift-invariant; equivalently, describable by a PSF |
| MI / NMI | mutual information; normalized mutual information |
| MIP | maximum-intensity projection |
| MTF | modulation transfer function, \(\lvert H(f)\rvert/\lvert H(0)\rvert\) |
| Nyquist frequency | \(1/(2\Delta x)\), the highest representable frequency |
| PSF | point spread function |
| PVE | partial-volume effect |
| RC | repeatability coefficient, \(2.77\,\sigma_w\) |
| SUV | standardized uptake value |
| TRE / FRE | target / fiducial registration error |
| wCV | within-subject coefficient of variation |
References and further reading
Textbooks and foundational references
- Prince, J.L. & Links, J.M. Medical Imaging Signals and Systems. Pearson. (Forward models, PSF/MTF, sampling, reconstruction.)
- Bushberg, J.T. et al. The Essential Physics of Medical Imaging. Wolters Kluwer. (Modality physics, dose–noise trade-offs, artifacts.)
- Barrett, H.H. & Myers, K.J. Foundations of Image Science. Wiley. (Rigorous treatment of imaging operators, noise, and observer performance.)
- Gonzalez, R.C. & Woods, R.E. Digital Image Processing. Pearson. (Filtering, morphology, Fourier methods, segmentation.)
- Hobbie, R.K. & Roth, B.J. Intermediate Physics for Medicine and Biology. Springer. (The 1D Fourier examples in this chapter follow their treatment.)
- Kak, A.C. & Slaney, M. Principles of Computerized Tomographic Imaging. IEEE Press / SIAM. (Central-slice theorem and filtered back-projection.)
Standards and metrology
- JCGM 100:2008, Evaluation of measurement data, Guide to the expression of uncertainty in measurement (GUM).
- Zwanenburg, A. et al. The Image Biomarker Standardisation Initiative: standardized quantitative radiomics. Radiology, 2020.
- Sullivan, D.C. et al. Metrology standards for quantitative imaging biomarkers. Radiology, 2015. (QIBA framework: bias, repeatability, reproducibility.)
- Bland, J.M. & Altman, D.G. Statistical methods for assessing agreement between two methods of clinical measurement. Lancet, 1986.
- Eisenhauer, E.A. et al. New response evaluation criteria in solid tumours: revised RECIST guideline (version 1.1). Eur. J. Cancer, 2009.
Algorithms
- Otsu, N. A threshold selection method from gray-level histograms. IEEE Trans. Syst. Man Cybern., 1979.
- Tustison, N.J. et al. N4ITK: improved N3 bias correction. IEEE TMI, 2010.
- Maes, F. et al. Multimodality image registration by maximization of mutual information. IEEE TMI, 1997.
- Studholme, C. et al. An overlap invariant entropy measure of 3D medical image alignment. Pattern Recognition, 1999. (Normalized mutual information.)
- Tomasi, C. & Manduchi, R. Bilateral filtering for gray and color images. ICCV, 1998.
- Ronneberger, O. et al. U-Net: convolutional networks for biomedical image segmentation. MICCAI, 2015.
Software
EBImage(Bioconductor): image processing and analysis for R — https://bioconductor.org/packages/EBImage.RNifti,oro.nifti,oro.dicom: NIfTI and DICOM I/O in R.brainR,plotly,r3dmol: 3D and interactive rendering.Rtsne,uwot: t-SNE and UMAP embeddings.- ITK / SimpleITK, ANTs, elastix, 3D Slicer, ITK-SNAP, FSL, FreeSurfer: production registration, segmentation, and analysis toolkits.
- PyRadiomics: an IBSI-aligned radiomic feature library.
SOCR / DSPA companions
- SOCR: Statistics Online Computational Resource, University of Michigan — https://socr.umich.edu.
- Dinov, I.D. Data Science and Predictive Analytics; dimensionality reduction (PCA/ICA, t-SNE, UMAP): DSPA2 Chapter 4.
- Companion chapters in this book: X-ray and CT, MRI, ultrasound, and nuclear medicine (PET/SPECT) for modality physics; the signal-analysis and least-squares chapter for the estimation methods used here.
Closing thought. Across every modality the workflow is the same: turn a physical signal into a calibrated tensor, recover structure from noise, label the structures of interest, measure them with an honest uncertainty, and let those measurements drive a model. The algorithms are general; the physics gives them meaning; the uncertainty makes them usable; and reproducibility makes them trustworthy.