SOCR ≫

BPAD1 Website ≫

BPAD GitHub ≫

Biomedical Physics with Applications to Disease (BPAD1)

Chapter 1: Mathematical and Statistical Foundations of Biomedical Physics

Chapter 1: Mathematical and Statistical Foundations of Biomedical Physics

“The book of nature is written in the language of mathematics.” Galileo Galilei’s claim (The Assayer (Il Saggiatore), 1623) is nowhere more literally true than in biomedical physics, where a photon count, a spin echo, a tracer washout curve, and a survival probability are all read from the same small library of equations.

How to Use This Chapter

This chapter is the mathematical toolbox for the entire textbook. It is designed for three complementary uses.

1. As a course text. Each section develops a tool from its definition, states the properties that matter in practice, and immediately anchors it in a biomedical measurement.
2. As a reference. Sections are self-contained. A reader who arrives from Chapter 4 needing the \(2\)-D Fourier transform can read Section 1.9 alone; a reader arriving from Chapter 6 needing Poisson statistics can read Section 1.15 alone.
3. As classroom material. Every section ends with a Section summary suitable for a lecture slide, and a Checkpoint question suitable for a think-pair-share prompt. Worked examples are set off in their own boxes so they can be projected or assigned.

The conceptual sequence that organizes the chapter is

\[\underbrace{\text{calculus}}_{\text{change}} \rightarrow \underbrace{\text{complex analysis}}_{\text{phase}} \rightarrow \underbrace{\text{linear algebra}}_{\text{structure}} \rightarrow \underbrace{\text{Fourier analysis}}_{\text{frequency}} \rightarrow \underbrace{\text{differential equations}}_{\text{dynamics}} \rightarrow \underbrace{\text{probability and inference}}_{\text{uncertainty}} \rightarrow \underbrace{\text{dimensionality reduction}}_{\text{representation}}.\]

Each arrow is a genuine dependency rather than a curricular convention. Phasors (Section 1.6) require the exponential function of Section 1.5; the Fourier transform (Section 1.9) is built on the complex exponential of Section 1.6 and the integral of Section 1.1; the Nyquist criterion (Section 1.10) is a statement about the Fourier transform; the diffusion equation (Section 1.12) is a partial differential equation whose Green’s function is the Gaussian of Section 1.14; and principal component analysis (Section 1.18) is the eigendecomposition of Section 1.17 applied to the covariance matrix of Section 1.14.

Chapter-level clinical puzzle. A radiologist wants to know how a tracer redistributes in tissue over 60 minutes; a cardiologist wants to know why beat-to-beat ECG intervals fluctuate; an imaging physicist must reconstruct a CT image from projections acquired at 1000 angles; and a data scientist must compress a \(10^5\)-voxel fMRI time series into a handful of interpretable components. What do the four problems have in common? Each reduces to a compartmental ODE, a power spectrum, a Fourier inversion, and an eigendecomposition, respectively, all four of which are developed in this chapter from a common set of principles.

Key idea. The tools in this chapter are not a list of isolated techniques. They are the shared language in which physical, physiological, and statistical phenomena are written. Fluency in that language is what allows a student to move between modalities (optical, acoustic, magnetic, ionizing) and between scales (molecule \(\rightarrow\) cell \(\rightarrow\) tissue \(\rightarrow\) organ \(\rightarrow\) population) without relearning the mathematics each time.

Prerequisites and computational environment

The chapter assumes a first course in single-variable calculus and comfort with elementary algebra. Multivariable calculus, linear algebra, complex analysis, and probability are developed here from their definitions, at the level of rigor needed for the physics that follows, and no further.

All computations use R and are reproducible without any external data or network access. The core figures rely only on base R and ggplot2; a small number of demonstrations optionally use gridExtra (multi-panel layout) and deSolve (numerical ODE solution). Every such use is guarded, so the chapter knits even when the optional packages are absent. Two flags at the top of the source, run_network and run_heavy, control the demonstrations that download real SOCR neuroimaging volumes or run longer computations; both default to FALSE, and the corresponding synthetic demonstrations are mathematically equivalent.

Learning Objectives

By the end of this chapter, the reader should be able to:

#	Objective	Section
1	State the \(\epsilon\)–\(\delta\) definition of a limit and evaluate limits at removable and essential singularities, including by L’Hôpital’s rule.	1.1
2	Differentiate and integrate symbolically and numerically, apply the fundamental theorem of calculus, and locate and classify critical and inflection points.	1.1
3	Expand a function in a Taylor series, bound the truncation error, and derive forward and central finite-difference derivative approximations with their orders of accuracy.	1.2
4	Manipulate scalars, vectors, matrices, and tensors; compute inner products, transposes, traces, determinants, and inverses; and identify the tensor structure of \(2\)-D, \(3\)-D, and \(4\)-D biomedical images.	1.3, 1.17
5	Relate displacement, velocity, and acceleration by differentiation and integration, and solve the projectile problem in closed form.	1.4
6	Compare the asymptotic growth of logarithmic, polynomial, and exponential functions and identify which describes a given biological process.	1.5
7	Perform complex arithmetic in rectangular and polar form, apply Euler’s identity and De Moivre’s theorem, and use phasors to analyze bioimpedance.	1.6
8	Define a linear time-invariant system by its impulse response, compute \(1\)-D and \(2\)-D convolutions, and apply the convolution theorem.	1.7
9	Expand a periodic signal in a Fourier series, interpret the harmonic spectrum of an ECG, and explain the Gibbs phenomenon.	1.8
10	Compute the discrete, continuous, and \(2\)-D Fourier transforms; interpret magnitude and phase spectra; and state Parseval’s and the Wiener–Khinchin theorems.	1.9
11	Apply the Nyquist criterion, predict the aliased frequency of an undersampled signal, and specify anti-aliasing requirements in time and in space.	1.10
12	Convert among Cartesian, polar, cylindrical, and spherical coordinates and select the system that simplifies a given geometry.	1.11
13	Compute partial derivatives, gradients, divergence, curl, and the Laplacian; state the diffusion equation; and verify the divergence theorem on a worked example.	1.12
14	Solve first-order linear ODEs; distinguish exponential, logistic, and Gompertz growth; and analyze one- and two-compartment pharmacokinetic models.	1.13
15	Apply the Kolmogorov axioms, conditional probability, and Bayes’ theorem; compute the first four moments; and select among binomial, Poisson, normal, and multivariate normal models.	1.14
16	Distinguish Gaussian, Poisson, and Rician noise; compute SNR and CNR; and explain the \(\sqrt{N}\) law and its dose implications.	1.15
17	Derive maximum-likelihood estimators, propagate uncertainty through a nonlinear function, construct confidence intervals, and interpret hypothesis tests and their errors.	1.16
18	Solve \(A\mathbf{x} = \mathbf{b}\); compute eigenvalues, eigenvectors, and the SVD; and apply them to diffusion tensors, low-rank denoising, and image registration.	1.17
19	Contrast PCA, ICA, t-SNE, and UMAP and choose an appropriate dimensionality-reduction method for a given biomedical dataset.	1.18
20	Given a biomedical measurement problem, identify the appropriate mathematical tool and justify the choice.	1.19

Notation and Conventions

Consistent notation is used throughout the textbook. Deviations from the table below are always stated explicitly where they occur.

Symbol	Meaning	First used
\(x, y, z, t\)	Cartesian coordinates and time	1.1
\(f'(x),\ \dot{x}\)	derivative with respect to the argument; overdot denotes \(d/dt\)	1.1
\(\partial f/\partial x\)	partial derivative	1.12
\(\nabla f,\ \nabla\cdot\mathbf{F},\ \nabla\times\mathbf{F},\ \nabla^2 f\)	gradient, divergence, curl, Laplacian	1.12
\(\mathbf{a}, \mathbf{b}\) (bold lowercase)	vectors	1.3
\(A, \Sigma\) (uppercase)	matrices; \(\Sigma\) is reserved for covariance	1.3, 1.14
\(A^{T},\ A^{-1},\ \mathrm{tr}(A),\ \det(A)\)	transpose, inverse, trace, determinant	1.3
\(\lambda_k,\ \mathbf{v}_k\)	\(k\)-th eigenvalue and eigenvector	1.17
\(\sigma_k\)	\(k\)-th singular value	1.17
\(i\)	imaginary unit, \(i^2 = -1\)	1.6
\(z = x + iy = \lvert z\rvert e^{i\theta}\)	complex number in rectangular and polar form	1.6
\(\tilde{A} = A e^{i\phi}\)	phasor (amplitude \(A\), phase \(\phi\))	1.6
\((x * h)(t)\)	convolution	1.7
\(\mathcal{F}\{\cdot\},\ \hat{f}(\omega)\)	Fourier transform	1.9
\(f_s,\ f_{\mathrm{Nyq}} = f_s/2\)	sampling frequency, Nyquist frequency	1.10
\(k\)-space	spatial-frequency domain (MRI usage)	1.9
\(P(A),\ P(A\mid B)\)	probability, conditional probability	1.14
\(\mathbb{E}[X],\ \mathrm{Var}(X)\)	expectation, variance	1.14
\(\mathcal{N}(\mu,\sigma^2)\)	normal distribution	1.14
\(\hat{\theta}\)	estimator of the parameter \(\theta\)	1.16
\(\mathrm{SNR},\ \mathrm{CNR}\)	signal-to-noise, contrast-to-noise ratio	1.15
\(\mathbb{R}, \mathbb{C}, \mathbb{Z}, \mathbb{N}_0\)	reals, complex numbers, integers, non-negative integers	1.1

Equations are numbered \((1.n)\) within this chapter and are referenced across the book as “Eq. (1.n)”. Figures are numbered “Figure 1.n”.

Roadmap: Where Each Tool Reappears

The table below is the chapter’s contract with the rest of the textbook. Each row states a mathematical construct developed here, the chapter in which it becomes physics, and a concrete biomedical instance. It is worth returning to this table after each subsequent chapter.

Foundation (this chapter)	Where it becomes physics	Concrete biomedical instance
Exponential decay, first-order ODE ((13))	Ch. 2, 5, 6	Beer–Lambert attenuation; fluorescence lifetime; radiotracer decay; drug clearance
Definite integration and optimization ((1))	Ch. 2	Stefan–Boltzmann law as an integral of the Planck function; Wien’s law as a critical point
Fourier transform ((9))	Ch. 2, 4, 5, 6, 7	FTIR interferogram \(\rightarrow\) spectrum; MRI \(k\)-space \(\rightarrow\) image; central slice theorem in CT
Wave equation and diffusion equation ((12))	Ch. 3, 4	Photoacoustic pressure waves; apparent diffusion coefficient in stroke imaging
Convolution and the PSF ((7))	Ch. 3, 5, 6, 7, 8	Ultrasound beam profile; CT detector blur; CNN convolution layers
Nyquist sampling ((10))	Ch. 3, 4, 7	Doppler velocity aliasing; MRI wrap-around; digital pathology pixel size
Complex numbers and phasors ((6))	Ch. 2, 4	Quadrature MR signal; low-coherence interferometry in OCT; bioimpedance
Poisson statistics and the \(\sqrt{N}\) law ((15))	Ch. 5, 6	CT dose–noise trade-off; PET count statistics and NECR
Rician statistics ((15))	Ch. 4	Magnitude-image noise bias at low SNR in diffusion MRI
Eigenvalues and eigenvectors ((17))	Ch. 4, 7, 8	Diffusion-tensor eigenvalues and fractional anisotropy; shape modes
SVD and low-rank structure ((17))	Ch. 6, 7, 8	Image denoising; accelerated reconstruction; radiomic feature compression
Bayes’ theorem ((14))	Ch. 6, 8	MAP reconstruction in PET; converting sensitivity and specificity into post-test probability
Maximum likelihood ((16))	Ch. 6, 8	MLEM reconstruction in emission tomography; logistic regression
Error propagation ((16))	Ch. 7, 8	Uncertainty in a segmented tumor volume; radiomic feature stability
Compartmental ODEs ((13))	Ch. 6, 8	FDG two-tissue kinetic model; Patlak analysis; longitudinal tumor growth
PCA / ICA / t-SNE / UMAP ((18))	Ch. 7, 8	fMRI network separation; radiomic feature reduction; cohort visualization

1.1 Single-Variable Calculus

Calculus is the mathematics of continuous change. The derivative measures an instantaneous rate of change; the integral measures an accumulated total change. The fundamental theorem of calculus states that the two operations are inverses of one another. Together they supply the grammar for every physical law in this textbook: Newton’s second law, the Bloch equations of Chapter 4, the Beer–Lambert law of Chapter 2, the decay law of Chapter 6, and the diffusion equation of Section 1.12 are all statements about derivatives, and every quantity that a detector actually reports (a photon count, an area under a curve, an average intensity) is an integral.

1.1.1 Limits and continuity

Definition 1.1 (Limit). The function \(f\) has limit \(L\) as \(x \to a\), written \(\lim_{x\to a} f(x) = L\), if for every \(\epsilon > 0\) there exists \(\delta_\epsilon > 0\) such that

\[\begin{equation} 0 < \lvert x - a\rvert < \delta_\epsilon \quad \Longrightarrow \quad \lvert f(x) - L\rvert < \epsilon. \tag{1} \end{equation}\]

Read operationally: whatever tolerance \(\epsilon\) you demand on the output, I can find a tolerance \(\delta_\epsilon\) on the input that guarantees it. The strict inequality \(0 < \lvert x-a\rvert\) excludes the point \(a\) itself, which is precisely why a limit can exist where the function is undefined.

Two consequences are used constantly:

• Linearity. \(\lim_{x\to a}\left[\lambda f(x) + \gamma g(x)\right] = \lambda \lim_{x\to a} f(x) + \gamma \lim_{x\to a} g(x)\) whenever both limits exist.
• Continuity is an extra assumption. \(f\) is continuous at \(a\) when \(\lim_{x\to a} f(x) = f(a)\). Without continuity, \(f(a)\) may be undefined or unrelated to the limit.

The distinction is not pedantic. Consider

\[\begin{equation} f_1(t) = \frac{1}{t}, \qquad f_2(t) = \begin{cases} 20\left(1 + \dfrac{1}{t+1}\right), & \lvert t+1\rvert \ge 0.1,\\[6pt] 0, & \lvert t+1\rvert < 0.1. \end{cases} \tag{2} \end{equation}\]

Both have a pole, \(f_1\) at \(t = 0\) and \(f_2\) at \(t = -1\). The function \(f_1\) is genuinely unbounded near the origin: the two one-sided limits are \(\pm\infty\) and no limit exists. The piecewise definition of \(f_2\) regularizes the singularity by replacing the divergent neighborhood with a finite value. This is exactly what a detector does when its output saturates, and exactly what a reconstruction algorithm does when it clips a divergent ratio. The regularized function is bounded, but it is still discontinuous, and the discontinuity has consequences downstream (Section 1.8, Gibbs phenomenon).

t  <- seq(-pi, pi, length.out = 1201)
f1 <- 1/t
f2 <- ifelse(abs(t + 1) < 0.1, 0, 20 * (1 + 1/(t + 1)))
f1[abs(f1) > 60] <- NA          # suppress the divergent branch for plotting
f2[abs(f2) > 60] <- NA

lim_df <- rbind(
  data.frame(t = t, y = f1, fn = "f1(t) = 1/t"),
  data.frame(t = t, y = f2, fn = "f2(t): regularized at t = -1"))

ggplot(lim_df, aes(t, y, colour = fn)) +
  geom_line(linewidth = 0.7, na.rm = TRUE) +
  geom_vline(xintercept = c(-1, 0), linetype = "dashed", colour = "grey55") +
  scale_colour_manual(values = bpad_pal[1:2]) +
  coord_cartesian(ylim = c(-60, 60)) +
  labs(title = "Behaviour of two functions near their singular points",
       subtitle = "Dashed vertical lines mark t = -1 and t = 0",
       x = "t", y = "function value")

Figure 1: Figure 1.1: Two functions near their singular points. The unregularized 1/t diverges at t = 0; the piecewise-defined f2 is bounded near t = -1 but remains discontinuous there. Values beyond +/-60 are suppressed for display.

Checkpoint 1.1. Is \(\lim_{t\to 0}\sin(t)/t = 1\) even though the expression is undefined at \(t = 0\)? Justify your answer using Definition 1.1, and then confirm it with L’Hôpital’s rule.

1.1.2 Derivatives

Definition 1.2 (Derivative). The derivative of \(f\) at \(x = a\) is

\[\begin{equation} f'(a) \equiv \lim_{x\to a}\frac{f(x) - f(a)}{x - a} = \lim_{h\to 0}\frac{f(a+h)-f(a)}{h}, \tag{3} \end{equation}\]

whenever the limit exists. Geometrically, \(f'(a)\) is the slope of the tangent line at \(a\). Physically, if \(f(t)\) is a position then \(f'(t)\) is a velocity; if \(f(t)\) is a concentration then \(f'(t)\) is a flux; if \(f(E)\) is a photon count then \(f'(E)\) is a spectral density.

The standard rules, with the argument suppressed, are

Rule	Statement
Constants	\(c' = 0\) for every constant \(c\in\mathbb{C}\)
Linearity	\((\alpha f + \beta g)' = \alpha f' + \beta g'\)
Power	\((x^n)' = n x^{n-1}\), \(n \in \mathbb{Z}\) (and, more generally, \(n\in\mathbb{R}\) for \(x>0\))
Product	\((fg)' = f'g + fg'\)
Quotient	\((f/g)' = (f'g - fg')/g^2\)
Chain	\(\big(f(g(x))\big)' = f'(g(x))\, g'(x)\)
Trigonometric	\(\sin' = \cos\), \(\cos' = -\sin\)
Exponential	\((e^x)' = e^x\)
Logarithm	\((\ln x)' = 1/x\), \(x>0\)
Higher order	\((f')' \equiv f''\), \((f^{(n-1)})' \equiv f^{(n)}\)

Two of these deserve emphasis. The exponential is the unique function (up to scaling) that equals its own derivative, which is the reason it appears as the solution of every first-order linear rate law in Section 1.13. And Euler’s number is itself defined by a limit,

\[e = \lim_{n\to\infty}\left(1 + \frac{1}{n}\right)^{n} = 2.718281828\ldots\]

L’Hôpital’s rule. If \(f\) and \(g\) are differentiable on an open interval \(I \ni a\), \(g'(x)\ne 0\) on \(I\), and either \(\lim_{x\to a}f = \lim_{x\to a}g = 0\) or \(\lim_{x\to a}f = \lim_{x\to a}g = \pm\infty\), then

\[\begin{equation} \lim_{x\to a}\frac{f(x)}{g(x)} = \lim_{x\to a}\frac{f'(x)}{g'(x)}, \tag{4} \end{equation}\]

provided the right-hand limit exists. This resolves the \(0/0\) and \(\infty/\infty\) forms that arise whenever a ratio of two physical quantities is evaluated at a boundary, for example the low-dose limit of a contrast-to-noise ratio, or \(\sin(x)/x\) at the origin (the sinc function of Section 1.10).

Critical points and concavity. Interior extrema of a differentiable function occur only where \(f'(x) = 0\); such points are critical points. The second derivative classifies them:

• \(f''(x) < 0\): locally concave down, relative maximum;
• \(f''(x) > 0\): locally concave up, relative minimum;
• \(f''(x) = 0\) with a sign change: inflection point, where the curvature reverses and where \(f'\) itself attains an extremum.

In optimization language, \(f\) is a cost or loss function, and in the multivariable case (Section 1.12) the role of \(f''\) is played by the positive-definiteness of the Hessian matrix. This is the machinery behind every model fit in Chapter 8.

Mean value theorem. If \(f\) is continuous on \([a,b]\) and differentiable on \((a,b)\), there exists \(c\in(a,b)\) with

\[\begin{equation} f'(c) = \frac{f(b)-f(a)}{b-a}. \tag{5} \end{equation}\]

The theorem is the formal guarantee that an average rate is attained as an instantaneous rate somewhere in the interval, the assumption implicitly made whenever a clinical measurement reports a mean washout rate or a mean growth rate.

Checkpoint 1.2. For \(f(x) = x^3 - 5x\), locate the critical points and classify each. Where is the inflection point, and why is it not a critical point?

1.1.3 Integrals

Antiderivatives. \(F\) is an antiderivative of \(f\) if \(F' = f\). All antiderivatives differ by a constant, so the indefinite integral is

\[\begin{equation} \int f(t)\,dt = F(t) + c. \tag{6} \end{equation}\]

Consider the elementary growth law \(dN/dt = kN\), in which the rate of change of a population, a bacterial colony, or a tumour cell count is proportional to the current size. Because \((e^{t})' = e^{t}\), the solution is

\[\begin{equation} N(t) = N_0 e^{kt}, \tag{7} \end{equation}\]

as verified directly: \(\left(N_0 e^{kt}\right)' = k N_0 e^{kt} = kN(t)\). Equation (7) recurs, with different names for \(k\), throughout the book: as the attenuation law in Chapter 5, as the decay law in Chapter 6, and as the relaxation law in Chapter 4.

Definite integrals. For \(f\) integrable on \([a,b]\),

\[\begin{equation} \int_a^b f(x)\,dx = \lim_{n\to\infty}\sum_{i=1}^{n} f(x_i^{*})\,\Delta x, \qquad \Delta x = \frac{b-a}{n}, \tag{8} \end{equation}\]

where \(x_i^{*}\) is any point in the \(i\)-th subinterval. The limit is the Riemann integral, and each finite sum is a Riemann sum. Different choices of \(x_i^{*}\) give different approximations that converge to the same limit but at very different rates, a fact of considerable practical importance because every numerical integration in imaging (an area under a time-activity curve, a total counted dose, a projection line integral) is a Riemann sum in disguise.

f <- function(x) sin(x) * (cos(2*x) + x)
a <- -pi/4; b <- 9*pi/8
exact <- integrate(f, a, b)$value

riemann <- function(n, rule = c("mid", "left", "right")) {
  rule <- match.arg(rule); h <- (b - a)/n
  xs <- switch(rule,
               left  = a + (0:(n-1))*h,
               mid   = a + (0:(n-1))*h + h/2,
               right = a + (1:n)*h)
  sum(f(xs)) * h
}

ns  <- c(6, 15, 51, 500)
tab <- data.frame(
  n              = ns,
  left           = sapply(ns, riemann, rule = "left"),
  midpoint       = sapply(ns, riemann, rule = "mid"),
  right          = sapply(ns, riemann, rule = "right"),
  midpoint_error = abs(sapply(ns, riemann, rule = "mid") - exact))
knitr::kable(tab, digits = 5,
             caption = paste0("Riemann-sum approximations to the exact value ",
                              round(exact, 5), "."))

Table 1: Riemann-sum approximations to the exact value 2.16476.

n	left	midpoint	right	midpoint_error
6	2.90365	2.18683	1.33527	0.02207
15	2.47098	2.16847	1.84363	0.00371
51	2.25637	2.16508	2.07185	0.00032
500	2.17416	2.16476	2.15534	0.00000

panels <- lapply(ns, function(n) {
  h    <- (b - a)/n
  left <- a + (0:(n-1))*h
  bars <- data.frame(xmin = left, xmax = left + h,
                     ymin = 0,    ymax = f(left + h/2))
  crv  <- data.frame(x = seq(a, b, length.out = 400))
  crv$y <- f(crv$x)
  ggplot() +
    geom_rect(data = bars, aes(xmin = xmin, xmax = xmax, ymin = ymin, ymax = ymax),
              fill = bpad_pal[1], alpha = 0.35, colour = "white", linewidth = 0.2) +
    geom_line(data = crv, aes(x, y), colour = bpad_pal[2], linewidth = 0.8) +
    labs(title = sprintf("n = %d,  sum = %.4f", n, riemann(n)), x = NULL, y = NULL)
})
bpad_grid(panels[[1]], panels[[2]], panels[[3]], panels[[4]], ncol = 2)

Figure 2: Figure 1.2: Riemann sums converging to a definite integral. Rectangle heights are taken at subinterval midpoints. The approximation error falls by a factor of about four each time n doubles, the signature of second-order accuracy.

Fundamental theorem of calculus. Differentiation and integration are inverse operations:

\[\begin{equation} \int_a^b f'(x)\,dx = f(b) - f(a), \tag{9} \end{equation}\]

\[\begin{equation} \frac{d}{dx}\int_a^{x} f(y)\,dy = f(x). \tag{10} \end{equation}\]

When the limits themselves are functions, the chain rule gives the Leibniz rule

\[\begin{equation} \frac{d}{dx}\int_{h(x)}^{g(x)} f(y)\,dy = f\big(g(x)\big)\,g'(x) - f\big(h(x)\big)\,h'(x). \tag{11} \end{equation}\]

Note the dummy variable \(y\) inside the integral: the definite integral is a function of its limits, not of the integration variable. Equation (11) is the tool used in Chapter 6 to differentiate a time-activity curve whose integration window itself moves.

Integrals as accumulation. If \(F'(t)\) is the rate of change of a quantity \(F\), then

\[\begin{equation} \int_a^b F'(t)\,dt = F(b) - F(a) \tag{12} \end{equation}\]

is the total change over \([a,b]\). This single identity underlies the area under the concentration–time curve (AUC) in pharmacokinetics, the cumulative dose in radiotherapy, and the total counts in a PET frame.

Geometric quantities from integrals. Beyond areas, integrals deliver lengths and volumes:

• arc length of \(y = f(x)\) over \([a,b]\): \(\displaystyle\int_a^b \sqrt{1 + \big(f'(x)\big)^2}\;dx\);
• volume of revolution about the \(x\)-axis: \(\displaystyle\int_a^b \pi \big(f(x)\big)^2\,dx\) (disk method);
• volume of revolution about the \(y\)-axis: \(\displaystyle\int_a^b 2\pi x\, f(x)\,dx\) (shell method).

The disk method is exactly how a tumour volume is estimated from a stack of segmented cross-sections in Chapter 7.

Improper integrals. When the interval is unbounded or the integrand is singular, the integral is defined as a limit, for instance \(\int_0^\infty e^{-x}dx = \lim_{b\to\infty}\int_0^b e^{-x}dx = 1\). Improper integrals appear whenever a total exposure is computed over all time, as in \(\mathrm{AUC}_{0\to\infty}\) in Section 1.13.

1.1.4 Worked example: growth rate of Coho salmon

Worked Example 1.1 (average value of a function). Fisher (1963) and Warren (1971) modelled the growth rate of juvenile Coho salmon as a function of dissolved oxygen concentration \(x\) (mg L\(^{-1}\)). Under an unrestricted ration the growth rate (mg g\(^{-1}\) day\(^{-1}\)) is

\[y_u(x) = 7.3\,(x + 3.5)\,e^{-0.05x},\]

while under a restricted ration it is essentially constant, \(y_r(x) = 28\). Compare the two diets by their average growth rates over the physiological range \(3 \le x \le 30\).

The average value of a function on \([a,b]\) is \(\bar{y} = \frac{1}{b-a}\int_a^b y(x)\,dx\). For the restricted ration this is immediate:

\[\bar{y}_r = \frac{1}{27}\int_3^{30} 28\,dx = \frac{28 \times 27}{27} = 28.\]

For the unrestricted ration, split the integrand and use \(7.3\times 3.5/27 = 0.946\) and \(7.3/27 = 0.270\):

\[\bar{y}_u = \underbrace{0.946\int_3^{30} e^{-0.05x}dx}_{I} + \underbrace{0.270\int_3^{30} x\,e^{-0.05x}dx}_{II}.\]

The first integral is elementary:

\[I = 0.946\left[-\frac{e^{-0.05x}}{0.05}\right]_3^{30} = 18.93\left(e^{-0.15} - e^{-1.5}\right) = 18.93(0.8607 - 0.2231) = 12.07.\]

The second requires integration by parts, \(\int u\,dv = uv - \int v\,du\), with \(u = x\) and \(dv = e^{-0.05x}dx\), so \(v = -e^{-0.05x}/0.05\):

\[\int_3^{30} x e^{-0.05x}dx = \left[-\frac{x e^{-0.05x}}{0.05}\right]_3^{30} + \frac{1}{0.05}\int_3^{30} e^{-0.05x}dx = \left[-(20x + 400)e^{-0.05x}\right]_3^{30} = 172.80,\]

hence \(II = 0.270 \times 172.80 = 46.72\) and

\[\bar{y}_u = I + II = 12.07 + 46.72 = 58.79\ \text{mg g}^{-1}\text{day}^{-1}.\]

The unrestricted diet therefore supports more than double the average growth rate of the restricted diet, even though the two curves nearly coincide at the low-oxygen end of the range. Averaging over the whole physiological range, rather than comparing at a single operating point, is what exposes the difference.

y_u <- function(x) 7.3*(x + 3.5)*exp(-0.05*x)
xs  <- seq(3, 30, length.out = 400)

avg_u <- integrate(y_u, 3, 30)$value / 27
avg_r <- 28
cat(sprintf("Average growth rate, unrestricted ration = %.4f mg/g/day\n", avg_u))

## Average growth rate, unrestricted ration = 58.7855 mg/g/day

cat(sprintf("Average growth rate, restricted ration   = %.4f mg/g/day\n", avg_r))

## Average growth rate, restricted ration   = 28.0000 mg/g/day

cat(sprintf("Hand calculation in Worked Example 1.1   = 58.79 (rounding only)\n"))

## Hand calculation in Worked Example 1.1   = 58.79 (rounding only)

df <- data.frame(x = xs, unrestricted = y_u(xs), restricted = 28)
ggplot(df, aes(x)) +
  geom_ribbon(aes(ymin = 0, ymax = unrestricted), fill = bpad_pal[1], alpha = 0.18) +
  geom_line(aes(y = unrestricted, colour = "unrestricted ration"), linewidth = 0.9) +
  geom_line(aes(y = restricted,   colour = "restricted ration"),   linewidth = 0.9) +
  geom_hline(yintercept = c(avg_u, avg_r), linetype = "dashed", colour = "grey40") +
  annotate("text", x = 28, y = avg_u + 3, label = sprintf("mean = %.1f", avg_u),
           size = 3.2, hjust = 1) +
  annotate("text", x = 28, y = avg_r - 3, label = sprintf("mean = %.1f", avg_r),
           size = 3.2, hjust = 1) +
  scale_colour_manual(values = bpad_pal[1:2]) +
  labs(title = "Growth rate of juvenile Coho salmon",
       x = "dissolved oxygen (mg/L)", y = "growth rate (mg/g/day)")

Figure 3: Figure 1.3: Coho salmon growth rate versus dissolved oxygen under unrestricted and restricted rations. Horizontal dashed lines mark the two average growth rates over the interval [3, 30] mg/L; the shaded region is the integral evaluated in Worked Example 1.1.

1.1.5 Worked example: distance as the integral of velocity

Worked Example 1.2. Velocity is the rate of change of position, so the distance travelled is the accumulated velocity:

\[D = \int_{t_0}^{t_1} v(t)\,dt = \int_{t_0}^{t_1} dx = x(t_1) - x(t_0).\]

For an object in free fall with \(v(t) = gt\) (starting from rest, \(g\) the gravitational acceleration),

\[\begin{equation} D = \int_0^{T} g t\,dt = g\left.\frac{t^2}{2}\right|_0^T = \tfrac{1}{2}gT^2. \tag{13} \end{equation}\]

The familiar quadratic displacement law is nothing more than the integral of a constant acceleration, performed twice. Section 1.4 makes this chain explicit.

1.1.6 Symbolic and numerical calculus in R

Base R provides symbolic differentiation through D() (single variable, returns an expression) and deriv() (multivariable, optionally returns an executable function with a gradient attribute). Neither requires an add-on package, so the code below runs in any R installation. The package Deriv generalizes both and can be substituted where available.

f1_expr <- expression(x^3 - 5*x)

d1 <- D(f1_expr, "x"); d1     # f'(x)  = 3x^2 - 5

## 3 * x^2 - 5

d2 <- D(d1,      "x"); d2     # f''(x) = 6x

## 3 * (2 * x)

## A general n-th order derivative by recursion
nth_derivative <- function(expr, var, n = 1) {
  if (n < 1)  stop("derivative order must be >= 1")
  if (n == 1) D(expr, var) else nth_derivative(D(expr, var), var, n - 1)
}
nth_derivative(f1_expr, "x", 3)          # third derivative: the constant 6

## 3 * 2

## Evaluate an expression at chosen points
x <- c(1, 2, 3)
data.frame(x = x, f = eval(f1_expr), f_prime = eval(d1), f_double_prime = eval(d2))

##   x  f f_prime f_double_prime
## 1 1 -4      -2              6
## 2 2 -2       7             12
## 3 3 12      22             18

## Partial derivatives of f2(x, y) = sin(cos(5x + y^3))
f2_expr <- expression(sin(cos(5*x + y^3)))
D(f2_expr, "x")

## -(cos(cos(5 * x + y^3)) * (sin(5 * x + y^3) * 5))

D(f2_expr, "y")

## -(cos(cos(5 * x + y^3)) * (sin(5 * x + y^3) * (3 * y^2)))

## deriv() returns a function whose "gradient" attribute holds both partials
g <- deriv(f2_expr, c("x", "y"), func = TRUE)
g(0.3, 0.7)

## [1] -0.2656274
## attr(,"gradient")
##              x         y
## [1,] -4.642896 -1.365011

f  <- function(x) x^3 - 5*x
fp <- function(x) 3*x^2 - 5
fq <- function(x) 6*x

## Definite integration, and comparison with a discrete trapezoid rule
integrate(f, lower = -5, upper = 8)

## 770.25 with absolute error < 1.1e-11

xs <- seq(-5, 8, by = 0.01); fx <- f(xs)
sum(diff(xs) * (head(fx, -1) + tail(fx, -1)) / 2)   # trapezoid rule

## [1] 770.251

crit <- c(-sqrt(5/3), sqrt(5/3))
xg   <- seq(-3, 3, length.out = 400)
dd   <- rbind(data.frame(x = xg, y = f(xg),  fn = "f(x) = x^3 - 5x"),
              data.frame(x = xg, y = fp(xg), fn = "f'(x) = 3x^2 - 5"),
              data.frame(x = xg, y = fq(xg), fn = "f''(x) = 6x"))
ggplot(dd, aes(x, y, colour = fn)) +
  geom_hline(yintercept = 0, colour = "grey70") +
  geom_line(linewidth = 0.8) +
  geom_point(data = data.frame(x = crit, y = f(crit)),
             aes(x, y), inherit.aes = FALSE, shape = 21, size = 3,
             fill = "white", colour = "black") +
  scale_colour_manual(values = bpad_pal[1:3]) +
  labs(title = "A cubic and its first two derivatives",
       x = "x", y = "value")

Figure 4: Figure 1.4: A cubic and its first two derivatives. Critical points of f (open circles) sit where f’ crosses zero; the inflection point sits where f’’ crosses zero.

Higher-dimensional integrals can be built by nesting integrate(). The helper below computes a double integral and is verified against a case with a known closed form, the total mass of a two-dimensional diffusion profile, which anticipates Section 1.12.

integrate2 <- function(f, xlim, ylim) {
  inner <- function(y) sapply(y, function(yy)
    integrate(function(x) f(x, yy), xlim[1], xlim[2])$value)
  integrate(inner, ylim[1], ylim[2])$value
}

D_coef <- 1e-3; t_obs <- 1
C <- function(x, y) exp(-(x^2 + y^2)/(4*D_coef*t_obs))   # unnormalized 2-D Gaussian

cat("numerical double integral :", integrate2(C, c(-1, 1), c(-1, 1)), "\n")

## numerical double integral : 0.01256637

cat("analytic value 4*pi*D*t   :", 4*pi*D_coef*t_obs, "\n")

## analytic value 4*pi*D*t   : 0.01256637

For genuinely multidimensional problems, cubature::adaptIntegrate() is the standard tool; for discrete data (as opposed to analytic functions), the discrete analogues of differentiation and integration are diff() and cumsum().

1.1.7 Optimization, critical points, and the logistic model

The first and second derivatives together locate three special classes of points: relative maxima, relative minima, and inflection points. The word relative matters, because a curve can rise or fall further elsewhere in its domain.

The logistic growth model provides the canonical biomedical illustration. Writing the population size as

\[N(t) = \frac{N_0(1+b)}{1 + b\,e^{-kt}} = N_0(1+b)\left(1 + b e^{-kt}\right)^{-1}, \qquad k > 0,\ b > 1,\]

the carrying capacity is \(N_\infty = N_0(1+b)\) and the growth rate is

\[\frac{dN}{dt} = N_0(1+b)\,b\,k\,e^{-kt}\left(1 + be^{-kt}\right)^{-2}.\]

The maximum growth rate occurs where the derivative of the growth rate vanishes, that is at the inflection point of \(N(t)\):

\[0 = \frac{d^2N}{dt^2} = \underbrace{N_0(1+b)bk^2 e^{-kt_1}\left(1+be^{-kt_1}\right)^{-3}}_{>\,0}\underbrace{\left(be^{-kt_1} - 1\right)}_{=\,0},\]

so \(b e^{-kt_1} = 1\), giving

\[\begin{equation} t_1 = \frac{\ln b}{k}, \qquad N_1 = N(t_1) = \frac{N_0(1+b)}{2} = \frac{N_\infty}{2}. \tag{14} \end{equation}\]

N0 <- 1; b <- 60; k <- 0.25
Ninf <- N0*(1 + b)
Nfun <- function(t) N0*(1 + b)/(1 + b*exp(-k*t))
dN   <- function(t) N0*(1 + b)*b*k*exp(-k*t)/(1 + b*exp(-k*t))^2
t1   <- log(b)/k

tg <- seq(0, 45, length.out = 600)
scale_rate <- max(Nfun(tg))/max(dN(tg))
ggplot(data.frame(t = tg), aes(t)) +
  geom_line(aes(y = Nfun(t),           colour = "N(t)"),               linewidth = 0.9) +
  geom_line(aes(y = dN(t)*scale_rate,  colour = "dN/dt (rescaled)"),   linewidth = 0.9) +
  geom_vline(xintercept = t1, linetype = "dashed", colour = "grey40") +
  geom_hline(yintercept = c(Ninf, Ninf/2), linetype = "dotted", colour = "grey60") +
  annotate("text", x = t1 + 0.6, y = Ninf*0.25, hjust = 0, size = 3.2,
           label = sprintf("inflection: t = ln(b)/k = %.1f,  N = N_inf/2 = %.0f", t1, Ninf/2)) +
  scale_colour_manual(values = bpad_pal[1:2]) +
  labs(title = "Logistic growth and its rate", x = "time", y = "population size")

Figure 5: Figure 1.5: Logistic growth (left axis) and its rate (right, rescaled). The population grows fastest exactly when it reaches half of the carrying capacity, the inflection point of N(t).

From equation to population biology and oncology. The result \(N_1 = N_\infty/2\), that the maximum growth rate occurs at half the carrying capacity, is the mathematical basis of maximum-sustainable-yield policy in fisheries management, and its oncological mirror image: a tumour held near half its carrying capacity is a tumour growing at its fastest possible rate. This is one reason cytostatic therapies that stabilize tumour size at an intermediate volume can be counterproductive, and why the Gompertz model of Section 1.13, whose inflection occurs at \(N_\infty/e \approx 0.37 N_\infty\), often fits observed tumour-volume curves better.

Section 1.1 summary.

• The derivative is an instantaneous rate; the integral is an accumulated total; the fundamental theorem makes them inverses.
• Limits, continuity, and L’Hôpital’s rule govern local behaviour, including at singular points that real detectors must regularize.
• Critical points (\(f' = 0\)) locate extrema; the sign of \(f''\) classifies them; inflection points (\(f'' = 0\) with a sign change) locate extrema of the rate.
• Every numerical integral in imaging is a Riemann sum: the choice of rule sets the convergence rate, not just the answer.
• In R: D(), deriv(), integrate(), diff(), cumsum().

1.2 Series Approximation: Taylor Expansions and Finite Differences

Almost no equation in biomedical physics is solved exactly. What makes the subject tractable is that smooth functions can be approximated by polynomials, with a controllable error. That single idea produces the small-angle approximations of optics, the linearized Bloch equations of MRI, the finite-difference solvers of Chapter 7, and the local linearizations that gradient descent uses in Chapter 8.

1.2.1 Taylor and Maclaurin series

Theorem 1.1 (Taylor’s theorem with remainder). If \(f\) is \((n+1)\)-times differentiable on an interval containing \(a\) and \(x\), then

\[\begin{equation} f(x) = \sum_{i=0}^{n}\frac{f^{(i)}(a)}{i!}(x-a)^{i} + R_n(x), \qquad R_n(x) = \frac{f^{(n+1)}(\xi)}{(n+1)!}(x-a)^{n+1} \tag{15} \end{equation}\]

for some \(\xi\) between \(a\) and \(x\). When \(f\) is infinitely differentiable and \(R_n \to 0\), the series converges:

\[f(x) = \sum_{i=0}^{\infty}\frac{f^{(i)}(a)}{i!}(x-a)^{i}.\]

An expansion about \(a = 0\) is a Maclaurin series. The two most important cases in this book are

\[\begin{equation} e^{x} = \sum_{n=0}^{\infty}\frac{x^n}{n!} = 1 + x + \frac{x^2}{2!} + \frac{x^3}{3!} + \cdots, \tag{16} \end{equation}\]

\[\begin{equation} \sin x = \sum_{k=0}^{\infty}\frac{(-1)^k x^{2k+1}}{(2k+1)!} = x - \frac{x^3}{3!} + \frac{x^5}{5!} - \cdots. \tag{17} \end{equation}\]

Equation (16) is used directly in Section 1.6 to prove Euler’s identity, and the truncation \(e^{-x}\approx 1 - x\) is the linearization behind the “small optical density” limit of the Beer–Lambert law in Chapter 2.

The remainder term makes the accuracy explicit and quantitative: the error grows like \(\lvert x-a\rvert^{n+1}\), so a truncated series is accurate near the expansion point and degrades away from it, while adding terms improves accuracy everywhere within the radius of convergence.

taylor_sin <- function(x, K) {                     # partial sum through x^(2K+1)
  s <- 0
  for (k in 0:K) s <- s + (-1)^k * x^(2*k + 1)/factorial(2*k + 1)
  s
}

xg <- seq(-6, 6, length.out = 600)
tay <- rbind(
  data.frame(x = xg, y = sin(xg),          series = "sin(x)"),
  data.frame(x = xg, y = taylor_sin(xg,0), series = "order 1"),
  data.frame(x = xg, y = taylor_sin(xg,1), series = "order 3"),
  data.frame(x = xg, y = taylor_sin(xg,2), series = "order 5"),
  data.frame(x = xg, y = taylor_sin(xg,4), series = "order 9"))
ggplot(tay, aes(x, y, colour = series)) +
  geom_line(linewidth = 0.8) +
  scale_colour_manual(values = c("sin(x)" = "black", "order 1" = bpad_pal[1],
                                 "order 3" = bpad_pal[2], "order 5" = bpad_pal[3],
                                 "order 9" = bpad_pal[4])) +
  coord_cartesian(ylim = c(-2, 2)) +
  labs(title = "Maclaurin approximations of sin(x)", x = "x", y = "value")

Figure 6: Figure 1.6: Maclaurin approximations of sin(x). Each additional term extends the interval over which the polynomial tracks the sine, but all approximations are anchored at the expansion point x = 0.

## Absolute truncation error, tabulated
xt  <- c(0.5, 1, 2, 3)
err <- sapply(c(0, 1, 2, 4), function(K) abs(taylor_sin(xt, K) - sin(xt)))
dimnames(err) <- list(sprintf("x = %.1f", xt),
                      sprintf("through x^%d", 2*c(0,1,2,4) + 1))
knitr::kable(signif(err, 3),
             caption = "Absolute error of truncated Maclaurin series for sin(x).")

Table 2: Absolute error of truncated Maclaurin series for sin(x).

	through x^1	through x^3	through x^5	through x^9
x = 0.5	0.0206	0.000259	1.50e-06	0.00000
x = 1.0	0.1590	0.008140	1.96e-04	0.00000
x = 2.0	1.0900	0.243000	2.40e-02	0.00005
x = 3.0	2.8600	1.640000	3.84e-01	0.00419

1.2.2 Finite differences and the two sources of numerical error

Taylor’s theorem converts derivatives into computable differences. Expanding \(f(x \pm h)\) about \(x\),

\[f(x+h) = f(x) + h f'(x) + \tfrac{1}{2}h^2 f''(x) + \tfrac{1}{6}h^3 f'''(x) + \cdots\]

and rearranging gives the forward difference

\[\begin{equation} f'(x) \approx \frac{f(x+h)-f(x)}{h} + O(h), \tag{18} \end{equation}\]

while subtracting the expansion of \(f(x-h)\) cancels the second-order term and gives the central difference

\[\begin{equation} f'(x) \approx \frac{f(x+h)-f(x-h)}{2h} + O(h^2). \tag{19} \end{equation}\]

The central difference is one order more accurate, which is why it is the default stencil in the image-gradient operators of Chapter 7.

There is, however, a competing error. As \(h\) shrinks, the truncation error falls but the round-off error rises, because \(f(x+h)\) and \(f(x)\) become nearly equal and their difference loses significant digits. The total error is minimized at an intermediate \(h\), approximately \(\sqrt{\varepsilon_{\text{mach}}} \approx 10^{-8}\) for the forward difference and \(\varepsilon_{\text{mach}}^{1/3}\approx 10^{-5}\) for the central difference, where \(\varepsilon_{\text{mach}}\approx 2.2\times 10^{-16}\) is double-precision machine epsilon.

f  <- sin; fp <- cos; x0 <- 1
hs <- 10^seq(-1, -12, by = -0.25)
fd <- data.frame(
  h      = rep(hs, 2),
  err    = c(abs((f(x0 + hs) - f(x0))/hs           - fp(x0)),
             abs((f(x0 + hs) - f(x0 - hs))/(2*hs)  - fp(x0))),
  scheme = rep(c("forward, O(h)", "central, O(h^2)"), each = length(hs)))

best <- do.call(rbind, lapply(split(fd, fd$scheme),
                              function(d) d[which.min(d$err), ]))
knitr::kable(data.frame(scheme = best$scheme,
                        optimal_h = signif(best$h, 3),
                        min_error = signif(best$err, 3)),
             row.names = FALSE,
             caption = "Step size minimizing the total finite-difference error.")

Table 3: Step size minimizing the total finite-difference error.

scheme	optimal_h	min_error
central, O(h^2)	5.6e-06	0
forward, O(h)	0.0e+00	0

ggplot(fd, aes(h, err, colour = scheme)) +
  geom_line(linewidth = 0.9) + geom_point(size = 1.1) +
  scale_x_log10() + scale_y_log10() +
  scale_colour_manual(values = bpad_pal[1:2]) +
  labs(title = "Truncation error versus round-off error",
       x = "step size h (log scale)", y = "absolute error (log scale)")

Figure 7: Figure 1.7: Error of finite-difference derivative estimates of sin(x) at x = 1, versus step size h. Both curves show the characteristic V shape: truncation error dominates on the right, round-off error on the left. The central difference is both more accurate and optimal at a larger h.

Common misconception. “Smaller \(h\) always gives a better numerical derivative.” It does not. Below the optimum, catastrophic cancellation destroys accuracy faster than truncation error improves it. The same phenomenon limits how finely a physiological signal can be differentiated in time and how sharply an image gradient can be estimated in space, which is why edge detectors in Chapter 7 smooth before they differentiate.

Section 1.2 summary.

• Taylor’s theorem approximates a smooth function by a polynomial with an explicit, bounded remainder.
• \(e^x\) and \(\sin x\) have Maclaurin series that recur throughout the book (Euler’s identity, small-angle optics, linearized attenuation).
• Finite differences follow directly from Taylor expansions: forward is \(O(h)\), central is \(O(h^2)\).
• Numerical differentiation balances truncation error against round-off; an intermediate step size is optimal.

Checkpoint 1.3. Adding terms to a Taylor series improves accuracy near the expansion point immediately, but far away only slowly. Explain this using the remainder term \(R_n\), and predict what happens to the order-9 sine approximation at \(x = 10\).

1.3 Scalars, Vectors, Matrices, and Tensors

Biomedical data are almost never scalars. A pixel is one number, but an image is a matrix, a volume is a third-order array, and a dynamic volume is a fourth-order array. This section fixes the vocabulary and the operations; Section 1.17 develops the deeper structure (eigenvalues, SVD) that these objects support.

1.3.1 The dimensional hierarchy

A scalar is a single number, a zero-dimensional object: \(1,\ 2.7,\ -30,\ \pi,\ e\).

A vector is a one-dimensional ordered array. In this textbook vectors are written as columns,

\[\mathbf{a} = \begin{bmatrix} 10 \\ 5\end{bmatrix}, \qquad \mathbf{b} = \begin{bmatrix} e \\ 25 \\ 77 \end{bmatrix},\]

with dimensions \((n\times 1)\); the row vector \(\mathbf{a}^{T}\) is \((1\times n)\). The terms list, array, and vector are used loosely across programming languages, but in biophysics “vector” carries the specific meaning of an element of a vector space, an object with magnitude and direction that transforms in a prescribed way under a change of basis.

A matrix is a two-dimensional \((n \times m)\) array, with \(n\) rows and \(m\) columns:

\[A = \begin{bmatrix} 10 & 4 & -1 \\ 5 & 0 & 3 \\ 6 & -200 & 3 \end{bmatrix}, \qquad B = \begin{bmatrix} 3 & 7 \\ 25 & -9 \\ 77 & -20 \end{bmatrix}.\]

A tensor is the general multi-dimensional object; its order (or rank) is the number of indices required to address an element. A scalar is order 0, a vector order 1, a matrix order 2.

A useful distinction, and a common conflation. In physics, a second-order tensor is defined by how its components transform under a change of coordinates, not merely by being a rectangular array of numbers. The Cauchy stress tensor

\[\sigma = \begin{bmatrix} \sigma_{xx} & \sigma_{xy} & \sigma_{xz} \\ \sigma_{yx} & \sigma_{yy} & \sigma_{yz} \\ \sigma_{zx} & \sigma_{zy} & \sigma_{zz} \end{bmatrix}\]

is a genuine tensor: it describes a physical state of a deformed material and its components change predictably when the axes rotate. The diffusion tensor of Chapter 4 is another. By contrast, a table of intensity values indexed by pixel row and column is a matrix that is not generally a tensor in this stricter sense, because rotating the patient does not transform pixel intensities the way tensor components transform. In machine learning (Chapter 8) the word “tensor” is used loosely to mean “multi-dimensional array”; both usages appear in this book, and the context disambiguates them.

1.3.2 Core operations

Addition and subtraction are defined element-wise and require identical dimensions. Multiplication is subtler and depends on the objects involved.

• Inner (dot) product. For \(\mathbf{a}, \mathbf{b}\in\mathbb{R}^n\),

\[\begin{equation} \mathbf{a}\cdot\mathbf{b} = \mathbf{a}^{T}\mathbf{b} = \sum_{i=1}^{n} a_i b_i = \lVert\mathbf{a}\rVert\,\lVert\mathbf{b}\rVert\cos\theta. \tag{20} \end{equation}\]

The dot product measures projection, and through the cosine it measures similarity. Correlation coefficients, matched filters, Fourier coefficients, and the kernels of Chapter 8 are all inner products. * Cross product. For \(\mathbf{a},\mathbf{b}\in\mathbb{R}^3\), \(\mathbf{a}\times\mathbf{b}\) is the vector orthogonal to both, with magnitude \(\lVert\mathbf{a}\rVert\lVert\mathbf{b}\rVert\sin\theta\). It supplies torque, angular momentum (hence the precession of Chapter 4), and surface normals (hence the flux integrals of Section 1.12). * Matrix product. \(AB\) is defined when the inner dimensions agree: \((m\times n)(n\times k) = (m\times k)\), with \((AB)_{ij} = \sum_{l} A_{il}B_{lj}\), computed row-by-column. Matrix multiplication is associative and distributive but not commutative. * Transpose. \((A^{T})_{ij} = A_{ji}\). A matrix is symmetric when \(A = A^{T}\). * Trace. \(\mathrm{tr}(A) = \sum_i A_{ii}\), defined for square matrices. It equals the sum of the eigenvalues (Section 1.17), and in DTI it equals three times the mean diffusivity. * Determinant. \(\det(A)\) is the factor by which the linear map \(A\) scales volume; \(\det(A) = 0\) signals a singular, non-invertible matrix, that is, a transformation that collapses volume to zero. * Inverse. \(A^{-1}\) satisfies \(AA^{-1} = A^{-1}A = I\) and exists precisely when \(\det(A)\ne 0\).

x <- 3; y <- -5                       # scalars
v <- c(x, y); w <- c(2*x, -14*y)      # vectors in R^2
m <- matrix(1:12,  nrow = 3, byrow = TRUE)   # 3 x 4
n <- matrix(24:35, nrow = 4, byrow = TRUE)   # 4 x 3

m %*% n                               # (3x4)(4x3) -> 3x3

##      [,1] [,2] [,3]
## [1,]  300  310  320
## [2,]  756  782  808
## [3,] 1212 1254 1296

dim(m %*% n)

## [1] 3 3

M <- cbind(v, w); M

##       v  w
## [1,]  3  6
## [2,] -5 70

solve(M)                              # inverse

##         [,1]    [,2]
## v 0.29166667 -0.0250
## w 0.02083333  0.0125

round(M %*% solve(M), 12)             # identity, to numerical precision

##      [,1] [,2]
## [1,]    1    0
## [2,]    0    1

c(dot_product = sum(v*w),
  determinant = det(M),
  trace       = sum(diag(M)),
  angle_deg   = acos(sum(v*w)/(sqrt(sum(v^2))*sqrt(sum(w^2))))*180/pi)

## dot_product determinant       trace   angle_deg 
##   -332.0000    240.0000     73.0000    144.1372

Two R conventions are worth memorizing now, because they are a common source of silent errors: %*% is matrix multiplication, whereas the bare * is element-wise (Hadamard) multiplication; and t() transposes.

1.3.3 Biomedical images as tensors

The dimensional hierarchy maps directly onto imaging (developed fully in Chapter 7):

Object	Order	Imaging instance
Scalar	0	one voxel intensity; a Hounsfield unit; an SUV
Vector	1	a time course at one voxel; an RGB triple; a \(b\)-vector in diffusion MRI
Matrix	2	one image slice; a design matrix of patients \(\times\) features
3rd-order array	3	a CT or MR volume \(I \in \mathbb{R}^{n_x\times n_y\times n_z}\); a diffusion tensor field
4th-order array	4	a dynamic (4-D) acquisition: fMRI, DCE-MRI, or a gated cardiac series

Slicing a higher-order array along any axis produces a lower-order one. Fixing \((x,y,z)\) in a 4-D fMRI dataset and letting \(t\) vary yields a vector, the voxel time course; fixing \((z,t)\) yields a matrix, one cross-sectional image; fixing \(t\) alone yields a volume.

The demonstration below constructs a synthetic 4-D dataset with a known activation, so the slicing operations are fully reproducible without a download. In the 5-D spacekime framework, the spacekime analytic representation generalizes 4-D spacetime modelling to 2-D complex time (kime), where longitudinal dynamics carry both a kime magnitude and a kime phase.

set.seed(2024)
nx <- 32; ny <- 32; nz <- 12; nt <- 40

g <- expand.grid(i = 1:nx, j = 1:ny)
r <- sqrt((g$i - nx/2)^2 + (g$j - ny/2)^2)
brain <- matrix(exp(-r^2/(2*8^2)), nx, ny)            # smooth "head" profile

hrf  <- function(tt) pmax(0, (tt/6)^2 * exp(-(tt - 6)/3))   # canonical HRF shape
stim <- rep(c(rep(1, 5), rep(0, 5)), length.out = nt)       # block design
resp <- stats::filter(stim, hrf(0:9), sides = 1)            # convolve (Section 1.7)
resp[is.na(resp)] <- 0

vol4d <- array(0, dim = c(nx, ny, nz, nt))
for (kz in 1:nz) for (tt in 1:nt)
  vol4d[, , kz, tt] <- brain*(1 + 0.05*resp[tt]*(kz == 6)) +
                       rnorm(nx*ny, 0, 0.02)

dim(vol4d)                          # 4th-order array: x, y, z, t

## [1] 32 32 12 40

length(vol4d[16, 16, 6, ])          # vector : one voxel time course

## [1] 40

dim(vol4d[, , 6, 20])               # matrix : one slice at one time point

## [1] 32 32

dim(vol4d[, , , 20])                # 3rd-order: one whole volume

## [1] 32 32 12

op <- par(mfrow = c(1, 2), mar = c(4, 4, 2.4, 1))
image(t(vol4d[, , 6, 20][nx:1, ]), col = gray.colors(256, 0, 1),
      axes = FALSE, asp = 1, main = "Slice z = 6, t = 20  (matrix)")
box(col = "grey70")
plot(1:nt, vol4d[16, 16, 6, ], type = "b", pch = 16, cex = 0.6,
     col = bpad_pal[1], xlab = "time frame", ylab = "intensity",
     main = "Voxel (16,16,6) time course  (vector)")
lines(1:nt, brain[16,16]*(1 + 0.05*resp), col = bpad_pal[2], lwd = 2)
legend("topright", c("measured", "noise-free model"), bty = "n", cex = 0.8,
       col = bpad_pal[1:2], lty = 1, pch = c(16, NA))

Figure 8: Figure 1.8: Slices of a synthetic 4-D dataset. Left: one axial slice at a fixed time (a matrix). Right: the time course of a single voxel (a vector), with the stimulus-locked haemodynamic response visible above the noise.

par(op)

The corresponding demonstration on real data is available and uses the same slicing logic. It is guarded by run_network because it requires a download; set the flag to TRUE at the top of the source to execute it.

## Real 4-D fMRI volume from the SOCR archive (requires network access and the
## brainR / oro.nifti packages). The slicing syntax is identical to the
## synthetic demonstration above.
library(brainR); library(httr)
httr::set_config(httr::config(ssl_verifypeer = 0L))

url  <- "https://socr.umich.edu/HTML5/BrainViewer/data/fMRI_FilteredData_4D.nii.gz"
dest <- file.path(tempdir(), "fMRI_FilteredData_4D.nii.gz")
invisible(httr::GET(url, httr::write_disk(dest, overwrite = TRUE)))
vol <- oro.nifti::readNIfTI(dest, reorient = FALSE)

d <- dim(vol); d
mid <- floor(d/2)
length(vol[mid[1], mid[2], mid[3], ])     # voxel time course (vector)
dim(vol[, , mid[3], mid[4]])              # one slice (matrix)
image(t(vol[, , mid[3], mid[4]]), col = gray.colors(256, 0, 1), asp = 1, axes = FALSE)

Section 1.3 summary.

• Scalar \(\rightarrow\) vector \(\rightarrow\) matrix \(\rightarrow\) higher-order array is a hierarchy of index counts, and slicing moves down it.
• The dot product measures projection and similarity; the determinant measures volume scaling; the trace sums the diagonal and the eigenvalues.
• Matrix multiplication is associative but not commutative; in R, %*% and * are different operations.
• A 4-D acquisition (fMRI, DCE-MRI, gated CT) is a fourth-order array whose slices are the images and time courses actually analyzed.

Checkpoint 1.4. A dynamic contrast-enhanced MRI study acquires a \(256\times 256\times 40\) volume every 5 s for 6 min. What is the order and shape of the full dataset? Which slice would you extract to plot a lesion’s enhancement curve, and which to display an image?

1.4 Kinematics: Displacement, Velocity, and Acceleration

The three fundamental kinematic quantities are related by one derivative and one integral each, which makes them the cleanest possible illustration of Section 1.1. They are also the basis for motion correction in imaging (Chapters 4 and 7), gait and biomechanics analysis, and the ballistic transport terms in Chapter 3.

Displacement \(\Delta \mathbf{x}\) is a change in position. Velocity \(\mathbf{v}\) is its time rate of change, and acceleration \(\mathbf{a}\) is the time rate of change of velocity. All three are vectors, carrying both magnitude and direction; the scalar magnitudes have separate names, distance for displacement and speed for velocity, while acceleration has no separate scalar name. The chain is

\[\begin{equation} \mathbf{v}(t) = \frac{d\mathbf{x}}{dt}, \qquad \mathbf{a}(t) = \frac{d\mathbf{v}}{dt} = \frac{d^{2}\mathbf{x}}{dt^{2}}, \qquad \mathbf{x}(t) = \mathbf{x}_0 + \int_{0}^{t}\mathbf{v}(\tau)\,d\tau. \tag{21} \end{equation}\]

Projectile motion. With constant gravitational acceleration \(g = 9.8\) m s\(^{-2}\) acting downward and no drag, integrating twice gives the vertical motion

\[\begin{equation} y(t) = y_0 + v_0 t - \tfrac{1}{2}g t^{2}, \qquad v(t) = v_0 - g t, \qquad a(t) = -g. \tag{22} \end{equation}\]

Because the horizontal component of acceleration is zero, the horizontal velocity is constant, and the trajectory is a parabola.

y0 <- 2; v0 <- 50; g <- 9.8
disp <- function(t) y0 + v0*t - 0.5*g*t^2
vel  <- function(t) v0 - g*t
acc  <- function(t) rep(-g, length(t))

t_apex <- v0/g
tg <- seq(0, 2*t_apex, length.out = 400)
km <- rbind(data.frame(t = tg, y = disp(tg), q = "displacement y(t)  [m]"),
            data.frame(t = tg, y = vel(tg),  q = "velocity v(t)  [m/s]"),
            data.frame(t = tg, y = acc(tg),  q = "acceleration a(t)  [m/s^2]"))

ggplot(km, aes(t, y, colour = q)) +
  geom_hline(yintercept = 0, colour = "grey70") +
  geom_vline(xintercept = t_apex, linetype = "dashed", colour = "grey40") +
  geom_line(linewidth = 0.9) +
  annotate("text", x = t_apex + 0.15, y = disp(t_apex), hjust = 0, size = 3.2,
           label = sprintf("apex: t = v0/g = %.2f s,  y = %.1f m",
                           t_apex, disp(t_apex))) +
  scale_colour_manual(values = bpad_pal[1:3]) +
  labs(title = "Projectile motion in the vertical direction",
       x = "time (s)", y = "value")

Figure 9: Figure 1.9: Vertical projectile motion. The quadratic displacement, linear velocity, and constant acceleration are successive derivatives of one another. The apex occurs where the velocity crosses zero.

Section 1.4 summary.

• Velocity is the first derivative of displacement and the first integral of acceleration; displacement is the second integral of acceleration.
• All three are vectors; the parabolic trajectory of a projectile is a direct consequence of integrating a constant acceleration twice.
• The apex of the trajectory is the critical point of \(y(t)\), where \(v(t) = 0\), an instance of Section 1.1.7.

Checkpoint 1.5. If gravity is the only force acting, why does the horizontal velocity of a projectile remain constant? Which of Eqs. (21) makes this immediate?

1.5 Polynomials, Exponentials, and Logarithms

Three families of elementary functions describe nearly every scalar relationship in this book, and the single most useful thing to know about them is how fast they grow.

Polynomials are finite sums of non-negative integer powers,

\[a(x) = 5x + 5x^2 + x^5, \qquad b(x,y) = x^2 + xy^3 - 9y^2 + 3.\]

Exponential functions carry the variable in the exponent; the sign of the exponent determines unbounded growth or asymptotic decay,

\[c(x) = 9e^{2x}, \qquad d(x) = 9e^{-8x+2}.\]

Logarithms are the inverses of exponentials, indexed by their base, with \(\log \equiv \log_{10}\) and \(\ln \equiv \log_e\). The change-of-base identity is

\[\begin{equation} \log_b(x) = \frac{\log_a(x)}{\log_a(b)}. \tag{23} \end{equation}\]

Their asymptotic ordering is the essential fact:

\[\begin{equation} \ln x \;\ll\; x \;\ll\; x^2 \;\ll\; x^n \;\ll\; e^{x} \qquad \text{as } x \to \infty. \tag{24} \end{equation}\]

xg <- seq(0.1, 8, length.out = 500)
gr <- rbind(data.frame(x = xg, y = log(xg), fn = "ln(x)"),
            data.frame(x = xg, y = xg,      fn = "x"),
            data.frame(x = xg, y = xg^2,    fn = "x^2"),
            data.frame(x = xg, y = xg^3,    fn = "x^3"),
            data.frame(x = xg, y = exp(xg), fn = "e^x"))
gr <- gr[gr$y > 0, ]

ggplot(gr, aes(x, y, colour = fn)) +
  geom_line(linewidth = 0.9) +
  scale_y_log10() +
  scale_colour_manual(values = bpad_pal[1:5]) +
  labs(title = "Asymptotic growth of elementary functions",
       subtitle = "logarithmic vertical axis",
       x = "x", y = "f(x)  (log scale)")

Figure 10: Figure 1.10: Growth rates of logarithmic, linear, quadratic, cubic, and exponential functions on a logarithmic vertical axis. On this scale an exponential is a straight line and a power law is a curve of slope equal to its exponent; the exponential overtakes every polynomial.

The relevance is not abstract. A single 24-hour bacterial culture traverses all three regimes: an initial lag phase that is roughly polynomial, a log phase that is exponential, and a stationary phase governed by the logistic saturation of Section 1.1.7. Recognizing which regime a measurement samples determines which model can legitimately be fitted to it.

Logarithms also serve a second, purely practical role: they linearize exponential relationships. Taking logarithms of \(N(t) = N_0 e^{kt}\) gives \(\ln N = \ln N_0 + kt\), converting a curve-fitting problem into a straight-line regression. This is precisely how the attenuation coefficient is extracted in Chapter 5, how the decay constant is extracted in Chapter 6, and how \(T_1\) and \(T_2\) are estimated in Chapter 4. Section 1.16 examines when this linearization is statistically appropriate and when it is not.

From equation to biomedical data. Pharmacokinetic clearance is exponential decay (Section 1.13). Receptor occupancy follows the Hill equation, a rational function of ligand concentration. Photon attenuation through tissue is exponential in depth (Chapter 5). Sound attenuation is exponential in depth and frequency-dependent (Chapter 3). In each case the same functional family appears with a different physical name for its rate constant, which is why identifying the family is often more informative than memorizing the individual law.

Section 1.5 summary.

• Polynomials, exponentials, and logarithms are ordered by asymptotic growth: \(\ln x \ll x^n \ll e^x\).
• Exponentials solve first-order rate laws; logarithms invert them and linearize them for fitting.
• Distinguishing polynomial, exponential, and saturating (logistic) regimes determines which model may legitimately be fitted to a dataset.

Checkpoint 1.6. Rank \(\log x\), \(x\), \(x^2\), \(e^x\) from slowest to fastest growth at large \(x\). Then explain why a semi-logarithmic plot of an exponential decay is a straight line, and what its slope represents physically.

1.6 Complex Numbers and Phasor Notation

Complex numbers are not an algebraic curiosity in biomedical physics; they are the natural language of anything that oscillates. Every alternating current, every propagating wave, every Fourier coefficient, and every MR signal is complex-valued, because two real numbers (amplitude and phase) are needed where one real number will not do.

1.6.1 Definition and arithmetic

Definition 1.3. A complex number is \(z = x + iy\) with \(x = \mathrm{Re}(z)\in\mathbb{R}\), \(y = \mathrm{Im}(z)\in\mathbb{R}\), and \(i^2 = -1\). The complex numbers \(\mathbb{C}\) form a field under the operations

\[\begin{equation} z_1 + z_2 = (x_1+x_2) + i(y_1+y_2), \qquad z_1 z_2 = (x_1x_2 - y_1y_2) + i(x_1y_2 + x_2y_1). \tag{25} \end{equation}\]

The conjugate is \(\bar{z} = z^{*} = x - iy\); the modulus is \(\lvert z\rvert = \sqrt{z\bar{z}} = \sqrt{x^2+y^2}\); the argument is \(\theta = \arg(z) = \mathrm{atan2}(y,x)\). The conjugate rationalizes denominators,

\[\frac{1}{z} = \frac{\bar{z}}{\lvert z\rvert^{2}} = \frac{x - iy}{x^2+y^2},\]

and, more importantly, it converts a complex amplitude into a real power: \(\lvert z\rvert^{2} = z\bar{z}\). That identity is why a magnitude image can be formed from a complex MR signal (Chapter 4) and why a power spectrum is real even when the spectrum is not (Section 1.9).

1.6.2 Polar form and Euler’s identity

Substituting \(x = i\theta\) into the exponential series ((16)) and separating even and odd powers reproduces the cosine and sine series, giving Euler’s formula

\[\begin{equation} e^{i\theta} = \cos\theta + i\sin\theta, \tag{26} \end{equation}\]

and therefore the polar representation

\[\begin{equation} z = \lvert z\rvert e^{i\theta} = \lvert z\rvert(\cos\theta + i\sin\theta). \tag{27} \end{equation}\]

The special case \(\theta = \pi\) yields Euler’s identity, \(e^{i\pi} + 1 = 0\).

Polar form turns multiplication into addition of angles, which is the entire computational advantage:

\[\begin{equation} z_1 z_2 = \lvert z_1\rvert\lvert z_2\rvert e^{i(\theta_1+\theta_2)}, \qquad \frac{z_1}{z_2} = \frac{\lvert z_1\rvert}{\lvert z_2\rvert}e^{i(\theta_1-\theta_2)}, \qquad z^{n} = \lvert z\rvert^{n}e^{in\theta}, \tag{28} \end{equation}\]

the last being De Moivre’s theorem.

1.6.3 Phasors

Definition 1.4 (Phasor). A time-harmonic real signal \(A\cos(\omega t + \phi)\) is represented by the complex constant

\[\begin{equation} \tilde{A} = A e^{i\phi} = A\angle\phi, \tag{29} \end{equation}\]

with the original signal recovered as \(A\cos(\omega t + \phi) = \mathrm{Re}\{\tilde{A}e^{i\omega t}\}\). The phasor strips away the common \(e^{i\omega t}\) factor shared by every signal at the same frequency and retains only what distinguishes them: amplitude and phase.

The payoff is that calculus becomes algebra. Because \(\frac{d}{dt}e^{i\omega t} = i\omega\, e^{i\omega t}\), differentiation of a phasor is multiplication by \(i\omega\) and integration is division by \(i\omega\). A linear differential equation in \(t\) becomes a linear algebraic equation in \(\tilde{A}\).

Phasors appear in four distinct places in this book:

1. Electrical impedance. Tissue impedance is frequency-dependent and complex, \(Z(\omega) = R + iX\), with \(R\) the resistance and \(X\) the reactance (Section 1.6.4).
2. Bioimpedance spectroscopy. The phase angle of \(Z\) tracks cell-membrane integrity and the balance between intra- and extracellular fluid, and is used clinically in body-composition and tissue-viability assessment.
3. Wave optics and interferometry. Optical fields in OCT and interferometric microscopy (Chapter 2) are phasors whose phase encodes optical path length.
4. Fourier analysis and MRI. The Fourier transform (Section 1.9) produces a phasor at every frequency; the MR receiver measures a quadrature (complex) signal directly, and the reconstructed image inherits both magnitude and phase (Chapter 4).

Worked Example 1.3 (phasor addition). Two voltages at the same frequency,

\[V_1(t) = 3\cos(\omega t + \pi/4), \qquad V_2(t) = 4\cos(\omega t - \pi/6),\]

are summed. Adding the time-domain expressions directly requires trigonometric identities; adding the phasors requires only complex arithmetic:

\[\tilde{V}_1 = 3e^{i\pi/4} = 3\left(\tfrac{\sqrt{2}}{2} + i\tfrac{\sqrt{2}}{2}\right) \approx 2.121 + 2.121\,i,\] \[\tilde{V}_2 = 4e^{-i\pi/6} = 4\left(\tfrac{\sqrt{3}}{2} - i\tfrac{1}{2}\right) \approx 3.464 - 2.000\,i,\] \[\tilde{V} = \tilde{V}_1 + \tilde{V}_2 \approx 5.585 + 0.121\,i.\]

Hence \(\lvert\tilde{V}\rvert \approx 5.587\) and \(\arg\tilde{V}\approx 0.0217\) rad, so

\[V(t) \approx 5.587\cos(\omega t + 0.0217).\]

Note that the sum’s amplitude, \(5.59\), is far below the arithmetic sum \(3 + 4 = 7\): the two components are nearly \(75^\circ\) out of phase and partially cancel. Phase, not just amplitude, determines what a detector measures.

V1 <- 3*exp(1i*pi/4)
V2 <- 4*exp(-1i*pi/6)
Vs <- V1 + V2
c(Mod_V1 = Mod(V1), Mod_V2 = Mod(V2),
  Mod_sum = Mod(Vs), Arg_sum_rad = Arg(Vs))

##      Mod_V1      Mod_V2     Mod_sum Arg_sum_rad 
##  3.00000000  4.00000000  5.58673940  0.02171748

seg <- data.frame(
  x0 = c(0, Re(V1), 0), y0 = c(0, Im(V1), 0),
  x1 = c(Re(V1), Re(Vs), Re(Vs)), y1 = c(Im(V1), Im(Vs), Im(Vs)),
  lab = c("V1", "V2 (head to tail)", "V1 + V2"))

p1 <- ggplot(seg) +
  geom_hline(yintercept = 0, colour = "grey75") +
  geom_vline(xintercept = 0, colour = "grey75") +
  geom_segment(aes(x = x0, y = y0, xend = x1, yend = y1, colour = lab),
               arrow = grid::arrow(length = grid::unit(0.18, "cm")), linewidth = 0.9) +
  scale_colour_manual(values = bpad_pal[1:3]) +
  coord_fixed() +
  labs(title = "Phasor addition", x = "Re", y = "Im")

tg <- seq(0, 2*pi, length.out = 500)
p2 <- ggplot(data.frame(t = tg), aes(t)) +
  geom_line(aes(y = Re(V1*exp(1i*t)), colour = "V1"),      linewidth = 0.8) +
  geom_line(aes(y = Re(V2*exp(1i*t)), colour = "V2"),      linewidth = 0.8) +
  geom_line(aes(y = Re(Vs*exp(1i*t)), colour = "V1 + V2"), linewidth = 1.0) +
  scale_colour_manual(values = bpad_pal[c(1,2,3)]) +
  labs(title = "Corresponding time-domain signals",
       x = "omega*t (rad)", y = "amplitude")

bpad_grid(p1, p2, ncol = 2)

Figure 11: Figure 1.11: Phasor addition in the complex plane (left) and the corresponding time-domain signals (right). The resultant amplitude is well below the arithmetic sum of the components because of the phase difference.

1.6.4 Impedance and bioimpedance spectroscopy

Impedance relates the voltage phasor to the current phasor, \(\tilde{V} = Z\tilde{I}\). Modelling a tissue as a resistor \(R\) (extracellular path) in parallel with a capacitor \(C\) (cell membranes) gives

\[\begin{equation} Z(\omega) = \frac{R}{1 + i\omega RC} = \frac{R\left(1 - i\omega RC\right)}{1 + (\omega RC)^2}, \tag{30} \end{equation}\]

\[\begin{equation} \lvert Z(\omega)\rvert = \frac{R}{\sqrt{1 + (\omega RC)^2}}, \qquad \phi(\omega) = -\arctan(\omega RC). \tag{31} \end{equation}\]

The physics is transparent in the two limits. At low frequency, \(\omega RC \ll 1\), the capacitive branch is effectively open and current flows only through the extracellular fluid, so \(Z \to R\) and \(\phi\to 0\). At high frequency, \(\omega RC \gg 1\), the membranes are short-circuited by capacitive coupling and current penetrates the cells, so \(\lvert Z\rvert \to 0\) and \(\phi \to -90^\circ\). The transition occurs at the characteristic frequency \(\omega_c = 1/RC\).

This frequency dependence is what makes bioimpedance spectroscopy clinically informative: sweeping \(\omega\) separates the intra- and extracellular water compartments, supporting assessment of hydration status, body composition, oedema, and tissue viability.

R_t <- 500; C_t <- 1e-6
omega <- 10^seq(0, 6, length.out = 600)
Z <- R_t/(1 + 1i*omega*R_t*C_t)
omega_c <- 1/(R_t*C_t)
cat(sprintf("characteristic frequency: %.0f rad/s  (%.1f Hz)\n",
            omega_c, omega_c/(2*pi)))

## characteristic frequency: 2000 rad/s  (318.3 Hz)

zd <- data.frame(omega = omega, mag = Mod(Z), phase = Arg(Z)*180/pi,
                 re = Re(Z), im = -Im(Z))

pa <- ggplot(zd, aes(omega, mag)) + geom_line(colour = bpad_pal[1], linewidth = 0.9) +
  geom_vline(xintercept = omega_c, linetype = "dashed", colour = "grey45") +
  scale_x_log10() + labs(title = "|Z| (Bode magnitude)",
                         x = "omega (rad/s)", y = "|Z| (ohm)")
pb <- ggplot(zd, aes(omega, phase)) + geom_line(colour = bpad_pal[2], linewidth = 0.9) +
  geom_vline(xintercept = omega_c, linetype = "dashed", colour = "grey45") +
  scale_x_log10() + labs(title = "phase of Z (Bode phase)",
                         x = "omega (rad/s)", y = "phase (degrees)")
pc <- ggplot(zd, aes(re, im)) + geom_path(colour = bpad_pal[3], linewidth = 0.9) +
  coord_fixed() + labs(title = "Cole-Cole plot",
                       x = "Re(Z) (ohm)", y = "-Im(Z) (ohm)")
bpad_grid(pa, pb, pc, ncol = 3)

Figure 12: Figure 1.12: Bioimpedance of a parallel RC tissue model (R = 500 ohm, C = 1 microfarad). Left and centre: Bode magnitude and phase. Right: the Cole-Cole plot, a semicircle in the complex impedance plane whose apex marks the characteristic frequency.

Section 1.6 summary.

• A complex number carries two real quantities, amplitude and phase, which is exactly what an oscillation requires.
• Euler’s formula \(e^{i\theta} = \cos\theta + i\sin\theta\) links exponentials to trigonometry; polar form turns multiplication into angle addition.
• A phasor represents a sinusoid by its amplitude and phase, converting differential equations into algebraic ones.
• Bioimpedance \(Z(\omega) = R/(1+i\omega RC)\) separates extracellular from intracellular conduction as frequency increases; the Cole-Cole plot displays this as a semicircle.

Checkpoint 1.7. Two phasors have equal magnitude and opposite phase. What is their sum, and what does this imply for a measurement that records only \(\lvert Z\rvert\) and discards phase?

1.7 Linear Systems, Convolution, and the Convolution Theorem

Almost every measurement instrument in this textbook blurs what it measures. A CT detector integrates over a finite aperture, an ultrasound transducer has a finite beam width, a microscope has a diffraction-limited point spread, and an amplifier has a finite bandwidth. The mathematics that describes all of these at once is the theory of linear systems, and its central operation is convolution.

Definition 1.5 (Linear system). An operator \(\mathcal{L}\) is linear if for all inputs \(x_1, x_2\) and scalars \(a\),

\[\begin{equation} \mathcal{L}\{x_1 + x_2\} = \mathcal{L}\{x_1\} + \mathcal{L}\{x_2\} \quad\text{(additivity)}, \qquad \mathcal{L}\{a\,x\} = a\,\mathcal{L}\{x\} \quad\text{(homogeneity)}. \tag{32} \end{equation}\]

It is time-invariant (LTI) if delaying the input merely delays the output. Most biomedical instruments are approximately LTI over their designed operating range; departures from linearity (detector saturation, tissue nonlinearity in harmonic imaging) are treated as corrections to this baseline.

1.7.1 The impulse response and convolution

The impulse response \(h(t)\) of an LTI system is its output when the input is a Dirac delta \(\delta(t)\). The delta is defined by its sifting property,

\[\int_{-\infty}^{\infty} f(t)\,\delta(t-t_0)\,dt = f(t_0),\]

which says that any input can be written as a superposition of shifted impulses. Linearity and time invariance then force the output to be the convolution

\[\begin{equation} y(t) = (x * h)(t) = \int_{-\infty}^{\infty} x(\tau)\,h(t-\tau)\,d\tau. \tag{33} \end{equation}\]

Read this as weighted superposition: the output at time \(t\) is a weighted sum of all input values, each weighted by the system’s response to an impulse delivered at that time. Convolution is commutative, associative, distributive, and has \(\delta\) as its identity, \(x * \delta = x\).

1.7.2 Discrete convolution

For sampled data the integral becomes a sum:

\[\begin{equation} y[n] = (x * h)[n] = \sum_{k=-\infty}^{\infty} x[k]\,h[n-k]. \tag{34} \end{equation}\]

Convolving finite sequences of lengths \(N\) and \(M\) produces a sequence of length \(N + M - 1\) (linear convolution); wrapping the result to length \(N\) produces circular convolution, which is what the FFT computes.

x <- c(rep(0, 20), rep(1, 20), rep(0, 20))
h <- dnorm(seq(-3, 3, length.out = 15)); h <- h/sum(h)   # normalized kernel
y <- convolve(x, rev(h), type = "open")                  # linear convolution

c(length_x = length(x), length_h = length(h), length_y = length(y),
  area_in = sum(x), area_out = sum(y))

## length_x length_h length_y  area_in area_out 
##       60       15       74       20       20

cv <- rbind(data.frame(n = seq_along(x), v = x,     s = "input x[n]"),
            data.frame(n = seq_along(h), v = h*10,  s = "kernel h[n] (x10)"),
            data.frame(n = seq_along(y), v = y,     s = "output y = x * h"))
ggplot(cv, aes(n, v, colour = s)) +
  geom_line(linewidth = 0.9) +
  scale_colour_manual(values = bpad_pal[1:3]) +
  labs(title = "Discrete convolution: Gaussian smoothing of a rectangular pulse",
       x = "sample index n", y = "amplitude")

Figure 13: Figure 1.13: One-dimensional convolution. A rectangular input convolved with a normalized Gaussian kernel produces a smoothed output with softened edges. Total area is preserved because the kernel sums to one.

1.7.3 Two-dimensional convolution and image filtering

For an image \(I(x,y)\) and kernel \(K(u,v)\),

\[\begin{equation} (I * K)(x,y) = \sum_{u}\sum_{v} I(x-u,\,y-v)\,K(u,v). \tag{35} \end{equation}\]

Three kernels recur throughout medical image processing (Chapter 7) and reappear as the learned filters of convolutional neural networks (Chapter 8):

• Gaussian, \(K \propto e^{-(u^2+v^2)/2\sigma^2}\): suppresses noise at the cost of resolution. It is the standard preprocessing step in essentially every imaging pipeline.
• Laplacian, \(\nabla^2\): enhances edges and fine detail by responding to curvature (Section 1.12).
• Sobel / Prewitt: directional first-derivative operators used for boundary delineation.

conv2 <- function(img, ker) {
  ni <- nrow(img); nj <- ncol(img); nk <- nrow(ker); nl <- ncol(ker)
  ph <- nk %/% 2; pw <- nl %/% 2
  pad <- matrix(0, ni + 2*ph, nj + 2*pw)
  pad[(ph+1):(ph+ni), (pw+1):(pw+nj)] <- img
  kf  <- ker[nk:1, nl:1]                 # flip the kernel: convolution, not correlation
  out <- matrix(0, ni, nj)
  for (i in 1:ni) for (j in 1:nj)
    out[i, j] <- sum(pad[i:(i+nk-1), j:(j+nl-1)] * kf)
  out
}

set.seed(42)
N   <- 64
rad <- outer(1:N, 1:N, function(i, j) sqrt((i - 32)^2 + (j - 32)^2))
img <- 0.3 + 0.5*exp(-rad^2/50)*(rad < 12) + matrix(rnorm(N*N, 0, 0.05), N)

gk  <- outer(seq(-2, 2, length.out = 5), seq(-2, 2, length.out = 5),
             function(a, b) exp(-(a^2 + b^2)/2)); gk <- gk/sum(gk)
lap <- matrix(c(0, 1, 0, 1, -4, 1, 0, 1, 0), 3, byrow = TRUE)

smoothed <- conv2(img, gk)
edges    <- conv2(img, lap)

bg <- 45:60                      # interior background: away from lesion and border
cat(sprintf("background noise SD: original %.4f -> smoothed %.4f (%.0f%% reduction)\n",
            sd(img[bg, bg]), sd(smoothed[bg, bg]),
            100*(1 - sd(smoothed[bg, bg])/sd(img[bg, bg]))))

## background noise SD: original 0.0532 -> smoothed 0.0153 (71% reduction)

op <- par(mfrow = c(1, 3), mar = c(0.6, 0.6, 2.2, 0.6))
show_img(img,       "noisy original")
show_img(smoothed,  "Gaussian smoothed")
show_img(abs(edges), "Laplacian edge response")

Figure 14: Figure 1.14: Two-dimensional convolution on a synthetic lesion image. Left: noisy original. Centre: Gaussian smoothing suppresses noise and blurs the lesion boundary. Right: the Laplacian responds to the boundary itself.

par(op)

1.7.4 The convolution theorem

Theorem 1.2 (Convolution theorem). With \(\mathcal{F}\) denoting the Fourier transform (Section 1.9),

\[\begin{equation} \mathcal{F}\{x * h\} = \mathcal{F}\{x\}\cdot\mathcal{F}\{h\}, \qquad \mathcal{F}\{x \cdot h\} = \mathcal{F}\{x\} * \mathcal{F}\{h\}. \tag{36} \end{equation}\]

Convolution in one domain is multiplication in the other. Four consequences drive practice:

1. Computational efficiency. Direct convolution of two \(N\)-point signals costs \(O(N^2)\); via the FFT it costs \(O(N\log N)\). For a \(512\times512\) image the saving is several orders of magnitude.
2. Filter design. Rather than hand-crafting a spatial kernel, specify the desired frequency response (low-pass to suppress noise, band-pass to isolate a physiological rhythm) and inverse-transform it.
3. Deconvolution. If the impulse response is known, \(\hat{x}(\omega) = \hat{y}(\omega)/\hat{h}(\omega)\) recovers the input. This is the basis of deconvolution microscopy and image deblurring, and it is ill-posed wherever \(\hat{h}(\omega)\) is small, which is why regularization is unavoidable (Chapters 7 and 8).
4. Transfer functions. \(H(\omega) = \mathcal{F}\{h(t)\}\) completely characterizes an LTI system: \(\lvert H(\omega)\rvert\) is the gain at each frequency and \(\arg H(\omega)\) the phase shift.

set.seed(42)
n   <- 256; tt <- 0:(n-1)
sig <- sin(2*pi*5*tt/n) + 0.5*sin(2*pi*20*tt/n) + 0.2*rnorm(n)
k   <- dnorm(seq(-3, 3, length.out = 21)); k <- k/sum(k)
kp  <- rep(0, n); kp[seq_along(k)] <- k

## Method 1: convolution by its definition, O(n^2)
circ_conv <- function(x, h) {
  n <- length(x); y <- numeric(n)
  for (i in 0:(n-1)) {
    s <- 0
    for (j in 0:(n-1)) s <- s + x[j+1] * h[((i - j) %% n) + 1]
    y[i+1] <- s
  }
  y
}

## Method 2: multiply the transforms, O(n log n)
y_direct <- circ_conv(sig, kp)
y_fft    <- Re(fft(fft(sig) * fft(kp), inverse = TRUE))/n

cat("maximum absolute difference:", format(max(abs(y_direct - y_fft)), digits = 3), "\n")

## maximum absolute difference: 1.33e-15

t_direct <- system.time(circ_conv(sig, kp))[["elapsed"]]
t_fft    <- system.time(for (i in 1:200) Re(fft(fft(sig)*fft(kp), inverse = TRUE))/n)[["elapsed"]]/200
cat(sprintf("direct: %.4f s per call;  FFT: %.6f s per call;  speed-up ~%.0fx\n",
            t_direct, t_fft, t_direct/t_fft))

## direct: 0.0000 s per call;  FFT: 0.000100 s per call;  speed-up ~0x

The two methods agree to machine precision, and the FFT route is faster by a factor that grows with \(N\).

1.7.5 The point spread function and instrument transfer

The point spread function (PSF) \(h(x,y)\) is the two-dimensional impulse response of an imaging system: the image it forms of an ideal point source. The measured image is therefore

\[\begin{equation} g(x,y) = (f * h)(x,y) + n(x,y), \tag{37} \end{equation}\]

where \(f\) is the true object and \(n\) is noise (Section 1.15). Equation (37) is the imaging equation, and it recurs in every modality chapter:

• Ultrasound (Chapter 3): the beam profile is a spatial PSF; axial resolution is set by pulse length, lateral resolution by beam width.
• X-ray and CT (Chapter 5): focal-spot size, detector aperture, and the reconstruction filter together define the PSF.
• PET and SPECT (Chapter 6): positron range, photon acollinearity, and collimator geometry blur the true activity distribution.
• MRI (Chapter 4): the \(k\)-space sampling extent and any apodization window define an effective PSF; truncation produces the ringing of Section 1.8.
• ECG amplifiers: the impulse response of the filter chain must be flat across roughly 0.05–150 Hz or the QRS morphology is distorted.

From equation to imaging. The imaging equation \(g = f * h + n\) says that what you see is what is there, blurred by the instrument, plus noise. Every resolution claim is a statement about \(h\); every denoising method is a statement about \(n\); and every attempt to recover \(f\) from \(g\) is an inverse problem whose difficulty is set by how small \(\hat h(\omega)\) becomes at the frequencies you care about. This single equation organizes Chapters 3 through 7.

Section 1.7 summary.

• A linear time-invariant system is completely characterized by its impulse response \(h\).
• The output is the convolution \(y = x * h\): a weighted superposition of the input.
• In two dimensions, convolution is image filtering; Gaussian, Laplacian, and Sobel kernels smooth, sharpen, and detect edges.
• The convolution theorem converts convolution into multiplication in the frequency domain, giving an \(O(N\log N)\) algorithm and a framework for filter design and deconvolution.
• The PSF is the imaging system’s two-dimensional impulse response, and \(g = f*h + n\) is the imaging equation.

Checkpoint 1.8. Two boxcar functions of unit height and unit width are convolved. Sketch the result without computing it, and explain how you knew its width and shape from the convolution property alone.

1.8 Fourier Series

Before generalizing to arbitrary signals, we treat the periodic case. Fourier series are historically first, conceptually simpler, and directly applicable to the many biological processes that repeat: the cardiac cycle, the respiratory cycle, circadian rhythms, gait, and periodic neural firing.

1.8.1 The trigonometric series

Theorem 1.3 (Fourier series). A sufficiently regular periodic function \(f(t)\) with period \(T\), fundamental frequency \(f_0 = 1/T\), and angular frequency \(\omega_0 = 2\pi/T\) can be expanded as

\[\begin{equation} f(t) = \frac{a_0}{2} + \sum_{n=1}^{\infty}\Big[a_n\cos(n\omega_0 t) + b_n\sin(n\omega_0 t)\Big], \tag{38} \end{equation}\]

with coefficients

\[\begin{equation} a_n = \frac{2}{T}\int_0^{T} f(t)\cos(n\omega_0 t)\,dt, \qquad b_n = \frac{2}{T}\int_0^{T} f(t)\sin(n\omega_0 t)\,dt. \tag{39} \end{equation}\]

The term \(a_0/2\) is the DC (mean) level. The \(n\)-th harmonic oscillates at \(n f_0\).

The coefficient formulas are not arbitrary: they are projections, and they work because the sinusoids form an orthogonal basis on \([0,T]\):

\[\begin{equation} \int_0^{T}\cos(n\omega_0 t)\cos(m\omega_0 t)\,dt = \begin{cases} 0, & n \ne m\\ T/2, & n = m \ne 0\\ T, & n = m = 0\end{cases} \qquad \int_0^{T}\cos(n\omega_0 t)\sin(m\omega_0 t)\,dt = 0 \ \ \forall\, n,m, \tag{40} \end{equation}\]

with the analogous relation for \(\sin\cdot\sin\). Orthogonality is the same idea as the dot product of Section 1.3: each coefficient is the inner product of the signal with one basis function, divided by that basis function’s squared norm. The Fourier coefficient is a projection in exactly the sense that a component of a vector is a projection.

1.8.2 The complex exponential form

Euler’s formula ((26)) compresses the two real series into one complex series:

\[\begin{equation} f(t) = \sum_{n=-\infty}^{\infty} c_n e^{i n\omega_0 t}, \qquad c_n = \frac{1}{T}\int_0^{T} f(t)\,e^{-i n \omega_0 t}\,dt, \tag{41} \end{equation}\]

with \(c_0 = a_0/2\), \(c_n = (a_n - i b_n)/2\) and \(c_{-n} = (a_n + i b_n)/2\) for \(n>0\). For real \(f\), \(c_{-n} = c_n^{*}\), which is why a real signal’s spectrum is conjugate-symmetric and why only half of it carries independent information (a fact exploited for acceleration in Chapter 4).

Parseval’s theorem states that total power is conserved between domains:

\[\begin{equation} \frac{1}{T}\int_0^{T}\lvert f(t)\rvert^{2}\,dt = \sum_{n=-\infty}^{\infty}\lvert c_n\rvert^{2} = \frac{a_0^2}{4} + \frac{1}{2}\sum_{n=1}^{\infty}\left(a_n^2 + b_n^2\right). \tag{42} \end{equation}\]

Power computed in time equals power computed in frequency. This is the licence for spectral analysis of physiological signals: partitioning the right-hand sum into frequency bands partitions the actual signal power, which is what heart-rate-variability and EEG band-power analyses assume.

1.8.3 Convergence and the Gibbs phenomenon

Truncating the series after \(K\) harmonics produces the best \(K\)-term least-squares approximation, but near a discontinuity the partial sums overshoot by approximately \(8.95\%\) of the jump, and this overshoot does not diminish as \(K\) increases; it merely narrows. This is the Gibbs phenomenon.

The consequence for imaging is direct and clinically visible. MRI acquires a finite region of \(k\)-space, which is a truncation of a Fourier series, so sharp tissue boundaries acquire parallel ringing bands, the Gibbs (truncation) artifact of Chapter 4. It is mitigated, never removed, by apodization with a smooth window (Hamming, Hanning, Kaiser), at the cost of resolution.

tt <- seq(0, 2, length.out = 4000)
partial <- function(K) {
  s <- 0
  for (m in seq(1, K, by = 2)) s <- s + (4/(pi*m))*sin(2*pi*m*tt)
  s
}
ov <- sapply(c(9, 49, 199), function(K) max(partial(K)) - 1)
knitr::kable(data.frame(harmonics = c(9, 49, 199),
                        max_overshoot = round(ov, 4),
                        percent = paste0(round(100*ov, 2), "%")),
             caption = "Gibbs overshoot is nearly independent of the number of harmonics.")

Table 4: Gibbs overshoot is nearly independent of the number of harmonics.

harmonics	max_overshoot	percent
9	0.1823	18.23%
49	0.1791	17.91%
199	0.1790	17.9%

gp <- rbind(data.frame(t = tt, y = sign(sin(2*pi*tt)), s = "square wave"),
            data.frame(t = tt, y = partial(9),   s = "K = 9"),
            data.frame(t = tt, y = partial(49),  s = "K = 49"),
            data.frame(t = tt, y = partial(199), s = "K = 199"))
ggplot(gp, aes(t, y, colour = s)) +
  geom_line(linewidth = 0.7) +
  scale_colour_manual(values = c("square wave" = "black", "K = 9" = bpad_pal[1],
                                 "K = 49" = bpad_pal[2], "K = 199" = bpad_pal[3])) +
  coord_cartesian(xlim = c(0.35, 0.65), ylim = c(-1.3, 1.3)) +
  labs(title = "Gibbs overshoot near a discontinuity",
       x = "t", y = "amplitude")

Figure 15: Figure 1.15: The Gibbs phenomenon. Partial Fourier sums of a square wave overshoot the discontinuity by about 8.95 percent regardless of how many harmonics are retained; the overshoot narrows but does not shrink.

1.8.4 Biomedical application: harmonic analysis of the ECG

The ECG is approximately periodic at the heart rate. Its Fourier decomposition separates the information carried by the rate from the information carried by the waveform shape:

• \(c_{\pm 1}\) corresponds to the fundamental, the heart rate itself;
• higher harmonics \(c_{\pm 2}, c_{\pm 3},\ldots\) encode the morphology: the sharp QRS complex requires many harmonics, whereas the broad T wave requires few;
• in heart-rate variability analysis, power in the 0.04–0.15 Hz band (low frequency) and the 0.15–0.40 Hz band (high frequency, respiratory) provides a non-invasive probe of autonomic cardiovascular control.

trapz <- function(x, y) sum(diff(x) * (head(y, -1) + tail(y, -1))/2)

f_cardiac <- 1.2                    # Hz, that is 72 beats per minute
Tp <- 1/f_cardiac
ecg_fun <- function(t, Tp) {        # P wave, QRS complex, T wave as Gaussians
  tn  <- t %% Tp
  gau <- function(mu, s, a) a*exp(-((tn - mu*Tp)^2)/(2*(s*Tp)^2))
  gau(0.15, 0.030,  0.10) +                              # P  (atrial depolarization)
  gau(0.35, 0.015, -0.15) + gau(0.38, 0.020, 1.00) +
  gau(0.42, 0.018, -0.25) +                              # QRS (ventricular depol.)
  gau(0.65, 0.060,  0.20)                                # T  (ventricular repol.)
}

tg  <- seq(0, 3*Tp, length.out = 1500); ecg <- ecg_fun(tg, Tp)
Nh  <- 20; w0 <- 2*pi/Tp
tp  <- seq(0, Tp, length.out = 4000); yp <- ecg_fun(tp, Tp)
nn  <- -Nh:Nh
cn  <- sapply(nn, function(m) trapz(tp, yp*exp(-1i*m*w0*tp))/Tp)   # Eq. \@ref(eq:eq1-41)

recon <- function(K) {
  idx <- which(abs(nn) <= K)
  Re(colSums(cn[idx] * exp(1i * outer(nn[idx], tg) * w0)))
}

cat(sprintf("Parseval check:  time-domain power = %.6f,  spectral power = %.6f\n",
            trapz(tp, yp^2)/Tp, sum(Mod(cn)^2)))

## Parseval check:  time-domain power = 0.032902,  spectral power = 0.032892

for (K in c(3, 7, 15, 20))
  cat(sprintf("  K = %2d harmonics -> RMSE = %.4f\n", K,
              sqrt(mean((recon(K) - ecg)^2))))

##   K =  3 harmonics -> RMSE = 0.1476
##   K =  7 harmonics -> RMSE = 0.1021
##   K = 15 harmonics -> RMSE = 0.0208
##   K = 20 harmonics -> RMSE = 0.0033

rec_df <- rbind(data.frame(t = tg, y = ecg,       s = "original"),
                data.frame(t = tg, y = recon(3),  s = "K = 3"),
                data.frame(t = tg, y = recon(7),  s = "K = 7"),
                data.frame(t = tg, y = recon(15), s = "K = 15"))
p1 <- ggplot(rec_df, aes(t, y, colour = s)) +
  geom_line(linewidth = 0.75) +
  scale_colour_manual(values = c("original" = "black", "K = 3" = bpad_pal[1],
                                 "K = 7" = bpad_pal[2], "K = 15" = bpad_pal[3])) +
  labs(title = "Fourier-series reconstruction of a synthetic ECG",
       x = "time (s)", y = "amplitude (mV)")

sp <- data.frame(harmonic = 0:Nh, magnitude = Mod(cn[(Nh+1):(2*Nh+1)]))
p2 <- ggplot(sp, aes(harmonic, magnitude)) +
  geom_col(fill = bpad_pal[1], width = 0.65) +
  labs(title = "Harmonic magnitude spectrum |c_n|",
       x = "harmonic number n  (frequency = n x 1.2 Hz)", y = "|c_n|")

bpad_grid(p1, p2, ncol = 1)

Figure 16: Figure 1.16: Fourier-series reconstruction of a synthetic ECG. Few harmonics recover the slow P and T waves but blunt the QRS complex; recovering the sharp QRS requires the high harmonics, which is exactly why ECG amplifiers must pass frequencies well above the heart rate.

Section 1.8 summary.

• A periodic signal is a sum of harmonics at integer multiples of its fundamental frequency.
• Fourier coefficients are projections onto an orthogonal basis, the same operation as a vector dot product.
• Parseval’s theorem conserves power between the time and frequency domains, which licenses band-power analysis.
• Truncating the series produces Gibbs ringing of about 9% near discontinuities: the origin of MRI truncation artifact.
• Sharp features (the QRS complex) need high harmonics; smooth features (the T wave) do not.

Checkpoint 1.9. An ECG amplifier is band-limited to 0–40 Hz. Using Figure 1.16, predict qualitatively which parts of the waveform survive and which are distorted, and state the clinical risk.

1.9 The Fourier Transform

Letting the period \(T\to\infty\) in a Fourier series turns the discrete harmonics into a continuous spectrum, giving the Fourier transform. It decomposes (analyzes) a signal into frequency components, and reconstructs (synthesizes) the signal from them. It is the single most-used mathematical tool in the rest of this textbook: MRI acquires data directly in the Fourier domain (Chapter 4), CT reconstruction rests on a Fourier relationship between projections and the object (Chapter 5), FTIR recovers a spectrum from an interferogram (Chapter 2), and OCT recovers a depth profile from a spectral interferogram (Chapter 2).

1.9.1 Properties of sine and cosine, and orthogonality

The function \(\sin(2\pi\omega t/n)\) is periodic with period \(n/\omega\) and frequency \(\omega/n\); it is odd, \(\sin(-\theta)=-\sin\theta\), with \(\sin 0 = 0\). The cosine is even, \(\cos(-\theta)=\cos\theta\), with \(\cos 0 = 1\), and the two are related by a quarter-period shift, \(\sin(\theta + \pi/2) = \cos\theta\). Together with \(\sin^2\theta + \cos^2\theta = 1\) and the angle-addition identities, these are all the trigonometry the chapter requires.

Two summation identities make the discrete theory work. For integer \(\omega\),

\[\begin{equation} \sum_{t=0}^{n-1}\sin\!\left(2\pi\omega\frac{t}{n}\right) = 0, \qquad \sum_{t=0}^{n-1}\cos\!\left(2\pi\omega\frac{t}{n}\right) = \begin{cases} n, & \omega \equiv 0 \ (\mathrm{mod}\ n)\\ 0, & \text{otherwise}\end{cases} \tag{43} \end{equation}\]

and sinusoids at different frequencies are orthogonal:

\[\begin{equation} \sum_{t=1}^{n}\sin\!\left(2\pi\omega_1\frac{t}{n}\right)\sin\!\left(2\pi\omega_2\frac{t}{n}\right) = \begin{cases} n/2, & \omega_1 = \omega_2 \ne 0\\ -n/2, & \omega_1 = -\omega_2 \ne 0\\ 0, & \text{otherwise.}\end{cases} \tag{44} \end{equation}\]

Orthogonality is also where the limits of resolution originate. A very high-frequency sinusoid completes a full cycle between consecutive samples and becomes indistinguishable from a lower-frequency one (this is aliasing, Section 1.10). A very low-frequency sinusoid does not complete even half a cycle within the record and becomes indistinguishable from a constant plus drift (this is why long-duration physiological trends require long recordings).

1.9.2 The discrete Fourier transform

Definition 1.6 (DFT and inverse DFT). For a sampled signal \(x(t)\), \(t = 0,\ldots,n-1\),

\[\begin{equation} \hat{x}(\omega) = \sum_{t=0}^{n-1} x(t)\,e^{-2\pi i t\omega/n}, \qquad \omega = 0,\ldots,n-1, \tag{45} \end{equation}\]

\[\begin{equation} x(t) = \frac{1}{n}\sum_{\omega=0}^{n-1}\hat{x}(\omega)\,e^{+2\pi i t\omega/n}. \tag{46} \end{equation}\]

Each \(\hat{x}(\omega) = r e^{i\phi}\) is a phasor: \(r = \lvert\hat x(\omega)\rvert\) is the amplitude at that frequency and \(\phi = \arg\hat{x}(\omega)\) is its phase. The fast Fourier transform (FFT) computes ((45)) in \(O(n\log n)\) instead of \(O(n^2)\); in R both directions are provided by fft(), with inverse = TRUE supplying ((46)) up to the factor \(1/n\), which R does not apply automatically.

Two properties are used repeatedly:

• Differentiation becomes multiplication. Since \(\frac{d}{dt}e^{2\pi i\omega t/n} = \frac{2\pi i\omega}{n}e^{2\pi i \omega t/n}\), differentiating a signal multiplies its transform by \(2\pi i \omega/n\). Differential equations become algebraic equations, which is how Fourier methods solve PDEs.
• Shifts become phase ramps. If \(y(t) = x(t - t_0)\) then \(\hat{y}(\omega) = \hat{x}(\omega)e^{-2\pi i t_0 \omega/n}\): a translation changes only the phase, never the magnitude. This is why magnitude-only images lose position information, and why phase-based registration (Chapter 7) works.

n <- 128; t <- 0:(n - 1)
y <- 1.5*sin(2*pi*1*t/n) + 0.5*sin(2*pi*5*t/n)

Y   <- fft(y)
mag <- Mod(Y); phase <- Arg(Y)
rec <- Re(fft(Y, inverse = TRUE))/n

cat("spectral peaks at k =", which(mag[1:(n/2)] > 1) - 1, "\n")

## spectral peaks at k = 1 5

cat("maximum reconstruction error:", format(max(abs(rec - y)), digits = 3), "\n")

## maximum reconstruction error: 6.66e-16

cat(sprintf("Parseval: sum|x|^2 = %.4f,  sum|X|^2/n = %.4f\n",
            sum(y^2), sum(mag^2)/n))

## Parseval: sum|x|^2 = 160.0000,  sum|X|^2/n = 160.0000

p1 <- ggplot(data.frame(t = t, y = y), aes(t, y)) +
  geom_line(colour = bpad_pal[1], linewidth = 0.8) +
  labs(title = "signal x(t)", x = "sample t", y = "amplitude")
p2 <- ggplot(data.frame(k = 0:(n/2 - 1), m = mag[1:(n/2)]), aes(k, m)) +
  geom_segment(aes(xend = k, yend = 0), colour = bpad_pal[2], linewidth = 0.6) +
  coord_cartesian(xlim = c(0, 20)) +
  labs(title = "magnitude spectrum |X(k)|", x = "frequency index k", y = "|X(k)|")
p3 <- ggplot(data.frame(t = t, y = y, r = rec), aes(t)) +
  geom_line(aes(y = y), colour = "black", linewidth = 0.9) +
  geom_point(aes(y = r), colour = bpad_pal[3], size = 0.9) +
  labs(title = "inverse DFT (points) over original (line)",
       x = "sample t", y = "amplitude")
bpad_grid(p1, p2, p3, ncol = 3)

Figure 17: Figure 1.17: Discrete Fourier analysis and synthesis. Left: a two-tone signal. Centre: its magnitude spectrum, with peaks at exactly the two constituent frequencies. Right: the inverse transform recovers the original signal to machine precision.

1.9.3 The power spectrum and the Wiener–Khinchin theorem

For a stationary stochastic process \(x(t)\) with mean \(\mu\), the autocovariance at lag \(h\) is

\[\gamma(h) = \mathrm{cov}\big(x(t),\,x(t+h)\big) = \mathbb{E}\big[(x(t)-\mu)(x(t+h)-\mu)\big].\]

Because \(x\) is random, so are \(\hat{x}(\omega)\) and \(\lvert\hat{x}(\omega)\rvert^{2}\). Their expectation, however, is a deterministic real function, the power spectrum

\[S(\omega) = \mathbb{E}\big[\lvert\hat{x}(\omega)\rvert^{2}\big].\]

Theorem 1.4 (Wiener–Khinchin). The power spectrum is the Fourier transform of the autocovariance,

\[\begin{equation} S(\omega) = \hat{\gamma}(\omega) = \sum_{h}\gamma(h)\,e^{-2\pi i h\omega/n}. \tag{47} \end{equation}\]

Estimating either quantity therefore determines the other. This equivalence underlies spectral estimation in EEG, resting-state fMRI connectivity, heart-rate variability, and the noise-power-spectrum characterization of CT and MRI detectors in Chapter 7.

1.9.4 The continuous Fourier transform

For an integrable \(f:\mathbb{R}\to\mathbb{C}\),

\[\begin{equation} \hat{f}(\omega) = \mathcal{F}\{f\}(\omega) = \int_{-\infty}^{\infty} f(x)\,e^{-i2\pi\omega x}\,dx, \qquad f(x) = \int_{-\infty}^{\infty}\hat{f}(\omega)\,e^{+i2\pi\omega x}\,d\omega, \tag{48} \end{equation}\]

with the multivariate generalization \(\hat{g}(\boldsymbol{\omega}) = \int_{\mathbb{R}^n} g(\mathbf{x})e^{-i2\pi\langle\boldsymbol{\omega},\mathbf{x}\rangle}d\mathbf{x}\).

The properties below are used, explicitly or implicitly, in every subsequent chapter.

Property	Statement	Where it is used
Linearity	\(\mathcal{F}\{\alpha f + \beta g\} = \alpha\hat f + \beta\hat g\)	superposition of sources
Translation	\(h(x)=f(x-a) \Rightarrow \hat h(\omega) = e^{-2\pi i a\omega}\hat f(\omega)\)	motion, registration
Scaling	\(h(x)=f(ax) \Rightarrow \hat h(\omega) = \frac{1}{\lvert a\rvert}\hat f(\omega/a)\)	resolution–bandwidth trade-off
Differentiation	\(\mathcal{F}\{f^{(k)}\}(\omega) = (2\pi i\omega)^{k}\hat f(\omega)\)	solving PDEs; ramp filter in CT
Convolution	\(\widehat{f*g} = \hat f\cdot\hat g\) and \(\widehat{f\cdot g} = \hat f * \hat g\)	filtering, PSF, deconvolution
Parseval–Plancherel	\(\int\lvert f\rvert^{2}dx = \int\lvert\hat f\rvert^{2}d\omega\)	energy and SNR budgets
Duality	\(\mathcal{F}\{\mathcal{F}\{f\}\}(\omega) = f(-\omega)\)	symmetry arguments
Impulse	\(\mathcal{F}\{\delta(x - a)\}(\omega) = e^{-2\pi i a\omega}\)	point sources, PSF measurement
Gaussian	\(\mathcal{F}\{e^{-\pi x^2}\} = e^{-\pi\omega^2}\)	the Gaussian is its own transform

The last entry is worth dwelling on: the Gaussian is an eigenfunction of the Fourier transform, and a narrow Gaussian in space transforms to a wide Gaussian in frequency. That reciprocal width relation is the mathematical form of the uncertainty principle, and it is why a short ultrasound pulse necessarily has a wide bandwidth (Chapter 3), why a broad optical spectrum gives fine OCT axial resolution (Chapter 2), and why a sharply localized PSF requires a wide \(k\)-space coverage (Chapter 4).

1.9.5 The two-dimensional Fourier transform

For an image \(f(x,y)\),

\[\begin{equation} F(u,v) = \int\!\!\int f(x,y)\,e^{-i2\pi(ux+vy)}\,dx\,dy, \qquad f(x,y) = \int\!\!\int F(u,v)\,e^{+i2\pi(ux+vy)}\,du\,dv, \tag{49} \end{equation}\]

where \(u,v\) are spatial frequencies. Writing \(F = F_R + iF_I\), the magnitude spectrum is \(\lvert F\rvert = \sqrt{F_R^2+F_I^2}\) and the phase spectrum is \(\arctan(F_I/F_R)\).

Geometrically, the kernel \(e^{-i2\pi(ux+vy)} = \cos\big(2\pi\langle \mathbf{U},\mathbf{X}\rangle\big) + i\sin\big(2\pi\langle\mathbf{U},\mathbf{X}\rangle\big)\), with \(\mathbf{U} = (u,v)^{T}\) and \(\mathbf{X}=(x,y)^{T}\), describes a set of parallel plane waves with normal direction \(\mathbf{U}\) and wavelength \(1/\sqrt{u^2+v^2}\). Every image is a superposition of such gratings.

Four applications organize later chapters:

1. MRI reconstruction (Chapter 4). The scanner samples \(k\)-space, the 2-D Fourier domain, directly; the image is the 2-D inverse FFT. The centre of \(k\)-space carries contrast, the periphery carries edges.
2. CT reconstruction (Chapter 5). The central slice theorem states that the 1-D Fourier transform of a projection at angle \(\theta\) equals a radial line at angle \(\theta\) through the 2-D Fourier transform of the object, which is what makes filtered back-projection work.
3. Spatial-frequency filtering (Chapter 7). Low frequencies carry gross structure and contrast; high frequencies carry edges, texture, and noise. Every denoising/sharpening trade-off is a choice about which to keep.
4. Phase contrast (Chapter 2). Phase encodes refractive-index variation, enabling visualization of weakly absorbing structures that produce almost no amplitude contrast.

fftshift2 <- function(M) {                  # move zero frequency to the centre
  d <- dim(M); r <- d[1]; cc <- d[2]
  M[c((r %/% 2 + 1):r, 1:(r %/% 2)), c((cc %/% 2 + 1):cc, 1:(cc %/% 2))]
}

set.seed(42)
N   <- 128
rad <- outer(1:N, 1:N, function(i, j) sqrt((i - N/2)^2 + (j - N/2)^2))
cell <- ifelse(rad < 40, 0.6 + 0.3*sin(2*pi*rad/15), 0)          # "cell body"
nuc  <- outer(1:N, 1:N, function(i, j) sqrt((i - N/2)^2 + (j - N/2 - 8)^2))
cell <- cell + 0.4*exp(-nuc^2/50)*(nuc < 15)                     # "nucleus"
cell <- cell + matrix(rnorm(N*N, 0, 0.03), N)

Fc  <- fft(cell); Fs <- fftshift2(Fc)
mag <- log1p(Mod(Fs)); pha <- Arg(Fs)

keep <- 20                                   # central block of frequencies
mask <- matrix(0, N, N)
ctr  <- (N/2 - keep/2 + 1):(N/2 + keep/2)
mask[ctr, ctr] <- 1
low  <- Re(fft(fftshift2(Fs*mask), inverse = TRUE))/(N*N)

cat(sprintf("fraction of k-space retained: %.2f%%;  correlation with original: %.3f\n",
            100*sum(mask)/N^2, cor(c(cell), c(low))))

## fraction of k-space retained: 2.44%;  correlation with original: 0.990

op <- par(mfrow = c(2, 2), mar = c(0.6, 0.6, 2.2, 0.6))
show_img(cell, "object f(x,y)")
show_img(mag,  "log(1 + |F(u,v)|)")
show_img(pha,  "phase of F(u,v)")
show_img(low,  sprintf("reconstruction from %d x %d central frequencies", keep, keep))

Figure 18: Figure 1.18: Two-dimensional Fourier analysis of a synthetic cell image. The concentric rings in the object appear as a ring in the magnitude spectrum. Retaining only a 20 x 20 central block of frequencies (low-pass filtering) preserves gross shape while discarding fine texture.

par(op)

Section 1.9 summary.

• The Fourier transform extends Fourier series to non-periodic signals; the DFT is its sampled version and the FFT its fast algorithm.
• Every transform value is a phasor: magnitude and phase. Translation alters phase only, never magnitude.
• Differentiation in one domain is multiplication by \(2\pi i\omega\) in the other, which turns differential equations into algebra.
• Wiener–Khinchin identifies the power spectrum with the transform of the autocovariance.
• The 2-D transform is the reconstruction engine for MRI (\(k\)-space) and CT (central slice theorem), and the basis of spatial-frequency filtering.
• A narrow function has a wide transform: the reciprocal-width relation behind every resolution–bandwidth trade-off in the book.

Checkpoint 1.10. In Figure 1.18, only 2.4% of the \(k\)-space coefficients were retained, yet the reconstruction still resembles the object. Which image features were lost, and what does this predict about accelerated MRI acquisition?

1.10 Sampling, Aliasing, and the Nyquist Criterion

Every digital biomedical measurement is a sampled measurement. The sampling theorem states exactly when sampling loses nothing, and aliasing is what happens when the condition is violated. Because the criterion is a statement about the Fourier transform, it applies identically to sampling in time (ECG, EEG, Doppler) and in space (CT, MRI, microscopy, digital pathology).

1.10.1 The sampling theorem

Theorem 1.5 (Nyquist–Shannon). A continuous signal \(x(t)\) that is bandlimited, meaning \(\hat{x}(f) = 0\) for \(\lvert f\rvert > f_{\max}\), is completely determined by its samples taken at uniform intervals \(T_s = 1/f_s\) provided

\[\begin{equation} f_s > 2 f_{\max}. \tag{50} \end{equation}\]

The rate \(2f_{\max}\) is the Nyquist rate; \(f_{\mathrm{Nyq}} = f_s/2\) is the Nyquist frequency, the highest frequency that can be represented unambiguously.

Proof sketch. Sampling multiplies the signal by an impulse train. By the convolution theorem ((36)), multiplication in time is convolution in frequency, so the spectrum is replicated at intervals of \(f_s\). If \(f_s > 2f_{\max}\) the replicas do not overlap and an ideal low-pass filter recovers the original spectrum exactly. If \(f_s \le 2f_{\max}\) they overlap, and the overlap is irreversible. \(\square\)

1.10.2 Aliasing

When ((50)) fails, a component at \(f > f_{\mathrm{Nyq}}\) is folded into the observable band and appears at

\[\begin{equation} f_{\mathrm{alias}} = \left\lvert f - k f_s \right\rvert, \tag{51} \end{equation}\]

with \(k\) the integer that places \(f_{\mathrm{alias}}\) in \([0, f_s/2]\). The critical point is that an aliased component is indistinguishable from a genuine component at the alias frequency; no post-processing can separate them.

f0 <- 5; t_max <- 1
t_fine <- seq(0, t_max, length.out = 2000)

alias_of <- function(f, fs) abs(f - fs*round(f/fs))

panels <- lapply(c(20, 10, 8, 6), function(fs) {
  ts <- seq(0, t_max, by = 1/fs)
  fa <- alias_of(f0, fs)
  lab <- if (fs > 2*f0) "adequate" else if (fs == 2*f0) "critical" else "aliased"
  gg <- ggplot() +
    geom_line(data = data.frame(t = t_fine, y = cos(2*pi*f0*t_fine)),
              aes(t, y), colour = "grey35", linewidth = 0.6) +
    geom_point(data = data.frame(t = ts, y = cos(2*pi*f0*ts)),
               aes(t, y), colour = bpad_pal[2], size = 2.2) +
    labs(title = sprintf("fs = %d Hz (%s), apparent f = %.0f Hz", fs, lab, fa),
         x = "time (s)", y = "amplitude") +
    coord_cartesian(ylim = c(-1.4, 1.4))
  if (fs < 2*f0)
    gg <- gg + geom_line(data = data.frame(t = t_fine, y = cos(2*pi*fa*t_fine)),
                         aes(t, y), colour = bpad_pal[1],
                         linetype = "dashed", linewidth = 0.8)
  gg
})

knitr::kable(data.frame(sampling_rate_Hz = c(20, 10, 8, 6),
                        Nyquist_Hz = c(20, 10, 8, 6)/2,
                        apparent_frequency_Hz =
                          sapply(c(20, 10, 8, 6), function(fs) alias_of(f0, fs))),
             caption = "Apparent frequency of a 5 Hz signal at four sampling rates.")

Table 5: Apparent frequency of a 5 Hz signal at four sampling rates.

sampling_rate_Hz	Nyquist_Hz	apparent_frequency_Hz
20	10	5
10	5	5
8	4	3
6	3	1

bpad_grid(panels[[1]], panels[[2]], panels[[3]], panels[[4]], ncol = 2)

Figure 19: Figure 1.19: Sampling a 5 Hz cosine at four rates. Above Nyquist the samples determine the signal. At exactly Nyquist the reconstruction is ambiguous. Below Nyquist the samples are consistent with a lower-frequency signal (dashed), and the true frequency is unrecoverable.

Aliasing across the modalities.

• ECG and EEG. Clinical ECG is sampled at 500–1000 Hz because clinically relevant QRS content extends to roughly 150 Hz. Sampling at 250 Hz would fold high-frequency content into the diagnostic band and could distort or mask ventricular arrhythmia morphology.
• Doppler ultrasound (Chapter 3). The pulse repetition frequency sets the Nyquist velocity; flow faster than this wraps around and reverses direction on the display, the familiar Doppler aliasing artifact.
• MRI (Chapter 4). Undersampling in the phase-encode direction folds anatomy outside the field of view back into the image (wrap-around). Parallel imaging (SENSE, GRAPPA) uses multiple receive coils to unfold it, deliberately trading a known aliasing pattern for speed.
• Optical microscopy (Chapter 2). The pixel pitch must satisfy \(\Delta x \le \lambda/(4\,\mathrm{NA})\) relative to the optical cutoff frequency; violating this discards resolution the optics already delivered.

1.10.3 Anti-aliasing and reconstruction

An anti-aliasing filter removes content above \(f_s/2\) before sampling. It must be analog (or applied at a higher intermediate rate), because once frequencies have folded, no digital filter can separate them from genuine signal. In practice:

• real filters have finite roll-off, so sampling at exactly \(2f_{\max}\) is inadequate; oversampling by a factor of 3–5 provides a guard band;
• oversampling followed by digital decimation improves effective resolution and shapes the noise spectrum, the principle of sigma-delta converters.

Perfect reconstruction of a bandlimited signal from adequate samples is given by Whittaker–Shannon interpolation,

\[\begin{equation} x(t) = \sum_{n=-\infty}^{\infty} x[n]\,\mathrm{sinc}\!\left(\frac{t - nT_s}{T_s}\right), \qquad \mathrm{sinc}(u) = \frac{\sin(\pi u)}{\pi u}, \tag{52} \end{equation}\]

the sinc being the impulse response of the ideal low-pass filter. Practical interpolators (linear, cubic, Lanczos) trade some fidelity for locality and speed; the choice matters in Chapter 7, where every resampling step in a registration pipeline applies one.

Spatial sampling rules of thumb.

• CT and MRI: with field of view FOV and matrix size \(N\), the sample spacing is \(\Delta x = \mathrm{FOV}/N\); structures finer than \(2\Delta x\) cannot be resolved.
• Digital pathology: nuclear detail is 5–10 µm, so whole-slide imaging typically requires \(\le 0.5\) µm per pixel.
• Micro-CT of trabecular bone (trabecular thickness \(\approx 100\) µm) requires voxels well below 50 µm.

Common misconception. “Aliasing can be removed by filtering the sampled data.” It cannot. Filtering after sampling removes the band containing the alias, and with it any real signal at that frequency. Anti-aliasing is a design decision made before the analog-to-digital converter, not a post-processing option.

Section 1.10 summary.

• A bandlimited signal is fully determined by samples taken faster than twice its highest frequency.
• Below the Nyquist rate, high frequencies fold to \(\lvert f - kf_s\rvert\) and become permanently indistinguishable from genuine low frequencies.
• Anti-aliasing filters must precede sampling; oversampling provides the necessary guard band.
• The criterion applies identically in time (ECG, Doppler) and space (CT, MRI, microscopy).
• Whittaker–Shannon sinc interpolation is the exact reconstruction; practical interpolators approximate it.

Checkpoint 1.11. A pediatric ECG contains QRS energy to 250 Hz and is digitized at 500 Hz. Compute the Nyquist frequency, state whether the QRS is faithfully captured, and propose a sampling rate and pre-filter with an adequate safety margin.

1.11 Polar, Cylindrical, and Spherical Coordinates

Cartesian coordinates are universal but rarely convenient. Whenever a problem has circular, cylindrical, or spherical symmetry, choosing coordinates aligned with that symmetry reduces the number of variables a solution depends on, and often converts an intractable partial differential equation into a solvable ordinary one.

1.11.1 Two dimensions: polar coordinates

A point is specified by a radius \(r \ge 0\) and an angle \(\theta\) measured counter-clockwise from the positive \(x\)-axis:

\[\begin{equation} x = r\cos\theta, \quad y = r\sin\theta; \qquad r = \sqrt{x^2+y^2}, \quad \theta = \mathrm{atan2}(y,x). \tag{53} \end{equation}\]

The area element is \(dA = r\,dr\,d\theta\); the extra factor of \(r\) is the Jacobian, and forgetting it is the single most common error in polar integration.

1.11.2 Three dimensions: spherical and cylindrical coordinates

Spherical coordinates use a radius \(r\ge 0\), a polar angle \(\theta\in[0,\pi]\) measured from the positive \(z\)-axis, and an azimuthal angle \(\phi\in[0,2\pi)\) measured in the \(xy\)-plane from the positive \(x\)-axis. This is the physics (ISO 80000-2) convention used throughout BPAD:

\[\begin{equation} x = r\sin\theta\cos\phi, \qquad y = r\sin\theta\sin\phi, \qquad z = r\cos\theta, \tag{54} \end{equation}\]

\[\begin{equation} r = \sqrt{x^2+y^2+z^2}, \qquad \theta = \mathrm{atan2}\!\left(\sqrt{x^2+y^2},\, z\right) = \arccos\!\frac{z}{r}, \qquad \phi = \mathrm{atan2}(y,x). \tag{55} \end{equation}\]

The volume element is \(dV = r^2\sin\theta\,dr\,d\theta\,d\phi\).

Notation warning. Many mathematics texts (and some software libraries) swap the roles of \(\theta\) and \(\phi\), using \(\theta\) for the azimuth and \(\phi\) for the polar angle. Always check the convention before using a formula or a library function. In this textbook, \(\theta\) is always the polar (colatitude) angle from \(+z\) and \(\phi\) is always the azimuth.

Cylindrical coordinates retain the polar description in the \(xy\)-plane and add \(z\) unchanged:

\[\begin{equation} x = \rho\cos\phi, \qquad y = \rho\sin\phi, \qquad z = z, \tag{56} \end{equation}\]

with volume element \(dV = \rho\,d\rho\,d\phi\,dz\).

The Laplacian (Section 1.12) takes the following forms, which is precisely why the choice of coordinates matters:

\[\begin{equation} \nabla^2 f = \frac{1}{\rho}\frac{\partial}{\partial\rho}\!\left(\rho\frac{\partial f}{\partial\rho}\right) + \frac{1}{\rho^2}\frac{\partial^2 f}{\partial\phi^2} + \frac{\partial^2 f}{\partial z^2} \qquad\text{(cylindrical)}, \tag{57} \end{equation}\]

\[\begin{equation} \nabla^2 f = \frac{1}{r^2}\frac{\partial}{\partial r}\!\left(r^2\frac{\partial f}{\partial r}\right) + \frac{1}{r^2\sin\theta}\frac{\partial}{\partial\theta}\!\left(\sin\theta\frac{\partial f}{\partial\theta}\right) + \frac{1}{r^2\sin^2\theta}\frac{\partial^2 f}{\partial\phi^2} \qquad\text{(spherical)}. \tag{58} \end{equation}\]

For a spherically symmetric problem, all \(\theta\) and \(\phi\) derivatives vanish and ((58)) collapses to a one-variable operator, which is what makes analytic solutions possible.

Where each system is natural in biomedicine.

System	Natural setting	Chapter
Polar	rotational scanning geometry; cardiac strain in short-axis views; the Radon transform’s \((\theta, s)\) sinogram	5, 7
Cylindrical	dose deposition around a radiotherapy beam axis; a blood vessel or catheter; the spinal cord; a fibre-optic probe	2, 5
Spherical	the eye and retina; the diffusion propagator and \(b\)-vector directions in diffusion MRI; isotropic point sources in emission tomography	2, 4, 6
Prolate spheroidal	ventricular geometry in cardiac modelling	4

cart2pol <- function(x, y)      c(r = sqrt(x^2 + y^2), theta = atan2(y, x))
pol2cart <- function(r, theta)  c(x = r*cos(theta), y = r*sin(theta))

cart2sph <- function(x, y, z)
  c(r = sqrt(x^2 + y^2 + z^2),
    theta = atan2(sqrt(x^2 + y^2), z),        # polar angle from +z
    phi   = atan2(y, x))                      # azimuth in the xy-plane
sph2cart <- function(r, theta, phi)
  c(x = r*sin(theta)*cos(phi), y = r*sin(theta)*sin(phi), z = r*cos(theta))

cart2cyl <- function(x, y, z) c(rho = sqrt(x^2 + y^2), phi = atan2(y, x), z = z)
cyl2cart <- function(rho, phi, z) c(x = rho*cos(phi), y = rho*sin(phi), z = z)

p  <- c(2, 8, -4)                              # a Cartesian point
sp <- cart2sph(p[1], p[2], p[3]); round(sp, 4)

##      r  theta    phi 
## 9.1652 2.0224 1.3258

cy <- cart2cyl(p[1], p[2], p[3]); round(cy, 4)

##     rho     phi       z 
##  8.2462  1.3258 -4.0000

## round-trip verification: conversions must be exact inverses
round(sph2cart(sp["r"], sp["theta"], sp["phi"]), 10)

## x.r y.r z.r 
##   2   8  -4

round(cyl2cart(cy["rho"], cy["phi"], cy["z"]),  10)

## x.rho y.rho   z.z 
##     2     8    -4

theta <- seq(0, 2*pi, length.out = 400)
circ  <- t(sapply(theta, function(a) pol2cart(1, a)))
ggplot(data.frame(x = circ[, 1], y = circ[, 2], theta = theta),
       aes(x, y, colour = theta)) +
  geom_path(linewidth = 1.1) +
  scale_colour_gradientn(colours = bpad_pal[c(1, 3, 5, 2)], name = "theta") +
  coord_fixed() +
  theme(legend.position = "right", legend.title = element_text()) +
  labs(title = "Unit circle traced in polar coordinates",
       x = "x", y = "y")

Figure 20: Figure 1.20: A circle is a one-parameter object in polar coordinates (constant r, varying theta) but requires a two-variable constraint in Cartesian coordinates. Aligning the coordinate system with the symmetry of the problem is what reduces its dimensionality.

Section 1.11 summary.

• Polar, cylindrical, and spherical coordinates are trigonometric re-parameterizations of the same space.
• Volume elements acquire Jacobian factors: \(r\,dr\,d\theta\), \(\rho\,d\rho\,d\phi\,dz\), \(r^2\sin\theta\,dr\,d\theta\,d\phi\).
• The Laplacian is coordinate-dependent; symmetry-aligned coordinates eliminate variables and make analytic solutions possible.
• Conventions for \(\theta\) and \(\phi\) differ between disciplines; always verify before applying a formula.

1.12 Partial Derivatives, Vector Calculus, and PDEs

Almost nothing in biology depends on a single variable. Concentration varies in space and time; temperature varies in three spatial directions; a magnetic field varies with position and with time. Partial derivatives isolate the effect of one variable while the others are held fixed, and partial differential equations relate them.

1.12.1 Partial derivatives and the gradient

For \(f(x,y,z) = 2x + 3xy^3 - 4z\), the partial derivative with respect to \(y\) treats \(x\) and \(z\) as constants, so \(\partial f/\partial y = 9xy^2\). Collecting all first partials gives the gradient

\[\begin{equation} \nabla f = \left(\frac{\partial f}{\partial x},\ \frac{\partial f}{\partial y},\ \frac{\partial f}{\partial z}\right), \tag{59} \end{equation}\]

a vector pointing in the direction of steepest increase, with magnitude equal to that rate. The gradient is the object that gradient descent follows in Chapter 8, and the object that edge detectors estimate in Chapter 7.

Collecting the second partials gives the Hessian \(H_{ij} = \partial^2 f/\partial x_i \partial x_j\), the multivariable generalization of \(f''\). At a critical point (\(\nabla f = 0\)), a positive-definite Hessian identifies a local minimum, a negative-definite Hessian a local maximum, and an indefinite Hessian a saddle point, extending the classification of Section 1.1.2.

1.12.2 The diffusion equation

The most important PDE in biomedical physics is the diffusion equation

\[\begin{equation} \frac{\partial C}{\partial t} = D\,\nabla^2 C, \tag{60} \end{equation}\]

with \(C(\mathbf{x},t)\) a concentration (or temperature, or probability density) and \(D\) the diffusion coefficient. In one dimension with a point release of total amount \(M\) at the origin, the solution (the Green’s function) is a spreading Gaussian:

\[\begin{equation} C(x,t) = \frac{M}{\sqrt{4\pi D t}}\exp\!\left(-\frac{x^{2}}{4Dt}\right), \qquad \langle x^{2}\rangle = 2Dt. \tag{61} \end{equation}\]

Two features of ((61)) do all the work. First, the profile is Gaussian at every time, so the machinery of Section 1.14 applies directly. Second, the mean squared displacement grows linearly in time, so the typical displacement grows only as \(\sqrt{t}\). Diffusion is efficient over micrometres and hopeless over centimetres, which is the physical reason multicellular organisms require circulatory systems.

Equation (60) governs Brownian motion, drug distribution in tissue, oxygen transport and the Krogh cylinder model, thermal transport in laser–tissue interaction and cryosurgery, and, most directly for imaging, the apparent diffusion coefficient (ADC) of diffusion-weighted MRI, whose reduction within minutes of ischaemic stroke is the basis of early stroke detection (Chapter 4).

Dc <- 1e-9                                     # m^2/s, small molecule in water
x  <- seq(-500, 500, length.out = 1201)*1e-6   # metres
ts <- c(0.05, 0.5, 2, 5)                       # seconds

Cfun <- function(x, t, M = 1) M/sqrt(4*pi*Dc*t)*exp(-x^2/(4*Dc*t))
trapz <- function(x, y) sum(diff(x)*(head(y, -1) + tail(y, -1))/2)

chk <- data.frame(
  time_s = ts,
  mass   = sapply(ts, function(tt) trapz(x, Cfun(x, tt))),
  rms_displacement_um = sqrt(2*Dc*ts)*1e6)
knitr::kable(chk, digits = 4,
             caption = "Mass is conserved; the root-mean-square displacement grows as sqrt(t).")

Table 6: Mass is conserved; the root-mean-square displacement grows as sqrt(t).

time_s	mass	rms_displacement_um
0.05	1	10.0000
0.50	1	31.6228
2.00	1	63.2456
5.00	1	100.0000

cat(sprintf("time to diffuse  10 um : %.3g s\n", (10e-6)^2/(2*Dc)))

## time to diffuse  10 um : 0.05 s

cat(sprintf("time to diffuse   1 cm : %.3g s  (about %.1f hours)\n",
            (1e-2)^2/(2*Dc), (1e-2)^2/(2*Dc)/3600))

## time to diffuse   1 cm : 5e+04 s  (about 13.9 hours)

dd <- do.call(rbind, lapply(ts, function(tt)
  data.frame(x = x*1e6, C = Cfun(x, tt), t = sprintf("t = %.2f s", tt))))
p1 <- ggplot(dd, aes(x, C, colour = t)) +
  geom_line(linewidth = 0.85) +
  scale_colour_manual(values = bpad_pal[1:4]) +
  labs(title = "Diffusive spreading of a point release",
       x = "position (micrometres)", y = "concentration")

tt2 <- seq(0.01, 5, length.out = 200)
p2 <- ggplot(data.frame(t = tt2, msd = 2*Dc*tt2*1e12), aes(t, msd)) +
  geom_line(colour = bpad_pal[2], linewidth = 0.9) +
  labs(title = "Mean squared displacement <x^2> = 2Dt",
       x = "time (s)", y = "MSD (square micrometres)")
bpad_grid(p1, p2, ncol = 2)

Figure 21: Figure 1.21: Solution of the one-dimensional diffusion equation from a point release. Left: the concentration profile spreads as a Gaussian whose width grows as the square root of time while the enclosed mass is conserved. Right: mean squared displacement grows linearly in time.

1.12.3 Other canonical PDEs

The wave equation, \(\partial^2 u/\partial t^2 = c^2\nabla^2 u\), governs acoustic propagation and is the starting point of Chapter 3; the photoacoustic equation is this equation with a heating source term. The Navier–Stokes equations for incompressible flow,

\[\begin{equation} \frac{\partial\mathbf{u}}{\partial t} + (\mathbf{u}\cdot\nabla)\mathbf{u} - \nu\nabla^2\mathbf{u} = -\nabla w + \mathbf{g}, \tag{62} \end{equation}\]

with \(\mathbf{u}\) the flow velocity, \(\nu\) the kinematic viscosity, \(w\) the specific thermodynamic work, and \(\mathbf{g}\) a body force, govern haemodynamics. They are famously unsolved in general: the existence and smoothness of three-dimensional solutions is one of the Clay Mathematics Institute Millennium Prize Problems. Practical haemodynamic modelling proceeds instead by simplification (Poiseuille flow, Windkessel models) or by numerical approximation.

1.12.4 Flux, divergence, curl, and the Laplacian

Let \(\mathbf{F}: V\to\mathbb{R}^{n}\) be a continuously differentiable vector field on a compact region \(V\) with piecewise-smooth boundary \(\partial V\) and outward unit normal \(\mathbf{n}\).

The divergence is the scalar field

\[\begin{equation} \nabla\cdot\mathbf{F} = \sum_{i=1}^{n}\frac{\partial F_i}{\partial x_i}, \tag{63} \end{equation}\]

measuring the local rate at which the field expands; a field with \(\nabla\cdot\mathbf{F} = 0\) is incompressible (blood is treated this way in Chapter 3).

The flux through the boundary measures the net amount transported per unit time:

\[\Phi = \oint_{\partial V}\mathbf{F}\cdot\mathbf{n}\,dS.\]

The curl (in \(\mathbb{R}^3\)) is the vector field

\[\begin{equation} \nabla\times\mathbf{F} = \left(\frac{\partial F_3}{\partial y} - \frac{\partial F_2}{\partial z},\ \ \frac{\partial F_1}{\partial z} - \frac{\partial F_3}{\partial x},\ \ \frac{\partial F_2}{\partial x} - \frac{\partial F_1}{\partial y}\right), \tag{64} \end{equation}\] F_2}{x} - ), \tag{64} \end{equation}

measuring local rotation. The curl is what makes Maxwell’s equations, and hence the RF physics of Chapter 4, vectorial.

The Laplacian is the divergence of the gradient,

\[\begin{equation} \nabla^2 f = \nabla\cdot(\nabla f) = \frac{\partial^2 f}{\partial x^2} + \frac{\partial^2 f}{\partial y^2} + \frac{\partial^2 f}{\partial z^2}. \tag{65} \end{equation}\]

It measures how much the value at a point differs from the average of its neighbours, which is exactly why the discrete Laplacian kernel of Section 1.7.3 detects edges, and why it appears in the diffusion equation: concentration flows toward regions where the local value is below the neighbourhood average.

Worked Example 1.4 (flux through a paraboloid). Compute the flux of \(\mathbf{F} = (1, 2, 4z)\) through the part of the paraboloid \(z = x^2+y^2\) lying above the unit disk \(x^2+y^2\le 1\).

Parametrize the surface as \(\mathbf{r}(r,\theta) = (r\cos\theta,\ r\sin\theta,\ r^2)\), so that

\[\mathbf{r}_r\times\mathbf{r}_\theta = (-2r^2\cos\theta,\ -2r^2\sin\theta,\ r),\]

and on the surface \(\mathbf{F} = (1,\,2,\,4r^2)\). Then

\[\Phi = \int_0^{2\pi}\!\!\int_0^{1}\left(-2r^2\cos\theta - 4r^2\sin\theta + 4r^3\right)dr\,d\theta.\]

The first two terms integrate to zero over a full revolution, leaving

\[\Phi = 2\pi\int_0^1 4r^3\,dr = 2\pi.\]

1.12.5 The divergence theorem

Theorem 1.6 (Gauss–Ostrogradsky). For a continuously differentiable field \(\mathbf{F}\) on a compact region \(V\) with piecewise-smooth boundary,

\[\begin{equation} \int_{V}\big(\nabla\cdot\mathbf{F}\big)\,dV = \oint_{\partial V}\mathbf{F}\cdot\mathbf{n}\,dS. \tag{66} \end{equation}\]

This is the higher-dimensional analogue of the fundamental theorem of calculus ((9)): what accumulates inside a region equals what crosses its boundary. Physically it is a conservation statement, and in that guise it generates the continuity equation, Fick’s laws, and the tracer mass-balance equations of Chapter 6.

Worked Example 1.5 (verifying the divergence theorem). Let \(V\) be the solid cone \(\{x^2+y^2\le z^2,\ 0\le z\le 1\}\) with height and base radius both equal to 1, and let \(\mathbf{F} = (x - y,\ x + z,\ z - y)\).

Left-hand side. The divergence is constant,

\[\nabla\cdot\mathbf{F} = \frac{\partial(x-y)}{\partial x} + \frac{\partial(x+z)}{\partial y} + \frac{\partial(z-y)}{\partial z} = 1 + 0 + 1 = 2,\]

and the cone’s volume is \(\tfrac{1}{3}\pi r^2 h = \pi/3\), so

\[\int_V(\nabla\cdot\mathbf{F})\,dV = 2\cdot\frac{\pi}{3} = \frac{2\pi}{3}.\]

Right-hand side. The boundary consists of the top disk \(S_d\) at \(z=1\) and the lateral cone surface \(S_c\). On \(S_d\) the outward normal is \(+\hat{z}\), so \(\mathbf{F}\cdot\mathbf{n} = z - y = 1 - y\), and

\[\oint_{S_d}\mathbf{F}\cdot\mathbf{n}\,dS = \int_{S_d}(1-y)\,dA = \pi\cdot 1^2 - 0 = \pi.\]

On \(S_c\), parametrized by \((\rho\cos\theta, \rho\sin\theta, \rho)\) with outward normal proportional to \((\rho\cos\theta, \rho\sin\theta, -\rho)\), the integrand reduces to \(-\rho^2\sin^2\theta + 2\rho^2\sin\theta\), giving

\[\oint_{S_c}\mathbf{F}\cdot\mathbf{n}\,dS = \int_0^{2\pi}\!\!\int_0^1\left(-\rho^2\sin^2\theta + 2\rho^2\sin\theta\right)d\rho\,d\theta = -\frac{\pi}{3}.\]

Summing, \(\pi - \pi/3 = 2\pi/3\), which matches the volume integral. \(\checkmark\)

## Monte Carlo confirmation of Worked Example 1.5
set.seed(5); n <- 2e6
xx <- runif(n, -1, 1); yy <- runif(n, -1, 1); zz <- runif(n, 0, 1)
inside  <- sqrt(xx^2 + yy^2) <= zz
vol_mc  <- mean(inside) * (2 * 2 * 1)          # box volume x hit fraction

cat(sprintf("cone volume: Monte Carlo %.4f, exact pi/3 = %.4f\n", vol_mc, pi/3))

## cone volume: Monte Carlo 1.0471, exact pi/3 = 1.0472

cat(sprintf("volume integral of div F = 2V: %.4f, exact 2*pi/3 = %.4f\n",
            2*vol_mc, 2*pi/3))

## volume integral of div F = 2V: 2.0941, exact 2*pi/3 = 2.0944

cat(sprintf("surface flux (analytic): pi - pi/3 = %.4f\n", pi - pi/3))

## surface flux (analytic): pi - pi/3 = 2.0944

Section 1.12 summary.

• Partial derivatives isolate one variable; the gradient assembles them into a direction of steepest ascent and the Hessian classifies critical points.
• The diffusion equation \(\partial C/\partial t = D\nabla^2 C\) has a spreading-Gaussian solution with \(\langle x^2\rangle = 2Dt\): diffusion is fast over micrometres, hopeless over centimetres.
• Divergence measures expansion, curl measures rotation, and the Laplacian measures departure from the local average, which is why it detects edges and drives diffusion.
• The divergence theorem equates interior accumulation with boundary flux: the multivariable fundamental theorem of calculus and the source of every conservation law used later.

Checkpoint 1.12. State the diffusion equation and identify the physical meaning of each term. Then explain, using \(\langle x^2\rangle = 2Dt\), why an ischaemic stroke reduces the apparent diffusion coefficient within minutes while conventional \(T_2\) contrast takes hours to change.

1.13 Ordinary Differential Equations in Biological Systems

Ordinary differential equations describe how a quantity evolves in time. They are the modelling language of pharmacokinetics, tracer kinetics, tumour growth, epidemiology, enzyme kinetics, and neuronal excitability. This section develops the cases that recur in the rest of the textbook.

1.13.1 First-order linear ODEs

The general first-order linear equation

\[\frac{dy}{dt} + p(t)\,y = q(t)\]

is solved by the integrating factor \(\mu(t) = \exp\!\left(\int p(t)\,dt\right)\):

\[\begin{equation} y(t) = \frac{1}{\mu(t)}\left[\int \mu(t)\,q(t)\,dt + C\right]. \tag{67} \end{equation}\]

For the homogeneous constant-coefficient case \(dy/dt = -ky\), the solution is pure exponential decay,

\[\begin{equation} y(t) = y_0 e^{-kt}, \tag{68} \end{equation}\]

with time constant \(\tau = 1/k\) (the time to fall to \(1/e\approx 37\%\)) and half-life \(t_{1/2} = \ln 2/k \approx 0.693/k\). If \(k<0\), ((68)) is exponential growth.

Equation (68) is arguably the single most reused formula in this textbook. It is the radioactive decay law of Chapter 6, the Beer–Lambert attenuation law of Chapters 2 and 5 (with depth replacing time), the \(T_1\) and \(T_2\) relaxation laws of Chapter 4, the fluorescence-lifetime law of Chapter 2, and the first-order drug elimination law below.

1.13.2 Exponential decay: radioactive tracers

\[\frac{dN}{dt} = -\lambda N, \qquad N(t) = N_0 e^{-\lambda t}, \qquad \lambda = \frac{\ln 2}{t_{1/2}}.\]

The choice of radionuclide in Chapter 6 is a direct application: \(^{18}\)F (\(t_{1/2}\approx 110\) min) must survive synthesis, transport, uptake, and a 20-minute acquisition; \(^{99m}\)Tc (\(t_{1/2}\approx 6\) h) matches a clinical workflow measured in hours; \(^{11}\)C (\(t_{1/2}\approx 20\) min) requires an on-site cyclotron.

1.13.3 Growth: exponential, logistic, and Gompertz

Unconstrained growth \(dN/dt = kN\) gives \(N(t) = N_0e^{kt}\) with doubling time \(t_d = \ln 2/k\). Real tumours encounter nutrient, oxygen, and space limits, which the logistic model captures:

\[\begin{equation} \frac{dN}{dt} = kN\left(1 - \frac{N}{K}\right), \qquad N(t) = \frac{K}{1 + \left(\frac{K}{N_0}-1\right)e^{-kt}}, \tag{69} \end{equation}\]

with carrying capacity \(K\) and inflection at \(N = K/2\) (Section 1.1.7). The Gompertz model,

\[\begin{equation} \frac{dN}{dt} = kN\ln\!\frac{K}{N}, \qquad N(t) = K\exp\!\left[\ln\!\left(\frac{N_0}{K}\right)e^{-kt}\right], \tag{70} \end{equation}\]

decelerates earlier, with inflection at \(N = K/e \approx 0.37K\), and empirically fits many measured tumour-volume trajectories better than the logistic model. All three are nonlinear first-order ODEs differing only in how the per-capita growth rate declines with size.

N0 <- 1e6; k <- 0.1; K <- 1e10
kG <- k/log(K/N0)            # matched initial per-capita rate for the Gompertz model
tg <- seq(0, 250, by = 0.5)
gm <- rbind(
  data.frame(t = tg, N = N0*exp(k*tg),                          model = "exponential"),
  data.frame(t = tg, N = K/(1 + (K/N0 - 1)*exp(-k*tg)),         model = "logistic"),
  data.frame(t = tg, N = K*exp(log(N0/K)*exp(-kG*tg)),          model = "Gompertz"))

cat(sprintf("doubling time (exponential): %.2f days\n", log(2)/k))

## doubling time (exponential): 6.93 days

cat(sprintf("logistic inflection: t = %.1f d at N = K/2 = %.3g cells\n",
            log(K/N0 - 1)/k, K/2))

## logistic inflection: t = 92.1 d at N = K/2 = 5e+09 cells

cat(sprintf("Gompertz inflection: t = %.1f d at N = K/e = %.3g cells\n",
            log(log(K/N0))/kG, K/exp(1)))

## Gompertz inflection: t = 204.5 d at N = K/e = 3.68e+09 cells

cat("\nmodel predictions at selected times (cells):\n")

## 
## model predictions at selected times (cells):

print(signif(data.frame(
  day         = c(5, 20, 50, 100),
  exponential = N0*exp(k*c(5, 20, 50, 100)),
  logistic    = K/(1 + (K/N0 - 1)*exp(-k*c(5, 20, 50, 100))),
  Gompertz    = K*exp(log(N0/K)*exp(-kG*c(5, 20, 50, 100)))), 3))

##   day exponential logistic Gompertz
## 1   5    1.65e+06 1.65e+06 1.63e+06
## 2  20    7.39e+06 7.38e+06 6.04e+06
## 3  50    1.48e+08 1.46e+08 4.74e+07
## 4 100    2.20e+10 6.88e+09 4.46e+08

p1 <- ggplot(gm, aes(t, N, colour = model)) +
  geom_line(linewidth = 0.9) +
  geom_hline(yintercept = K, linetype = "dotted", colour = "grey50") +
  scale_y_log10() +
  coord_cartesian(ylim = c(1e6, 2e10)) +
  scale_colour_manual(values = bpad_pal[1:3]) +
  labs(title = "Tumour-growth models", x = "time (days)", y = "cell count (log scale)")

nuc <- data.frame(nuclide = c("F-18 (110 min)", "Tc-99m (6.0 h)", "C-11 (20 min)"),
                  half_life_min = c(109.77, 360.6, 20.36))
tm  <- seq(0, 600, by = 1)
dec <- do.call(rbind, lapply(seq_len(nrow(nuc)), function(i)
  data.frame(t = tm, A = 100*exp(-log(2)/nuc$half_life_min[i]*tm),
             nuclide = nuc$nuclide[i])))
p2 <- ggplot(dec, aes(t, A, colour = nuclide)) +
  geom_line(linewidth = 0.9) +
  geom_hline(yintercept = 50, linetype = "dotted", colour = "grey50") +
  scale_colour_manual(values = bpad_pal[c(1, 3, 2)]) +
  labs(title = "Radionuclide decay", x = "time (minutes)",
       y = "remaining activity (% of initial)")
bpad_grid(p1, p2, ncol = 2)

Figure 22: Figure 1.22: Left: three tumour-growth models sharing the same initial size and the same initial per-capita growth rate, on a logarithmic axis. They are nearly indistinguishable for the first few weeks, the regime in which most imaging measurements are made, and diverge by orders of magnitude thereafter. Right: decay of three PET and SPECT radionuclides.

1.13.4 Compartmental pharmacokinetics

Pharmacokinetics describes absorption, distribution, metabolism, and excretion. Compartmental models express these with ODEs and yield the quantities that dosing decisions actually use.

One compartment, intravenous bolus. Treating the body as a single well-mixed volume with first-order elimination,

\[\begin{equation} \frac{dC}{dt} = -k_{el}C, \qquad C(t) = \frac{\mathrm{Dose}}{V_d}e^{-k_{el}t}, \tag{71} \end{equation}\]

where \(V_d\) is the volume of distribution (the apparent volume that would contain the whole dose at the measured concentration), \(t_{1/2} = \ln 2/k_{el}\), and the clearance is \(CL = k_{el}V_d\). Clinically, \(V_d\) reports on tissue distribution: a drug confined to plasma such as heparin has \(V_d \approx 3\) L, whereas an extensively tissue-bound drug such as digoxin may have \(V_d\approx 500\) L. Clearance sets the maintenance dose rate needed to hold a target concentration.

Two compartments. Many drugs show biexponential plasma decay: a rapid distribution (\(\alpha\)) phase followed by a slower elimination (\(\beta\)) phase. Separating a central compartment (blood and well-perfused organs) from a peripheral one (muscle, fat),

\[\begin{equation} \frac{dC_1}{dt} = -(k_{10}+k_{12})C_1 + k_{21}C_2, \qquad \frac{dC_2}{dt} = k_{12}C_1 - k_{21}C_2, \tag{72} \end{equation}\]

whose solution in the central compartment is \(C_1(t) = Ae^{-\alpha t} + Be^{-\beta t}\) with \(\alpha>\beta\). The area under the curve is

\[\begin{equation} \mathrm{AUC} = \int_0^{\infty}C_1(t)\,dt = \frac{A}{\alpha} + \frac{B}{\beta} = \frac{\mathrm{Dose}}{CL}, \tag{73} \end{equation}\]

an identity used directly in chemotherapy dosing (carboplatin, for example, is dosed to a target AUC using the Calvert formula).

The system ((72)) is linear, so it is also an eigenvalue problem: writing \(\dot{\mathbf{C}} = A\mathbf{C}\), the macro rate constants \(-\alpha\) and \(-\beta\) are precisely the eigenvalues of \(A\) (Section 1.17). This is the connection between compartmental kinetics and linear algebra, and it generalizes to the multi-tissue tracer models of Chapter 6.

k10 <- 0.2; k12 <- 0.4; k21 <- 0.3      # per hour
Dose <- 100; V1 <- 10                   # mg and litres
C0 <- Dose/V1

## Macro constants as eigenvalues of the system matrix (Section 1.17)
Amat <- matrix(c(-(k10 + k12), k21,
                   k12,       -k21), nrow = 2, byrow = TRUE)
ev <- eigen(Amat)$values
alpha <- -min(ev); beta <- -max(ev)
Acoef <- C0*(alpha - k21)/(alpha - beta)
Bcoef <- C0*(k21 - beta)/(alpha - beta)

AUCinf <- Acoef/alpha + Bcoef/beta
AUC24  <- Acoef/alpha*(1 - exp(-24*alpha)) + Bcoef/beta*(1 - exp(-24*beta))
cat(sprintf("alpha = %.4f /h (t1/2 = %.2f h);  beta = %.4f /h (t1/2 = %.2f h)\n",
            alpha, log(2)/alpha, beta, log(2)/beta))

## alpha = 0.8275 /h (t1/2 = 0.84 h);  beta = 0.0725 /h (t1/2 = 9.56 h)

cat(sprintf("A = %.3f, B = %.3f mg/L  (A + B = C0 = %.1f)\n", Acoef, Bcoef, C0))

## A = 6.987, B = 3.013 mg/L  (A + B = C0 = 10.0)

cat(sprintf("AUC(0-24) = %.2f,  AUC(0-inf) = %.2f mg*h/L\n", AUC24, AUCinf))

## AUC(0-24) = 42.71,  AUC(0-inf) = 50.00 mg*h/L

cat(sprintf("clearance from AUC = Dose/AUC = %.2f L/h;  k10*V1 = %.2f L/h\n",
            Dose/AUCinf, k10*V1))

## clearance from AUC = Dose/AUC = 2.00 L/h;  k10*V1 = 2.00 L/h

times <- seq(0, 24, by = 0.05)
C1 <- Acoef*exp(-alpha*times) + Bcoef*exp(-beta*times)
C2 <- C0*k12/(alpha - beta)*(exp(-beta*times) - exp(-alpha*times))

if (has_pkg("deSolve")) {                        # numerical check
  pk <- function(t, y, p) with(as.list(c(y, p)),
        list(c(-(k10 + k12)*C1 + k21*C2, k12*C1 - k21*C2)))
  num <- deSolve::ode(c(C1 = C0, C2 = 0), times, pk,
                      c(k10 = k10, k12 = k12, k21 = k21))
  cat(sprintf("max |numerical - analytic| = %.2e mg/L\n", max(abs(num[, 2] - C1))))
}

## max |numerical - analytic| = 2.16e-06 mg/L

pk_df <- rbind(data.frame(t = times, C = C1,  cmp = "central C1(t)"),
               data.frame(t = times, C = C2,  cmp = "peripheral C2(t)"),
               data.frame(t = times, C = C0*exp(-k10*times),
                          cmp = "one-compartment (k_el = k10)"))
ggplot(pk_df[pk_df$C > 1e-3, ], aes(t, C, colour = cmp)) +
  geom_line(linewidth = 0.9) +
  scale_y_log10() +
  scale_colour_manual(values = bpad_pal[1:3]) +
  labs(title = "Two-compartment kinetics after an IV bolus",
       x = "time (hours)", y = "concentration (mg/L, log scale)")

Figure 23: Figure 1.23: Two-compartment pharmacokinetics after an intravenous bolus, on a logarithmic concentration axis. The central-compartment curve is biexponential: the steep alpha phase is redistribution into tissue, the shallow beta phase is true elimination. A one-compartment model fitted to late data alone would badly misestimate early exposure.

Physiologically based models. PBPK models replace abstract compartments with anatomical organs connected by measured blood flows,

\[\begin{equation} V_i\frac{dC_i}{dt} = Q_i\left(C_{\text{art}} - C_{\text{ven},i}\right) - f_{\text{met},i} + R_i, \tag{74} \end{equation}\]

where \(V_i\) is the organ volume, \(Q_i\) its perfusion, \(f_{\text{met},i}\) its metabolic loss, and \(R_i\) any binding or reaction term. PBPK models predict tissue concentrations in organs that cannot be sampled and underpin quantitative molecular imaging in Chapter 6.

1.13.5 Nonlinear systems: Lotka–Volterra

Coupled nonlinear ODEs produce behaviour no single-compartment model can: oscillation, bistability, and thresholds. The Lotka–Volterra predator–prey system is the classical example,

\[\begin{equation} \frac{dx}{dt} = a\,x - b\,xy, \qquad \frac{dy}{dt} = -c\,y + d\,xy, \tag{75} \end{equation}\]

with prey \(x\) and predator \(y\). Each population grows or decays exponentially on its own; the interaction terms \(\pm xy\) couple them. The nontrivial equilibrium is \((x^{*}, y^{*}) = (c/d,\ a/b)\), around which trajectories form closed orbits in the phase plane: populations cycle indefinitely rather than settling.

The structure recurs well beyond ecology. The same competitive and predatory couplings model tumour–immune interaction, host–pathogen dynamics, and competing bacterial strains under antibiotic pressure; the same phase-plane analysis identifies the thresholds at which one population is driven to extinction.

## Fourth-order Runge-Kutta in base R (no package dependency)
rk4 <- function(f, y0, times, pars) {
  Y <- matrix(NA_real_, length(times), length(y0)); Y[1, ] <- y0
  for (i in 2:length(times)) {
    h <- times[i] - times[i-1]; y <- Y[i-1, ]; t <- times[i-1]
    k1 <- f(t, y, pars);            k2 <- f(t + h/2, y + h/2*k1, pars)
    k3 <- f(t + h/2, y + h/2*k2, pars); k4 <- f(t + h, y + h*k3, pars)
    Y[i, ] <- y + h/6*(k1 + 2*k2 + 2*k3 + k4)
  }
  Y
}

lv <- function(t, y, p) c(p["a"]*y[1] - p["b"]*y[1]*y[2],
                         -p["c"]*y[2] + p["d"]*y[1]*y[2])
pars <- c(a = 0.8, b = 0.04, c = 0.5, d = 0.02)
tv   <- seq(0, 60, by = 0.02)
Y    <- rk4(lv, c(20, 10), tv, pars)

cat(sprintf("equilibrium: prey* = c/d = %.1f, predator* = a/b = %.1f\n",
            pars["c"]/pars["d"], pars["a"]/pars["b"]))

## equilibrium: prey* = c/d = 25.0, predator* = a/b = 20.0

ts_df <- rbind(data.frame(t = tv, n = Y[, 1], sp = "prey"),
               data.frame(t = tv, n = Y[, 2], sp = "predator"))
p1 <- ggplot(ts_df, aes(t, n, colour = sp)) +
  geom_line(linewidth = 0.85) +
  scale_colour_manual(values = bpad_pal[1:2]) +
  labs(title = "Population time courses", x = "time", y = "population size")

p2 <- ggplot(data.frame(x = Y[, 1], y = Y[, 2]), aes(x, y)) +
  geom_path(colour = bpad_pal[3], linewidth = 0.8) +
  geom_point(data = data.frame(x = pars["c"]/pars["d"], y = pars["a"]/pars["b"]),
             aes(x, y), colour = bpad_pal[2], size = 3) +
  geom_point(data = data.frame(x = 20, y = 10), aes(x, y),
             colour = "black", size = 2.4, shape = 17) +
  labs(title = "Phase plane (equilibrium: circle; start: triangle)",
       x = "prey", y = "predator")
bpad_grid(p1, p2, ncol = 2)

Figure 24: Figure 1.24: The Lotka-Volterra system. Left: prey and predator populations oscillate with the predator peak lagging the prey peak. Right: the phase-plane trajectory is a closed orbit encircling the equilibrium, so the system neither converges to it nor diverges.

1.13.6 Higher-order equations and state-space form

Several biophysical systems require second-order or coupled equations:

• Harmonic oscillator \(m\ddot{x} + kx = 0\): mass–spring systems and molecular vibrations (Chapter 2, Raman).
• Damped oscillator \(m\ddot{x} + c\dot{x} + kx = 0\): underdamped when \(c^2 < 4mk\), critically damped when \(c^2 = 4mk\), overdamped when \(c^2 > 4mk\). The three regimes correspond to the three cases of the characteristic quadratic, and the underdamped case has complex roots, which is where the phasors of Section 1.6 re-enter.
• Driven oscillator \(m\ddot{x} + c\dot{x} + kx = F_0\cos\omega t\): resonance near \(\omega \approx \sqrt{k/m}\), the mechanism of both magnetic resonance (Chapter 4) and acoustic transducer design (Chapter 3).
• Hodgkin–Huxley model: four coupled nonlinear ODEs generating the neuronal action potential.

Any \(n\)-th order ODE can be rewritten as \(n\) coupled first-order equations (state-space form), which is what allows a single numerical solver (Euler, Runge–Kutta, or the deSolve package) to handle all of them.

From equation to oncology. Exponential, logistic, and Gompertz models are indistinguishable when a tumour is small, which is exactly the regime in which most imaging measurements are made. Distinguishing them requires either long longitudinal follow-up or an independent estimate of the carrying capacity. This is why Chapter 8 treats growth-rate estimation as a statistical problem with quantified uncertainty rather than a curve-fitting exercise: the model choice, not just the parameter, carries the clinical implication.

Section 1.13 summary.

• First-order linear ODEs have exponential solutions; the same \(y_0e^{-kt}\) describes radioactive decay, attenuation, relaxation, and drug clearance.
• Logistic and Gompertz models add saturation, with inflection at \(K/2\) and \(K/e\) respectively.
• Compartmental models convert dosing questions into linear systems; macro rate constants are eigenvalues, and \(\mathrm{AUC} = \mathrm{Dose}/CL\).
• Nonlinear coupled systems (Lotka–Volterra) generate oscillation and thresholds that no single compartment can.
• Any higher-order ODE reduces to a first-order system, which is what makes one numerical solver sufficient.

Checkpoint 1.13. Show that \(dN/dt = kN\) has solution \(N(t) = N_0e^{kt}\), and derive the doubling time. Then explain why a tumour whose volume doubles every 60 days is not necessarily growing exponentially, and what additional measurement would settle the question.

1.14 Probability and Statistics

Classical mechanics is deterministic, but biomedical measurement is not. Photon emission, thermal agitation, biological variability between patients, and measurement error are all irreducibly stochastic. Probability is the calculus of that uncertainty, and statistics is the machinery for reasoning backwards from data to the process that generated it. Additional depth is available in the SOCR Probability and Statistics EBook and in the Data Science and Predictive Analytics textbook.

1.14.1 The Kolmogorov axioms

A probability is a function mapping events to \([0,1]\subset\mathbb{R}\) satisfying the 1933 Kolmogorov axioms:

• Axiom 1. \(0 \le P(E) \le 1\) for every event \(E\).
• Axiom 2. \(P(\Omega) = 1\), where \(\Omega\) is the entire sample space.
• Axiom 3. For any countable sequence of mutually exclusive events, \(P\!\left(\bigcup_{i=1}^{\infty}E_i\right) = \sum_{i=1}^{\infty}P(E_i)\).

Everything else, including conditional probability, Bayes’ theorem, and the entire theory of estimation, follows from these three statements plus definitions.

1.14.2 Conditional probability and Bayes’ theorem

\[\begin{equation} P(A\mid B) = \frac{P(A\cap B)}{P(B)}, \quad P(B)>0. \tag{76} \end{equation}\]

Events are independent when \(P(A\cap B) = P(A)P(B)\), equivalently \(P(A\mid B) = P(A)\).

Theorem 1.7 (Bayes). For a hypothesis \(H\) and evidence \(E\),

\[\begin{equation} P(H\mid E) = \frac{P(E\mid H)\,P(H)}{P(E)} = \frac{P(E\mid H)P(H)}{\sum_i P(E\mid H_i)P(H_i)}. \tag{77} \end{equation}\]

Bayes’ theorem is the formal statement that a test result must be interpreted in light of the prior probability of disease. Sensitivity, \(P(+\mid D)\), and specificity, \(P(-\mid \bar{D})\), are properties of the test. Positive predictive value, \(P(D\mid +)\), is not: it depends on prevalence, and it is what the patient and clinician actually need.

Worked Example 1.6 (why screening a low-prevalence population is hard). A test has sensitivity 95% and specificity 90%. In a population with 1% prevalence,

\[\mathrm{PPV} = P(D\mid+) = \frac{0.95\times 0.01}{0.95\times 0.01 + 0.10\times 0.99} = \frac{0.0095}{0.1085} = 0.088.\]

Fewer than one in eleven positive results is a true positive, despite a test that sounds excellent. The reason is that the 10% false-positive rate acts on the 99% of the population who are healthy, generating far more false positives than the 95% sensitivity generates true positives among the 1% who are ill.

Raising prevalence to 10% (for example, by testing a symptomatic rather than a general population) gives \(\mathrm{PPV} = 0.095/0.185 = 0.514\), a nearly six-fold improvement obtained without changing the test at all. This is the quantitative argument for targeted rather than universal screening, and it recurs in Chapter 8 as the calibration problem in predictive modelling.

sens <- 0.95; spec <- 0.90
prev <- 10^seq(-3, log10(0.9), length.out = 400)
ppv  <- sens*prev/(sens*prev + (1 - spec)*(1 - prev))
npv  <- spec*(1 - prev)/(spec*(1 - prev) + (1 - sens)*prev)

knitr::kable(data.frame(
  prevalence = c(0.001, 0.01, 0.05, 0.10, 0.50),
  PPV = round(sens*c(0.001,0.01,0.05,0.10,0.50) /
        (sens*c(0.001,0.01,0.05,0.10,0.50) +
         (1-spec)*(1-c(0.001,0.01,0.05,0.10,0.50))), 4),
  NPV = round(spec*(1-c(0.001,0.01,0.05,0.10,0.50)) /
        (spec*(1-c(0.001,0.01,0.05,0.10,0.50)) +
         (1-sens)*c(0.001,0.01,0.05,0.10,0.50)), 5)),
  caption = "Predictive values at five prevalences (sensitivity 0.95, specificity 0.90).")

Table 7: Predictive values at five prevalences (sensitivity 0.95, specificity 0.90).

prevalence	PPV	NPV
0.001	0.0094	0.99994
0.010	0.0876	0.99944
0.050	0.3333	0.99708
0.100	0.5135	0.99387
0.500	0.9048	0.94737

bd <- rbind(data.frame(prev = prev, v = ppv, q = "positive predictive value"),
            data.frame(prev = prev, v = npv, q = "negative predictive value"))
ggplot(bd, aes(prev, v, colour = q)) +
  geom_line(linewidth = 0.95) +
  scale_x_log10() +
  scale_colour_manual(values = bpad_pal[1:2]) +
  labs(title = "Predictive value depends on prevalence, not only on the test",
       x = "disease prevalence (log scale)", y = "probability")

Figure 25: Figure 1.25: Positive and negative predictive value as functions of disease prevalence for a test with 95% sensitivity and 90% specificity. Sensitivity and specificity are fixed properties of the test; predictive values are not, and at low prevalence the PPV collapses.

Bayesian inference treats a parameter \(\theta\) as a random variable with prior \(P(\theta)\), updated by data \(D\) to the posterior \(P(\theta\mid D) \propto P(D\mid\theta)P(\theta)\). This is the framework behind maximum a posteriori reconstruction in PET (Chapter 6) and behind regularized estimation in Chapter 8, where a penalty term is exactly a log-prior.

1.14.3 Statistical moments

Four summary quantities characterize the shape of a distribution.

The mean locates its centre:

\[\begin{equation} \bar{x} = \frac{1}{N}\sum_i f_i x_i = \sum_i x_i p_i, \qquad p_i = f_i/N. \tag{78} \end{equation}\]

The standard deviation measures dispersion in the units of the data:

\[\begin{equation} s = \sqrt{\frac{1}{N-1}\sum_{i=1}^{N}\left(x_i - \bar{x}\right)^{2}}. \tag{79} \end{equation}\]

Skewness measures asymmetry:

\[\begin{equation} \mathrm{Skew} = \frac{m_3}{m_2^{3/2}}, \qquad m_r = \frac{1}{N}\sum_{i=1}^{N}(x_i - \bar{x})^{r}. \tag{80} \end{equation}\]

A positive skew means a long right tail (most data on the left); a negative skew means the opposite.

Kurtosis measures tail weight:

\[\begin{equation} \mathrm{Kurt} = \frac{m_4}{m_2^{2}}, \qquad \text{excess kurtosis} = \frac{m_4}{m_2^2} - 3, \tag{81} \end{equation}\]

with the subtraction of 3 making the excess kurtosis zero for a normal distribution. Positive excess kurtosis indicates heavier tails and therefore more frequent outliers.

These are not merely descriptive. A PET or SPECT intensity histogram is strongly right-skewed with heavy tails, because most voxels are background and a few are intensely avid. That shape is the reason the mean is a poor summary of tracer uptake, why SUV\(_{\max}\) and SUV\(_{\mathrm{peak}}\) are reported instead, and why Gaussian-based statistics must be applied cautiously to raw emission data (Chapter 6).

moments <- function(x) {
  n <- length(x); m <- mean(x)
  m2 <- mean((x - m)^2); m3 <- mean((x - m)^3); m4 <- mean((x - m)^4)
  c(mean = m, sd = sqrt(m2*n/(n - 1)),
    skewness = m3/m2^1.5, excess_kurtosis = m4/m2^2 - 3)
}

set.seed(7)
pet_like <- c(rgamma(20000, shape = 2, scale = 8),  # avid tissue
              rep(0, 12000) + abs(rnorm(12000, 0, 1.5)))  # background
round(moments(pet_like), 4)

##            mean              sd        skewness excess_kurtosis 
##         10.4384         11.4688          1.5845          3.1604

ggplot(data.frame(v = pet_like), aes(v)) +
  geom_histogram(bins = 80, fill = bpad_pal[1], colour = "white", linewidth = 0.15) +
  geom_vline(xintercept = mean(pet_like), colour = bpad_pal[2],
             linetype = "dashed", linewidth = 0.9) +
  geom_vline(xintercept = median(pet_like), colour = bpad_pal[3],
             linetype = "dotted", linewidth = 0.9) +
  annotate("text", x = mean(pet_like) + 6, y = Inf, vjust = 2, hjust = 0, size = 3.2,
           label = sprintf("mean = %.1f (dashed)\nmedian = %.1f (dotted)",
                           mean(pet_like), median(pet_like))) +
  labs(title = "Right-skewed intensity distribution of a synthetic emission image",
       x = "voxel intensity (arbitrary units)", y = "count")

Figure 26: Figure 1.26: A synthetic emission-image intensity histogram. The distribution is strongly right-skewed with a large background peak, so the mean sits well above the mode and well below the intense voxels that carry the diagnostic signal.

The real-data version of this computation, using an FDG PET volume from the SOCR archive, is guarded by run_network.

library(brainR); library(httr)
httr::set_config(httr::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(httr::GET(url, httr::write_disk(dest, overwrite = TRUE)))
pet <- oro.nifti::readNIfTI(dest, reorient = FALSE)

slice_vec <- c(pet[, , 40])
round(moments(slice_vec), 4)
image(t(pet[, , 40]), col = gray.colors(256, 0, 1), asp = 1, axes = FALSE)
hist(slice_vec, breaks = 80, main = "FDG PET slice intensity histogram",
     xlab = "intensity")

1.14.4 Common distributions

Binomial. For \(n\) independent trials each succeeding with probability \(p\),

\[\begin{equation} P(X = k) = \binom{n}{k}p^{k}(1-p)^{n-k}, \qquad \mathbb{E}[X] = np, \quad \mathrm{Var}(X) = np(1-p). \tag{82} \end{equation}\]

It models any counted proportion: responders in a trial arm, detected events among emitted ones, correctly classified cases in Chapter 8.

Normal (Gaussian).

\[\begin{equation} f(x\mid\mu,\sigma) = \frac{1}{\sigma\sqrt{2\pi}}\exp\!\left(-\frac{(x-\mu)^{2}}{2\sigma^{2}}\right), \tag{83} \end{equation}\]

with 68.27%, 95.45%, and 99.73% of the mass within one, two, and three standard deviations. Standardizing by \(Z = (X-\mu)/\sigma\) maps any normal to the standard normal. Its ubiquity is explained by the central limit theorem: a sum of many independent finite-variance contributions tends to a normal distribution regardless of the individual distributions, which is precisely the situation for thermal noise in an MR receiver or an ultrasound amplifier.

Poisson. For events occurring independently at a known mean rate \(\lambda\) per interval,

\[\begin{equation} P(X = k) = \frac{\lambda^{k}e^{-\lambda}}{k!},\quad k\in\mathbb{N}_0, \qquad \mathbb{E}[X] = \mathrm{Var}(X) = \lambda. \tag{84} \end{equation}\]

The equality of mean and variance is the defining property, and it is the origin of the \(\sqrt{N}\) noise law of Section 1.15. The Poisson arises as the limit of the binomial when \(n\to\infty\), \(p\to 0\) with \(np = \lambda\) fixed, which is exactly the physical situation of radioactive decay: an enormous number of nuclei, each with a tiny probability of decaying in the observation window. For \(\lambda \gtrsim 20\) the normal approximation \(\mathcal{N}(\lambda,\lambda)\) is usually adequate.

kb <- 50:90
p1 <- ggplot(data.frame(k = kb, p = dbinom(kb, 100, 0.7)), aes(k, p)) +
  geom_col(fill = bpad_pal[1], width = 0.75) +
  labs(title = "Binomial(n = 100, p = 0.7)", x = "successes k", y = "P(X = k)")

pd <- do.call(rbind, lapply(c(2, 10, 40), function(l)
  data.frame(k = 0:70, p = dpois(0:70, l), lam = sprintf("lambda = %d", l))))
p2 <- ggplot(pd, aes(k, p, colour = lam)) +
  geom_line(linewidth = 0.85) + geom_point(size = 0.7) +
  scale_colour_manual(values = bpad_pal[1:3]) +
  labs(title = "Poisson at three rates", x = "counts k", y = "P(X = k)")

kk <- 1500:2500
p3 <- ggplot(data.frame(k = kk, pois = dpois(kk, 2000),
                        norm = dnorm(kk, 2000, sqrt(2000))), aes(k)) +
  geom_line(aes(y = pois, colour = "Poisson(2000)"), linewidth = 1.1) +
  geom_line(aes(y = norm, colour = "Normal(2000, sqrt(2000))"),
            linewidth = 0.6, linetype = "dashed") +
  scale_colour_manual(values = bpad_pal[c(1, 2)]) +
  labs(title = "Poisson to normal at large lambda", x = "counts k", y = "density")

bpad_grid(p1, p2, p3, ncol = 3)

Figure 27: Figure 1.27: Three workhorse distributions. Left: the binomial for 100 trials at p = 0.7. Centre: the Poisson at three rates, converging in shape toward the normal as lambda grows. Right: the Poisson at lambda = 2 is markedly skewed, while at lambda = 2000 it is visually indistinguishable from the normal.

## A PET-detector question: P(15 counts in one minute) when the rate is 10 per 5 min
cat(sprintf("lambda = 2 counts/min;  P(X = 15) = %.3e\n", dpois(15, 2)))

## lambda = 2 counts/min;  P(X = 15) = 3.391e-09

1.14.5 Joint distributions and the multivariate normal

Real measurements are rarely univariate: a radiomic feature vector, a multi-echo MR signal, and a multi-wavelength optical spectrum are all vector-valued. The multivariate normal for \(\mathbf{x}\in\mathbb{R}^{d}\) is

\[\begin{equation} f(\mathbf{x}) = \frac{1}{(2\pi)^{d/2}\lvert\Sigma\rvert^{1/2}} \exp\!\left(-\tfrac{1}{2}(\mathbf{x}-\boldsymbol{\mu})^{T}\Sigma^{-1}(\mathbf{x}-\boldsymbol{\mu})\right), \tag{85} \end{equation}\]

with mean vector \(\boldsymbol{\mu}\) and covariance matrix \(\Sigma\). The quadratic form in the exponent is the squared Mahalanobis distance, the natural distance in a correlated space.

Three uses follow immediately: Gaussian mixture models for image segmentation (Chapter 7), Kalman filtering for motion tracking, and principal component analysis, which is nothing other than the eigendecomposition of \(\Sigma\) (Sections 1.17 and 1.18).

A random variable is discrete when it takes values in a countable set (an AP examination score, a photon count, a tumour stage) and continuous when it can take any value in an interval (a latitude, a concentration, a relaxation time). The distinction determines whether probabilities are summed over a mass function or integrated over a density.

Section 1.14 summary.

• Three Kolmogorov axioms generate the whole of probability theory.
• Bayes’ theorem converts test characteristics into patient-specific probabilities; predictive value depends on prevalence, sensitivity and specificity do not.
• The first four moments describe centre, spread, asymmetry, and tail weight; emission-imaging histograms are strongly right-skewed with heavy tails.
• Binomial for counted proportions, Poisson for counted events (mean = variance), normal for sums of many small effects (central limit theorem).
• The multivariate normal, parameterized by \(\boldsymbol{\mu}\) and \(\Sigma\), is the basis of segmentation, tracking, and PCA.

Checkpoint 1.14. A test has 99% sensitivity and 99% specificity, and disease prevalence is 1 in 10,000. Compute the PPV. Would you deploy this test for population screening? What would you change?

1.15 Noise Models, SNR, and CNR

Every measurement is corrupted. Understanding which noise distribution applies determines the correct estimator, the correct likelihood, the correct loss function in Chapter 8, and the correct answer to the practical question of how much longer to scan.

1.15.1 Gaussian (thermal) noise

Additive Gaussian noise arises when many independent small disturbances sum, by the central limit theorem:

\[\begin{equation} y = x + n, \qquad n \sim \mathcal{N}(0, \sigma^{2}). \tag{86} \end{equation}\]

Its defining feature is that \(\sigma\) is independent of the signal. Sources include Johnson–Nyquist thermal noise in MR receiver coils and electronic circuits, amplifier noise in ultrasound front-ends, and background fluctuation in optical detection. Under Gaussian noise the maximum-likelihood estimator is the least-squares estimator, which is why squared-error loss is the default for regression.

1.15.2 Poisson (quantum) noise

Poisson noise, also called shot or quantum noise, arises from counting discrete independent events. If the expected count is \(\lambda\), then \(N\sim\mathrm{Poisson}(\lambda)\) with \(\mathbb{E}[N] = \mathrm{Var}(N) = \lambda\), so

\[\begin{equation} \sigma = \sqrt{\lambda}, \qquad \mathrm{SNR} = \frac{\lambda}{\sqrt{\lambda}} = \sqrt{\lambda}. \tag{87} \end{equation}\]

Noise is therefore signal-dependent: bright regions are absolutely noisier but relatively cleaner. The \(\sqrt{N}\) law in ((87)) is the central practical fact of photon-limited imaging: to double the SNR you must quadruple the counts, hence quadruple the dose, the scan time, or the detector efficiency.

Poisson statistics dominate X-ray CT (detector photon flux), PET and SPECT (decay and detection are both Poisson processes), fluorescence microscopy, and OCT.

1.15.3 Rician noise

In coherent imaging (MRI, ultrasound, radar) the detected signal is a complex phasor \(S = S_{\mathrm{re}} + iS_{\mathrm{im}}\) with independent Gaussian noise in each channel:

\[S_{\mathrm{re}}\sim\mathcal{N}(A\cos\theta, \sigma^{2}), \qquad S_{\mathrm{im}}\sim\mathcal{N}(A\sin\theta,\sigma^{2}).\]

Forming a magnitude image is a nonlinear operation, \(M = \lvert S\rvert = \sqrt{S_{\mathrm{re}}^2 + S_{\mathrm{im}}^2}\), and the result follows a Rice distribution:

\[\begin{equation} p(M \mid A,\sigma) = \frac{M}{\sigma^{2}}\exp\!\left(-\frac{M^{2}+A^{2}}{2\sigma^{2}}\right) I_0\!\left(\frac{MA}{\sigma^{2}}\right), \tag{88} \end{equation}\]

with \(I_0\) the modified Bessel function of the first kind. Its behaviour has two regimes:

• High SNR (\(A/\sigma \gg 1\)): approximately \(\mathcal{N}(A,\sigma^2)\), and ordinary Gaussian methods apply.
• Low SNR (\(A/\sigma \ll 1\)): the distribution becomes Rayleigh, and the magnitude is biased upward: even when the true signal is zero, \(\mathbb{E}[M] = \sigma\sqrt{\pi/2}\approx 1.253\,\sigma > 0\).

This positive bias has real consequences. In diffusion MRI at high \(b\)-values and in fMRI at low signal, magnitude images systematically overestimate intensity, biasing the apparent diffusion coefficient and any quantitative parameter derived from them, which is why Chapter 4 discusses noise-floor correction rather than treating magnitude data as if it were Gaussian.

set.seed(42); sigma <- 5; nsamp <- 3e4

gauss <- 25 + rnorm(nsamp, 0, sigma)
pois  <- lapply(c(5, 20, 100), function(l) rpois(nsamp, l))
rice  <- lapply(c(0, 5, 30), function(A)
  sqrt(rnorm(nsamp, A, sigma)^2 + rnorm(nsamp, 0, sigma)^2))

cat("Poisson: sample SD vs sqrt(lambda):\n")

## Poisson: sample SD vs sqrt(lambda):

print(round(rbind(lambda = c(5, 20, 100),
                  sample_sd = sapply(pois, sd),
                  sqrt_lambda = sqrt(c(5, 20, 100))), 3))

##              [,1]   [,2]    [,3]
## lambda      5.000 20.000 100.000
## sample_sd   2.236  4.473  10.061
## sqrt_lambda 2.236  4.472  10.000

cat(sprintf("\nRician at A = 0: mean magnitude %.3f (theory sigma*sqrt(pi/2) = %.3f)\n",
            mean(rice[[1]]), sigma*sqrt(pi/2)))

## 
## Rician at A = 0: mean magnitude 6.251 (theory sigma*sqrt(pi/2) = 6.267)

cat(sprintf("Rician at A = 30: mean magnitude %.3f (true amplitude 30)\n",
            mean(rice[[3]])))

## Rician at A = 30: mean magnitude 30.392 (true amplitude 30)

pg <- ggplot(data.frame(v = gauss), aes(v)) +
  geom_histogram(bins = 60, fill = bpad_pal[1], colour = "white", linewidth = 0.1) +
  geom_vline(xintercept = 25, colour = bpad_pal[2], linetype = "dashed") +
  labs(title = "Gaussian: sigma independent of signal", x = "measured value", y = "count")

pp <- ggplot(do.call(rbind, lapply(seq_along(pois), function(i)
  data.frame(v = pois[[i]], lam = sprintf("lambda = %d", c(5,20,100)[i])))),
  aes(v, fill = lam)) +
  geom_histogram(bins = 60, position = "identity", alpha = 0.55) +
  scale_fill_manual(values = bpad_pal[1:3]) +
  labs(title = "Poisson: sigma = sqrt(lambda)", x = "counts", y = "count")

pr <- ggplot(do.call(rbind, lapply(seq_along(rice), function(i)
  data.frame(v = rice[[i]], A = sprintf("A = %d", c(0,5,30)[i])))),
  aes(v, fill = A)) +
  geom_histogram(bins = 60, position = "identity", alpha = 0.55) +
  scale_fill_manual(values = bpad_pal[1:3]) +
  labs(title = "Rician: biased high at low SNR", x = "magnitude", y = "count")

bpad_grid(pg, pp, pr, ncol = 3)

Figure 28: Figure 1.28: Three noise models. Left: Gaussian noise is signal-independent. Centre: Poisson noise widens as the square root of the mean, so the relative noise falls. Right: the Rice distribution at three signal levels; at A = 0 the mean magnitude is sigma*sqrt(pi/2), not zero, the low-SNR bias of magnitude images.

1.15.4 SNR, CNR, and the square-root law

\[\begin{equation} \mathrm{SNR} = \frac{\mu_{\text{signal}}}{\sigma_{\text{noise}}}, \qquad \mathrm{CNR} = \frac{\lvert\mu_A - \mu_B\rvert}{\sigma_{\text{noise}}}. \tag{89} \end{equation}\]

SNR measures detectability against background; CNR measures the distinguishability of two tissues, which is usually the clinically relevant quantity, since a lesion is found by its contrast with surrounding tissue rather than by its absolute intensity.

Each modality has a characteristic scaling:

• MRI: \(\mathrm{SNR}\propto (\text{voxel volume})\times\sqrt{N_{\mathrm{EX}}}\). Halving each linear voxel dimension reduces voxel volume eightfold, so 64 averages are needed to restore the original SNR.
• CT: \(\mathrm{SNR}\propto\sqrt{D}\) with \(D\) the dose. Halving the dose costs a factor \(\sqrt{2}\) in SNR, which can only be recovered by better detectors or better reconstruction, not by post-processing alone.
• Fluorescence microscopy: \(\mathrm{SNR}\propto\sqrt{n_{\text{photons}}}\), but longer exposure risks photobleaching and phototoxicity, so the trade-off is bounded on both sides.

lam <- 10^seq(1, 6, length.out = 400)
q1 <- ggplot(data.frame(lambda = lam, snr = sqrt(lam)), aes(lambda, snr)) +
  geom_line(colour = bpad_pal[1], linewidth = 0.95) +
  scale_x_log10() + scale_y_log10() +
  labs(title = "SNR = sqrt(lambda)", x = "expected counts (log)", y = "SNR (log)")

targ <- data.frame(SNR = c(10, 20, 40, 80))
targ$required_counts <- targ$SNR^2
q2 <- ggplot(targ, aes(factor(SNR), required_counts)) +
  geom_col(fill = bpad_pal[2], width = 0.6) +
  geom_text(aes(label = format(required_counts, big.mark = ",")),
            vjust = -0.4, size = 3.2) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.15))) +
  labs(title = "Counts required for a target SNR",
       x = "target SNR", y = "required counts")
bpad_grid(q1, q2, ncol = 2)

Figure 29: Figure 1.29: The square-root law. Left: SNR grows as the square root of the expected counts, so each doubling of SNR costs a fourfold increase in dose or time. Right: the same relation stated as the counts required to reach a target SNR.

From equation to dose. \(\mathrm{SNR}\propto\sqrt{N}\) is the single most consequential relation in low-dose imaging. It says that image quality cannot be bought cheaply: a 2-fold SNR gain costs a 4-fold dose increase. The only escapes are physical (more efficient detectors, better geometric efficiency) or algorithmic (priors and learned reconstruction that add information not present in the counts). Chapters 5, 6, and 8 are, in large part, an extended discussion of those two escapes.

Section 1.15 summary.

• Gaussian noise is additive and signal-independent; Poisson noise has \(\sigma = \sqrt{\lambda}\) and is signal-dependent; Rician noise arises from taking the magnitude of complex Gaussian data.
• Rician magnitudes are biased high at low SNR, with \(\mathbb{E}[M] = \sigma\sqrt{\pi/2}\) even when the true signal is zero.
• SNR measures detectability; CNR measures tissue distinguishability and is usually the clinically relevant quantity.
• \(\mathrm{SNR}\propto\sqrt{N}\): doubling SNR costs a fourfold increase in counts, dose, or time.

1.16 Statistical Estimation, Error Propagation, and Inference

Having a noise model is not enough; we must extract parameters from noisy data and state how confident we are. This section covers estimators and their properties, maximum likelihood, propagation of uncertainty, confidence intervals, and hypothesis testing.

1.16.1 Estimators and their properties

An estimator \(\hat{\theta}\) is a function of the data used to infer a parameter \(\theta\). Its quality is described by:

• Bias: \(\mathrm{Bias}(\hat\theta) = \mathbb{E}[\hat\theta] - \theta\); zero bias means correct on average.
• Variance: \(\mathrm{Var}(\hat\theta)\); lower variance means more reproducible.
• Mean squared error: \(\mathrm{MSE}(\hat\theta) = \mathrm{Bias}^2 + \mathrm{Var}(\hat\theta)\).
• Consistency: \(\hat\theta_n \to \theta\) in probability as \(n\to\infty\).
• Efficiency: among unbiased estimators, minimum variance, attaining the Cramér–Rao lower bound.

The MSE decomposition is the bias–variance trade-off, and it is worth internalizing here because Chapter 8 rests on it entirely: a slightly biased estimator with much lower variance can have smaller total error than the unbiased one. That is the entire justification for regularization, shrinkage, and smoothing.

1.16.2 Maximum likelihood estimation

Given data \(\{x_1,\ldots,x_n\}\) drawn from \(p(x\mid\theta)\), the likelihood is

\[\begin{equation} L(\theta) = \prod_{i=1}^{n}p(x_i\mid\theta), \tag{90} \end{equation}\]

and the maximum likelihood estimator maximizes it, or equivalently the log-likelihood,

\[\begin{equation} \hat{\theta}_{\mathrm{MLE}} = \arg\max_{\theta}\ \sum_{i=1}^{n}\log p(x_i\mid\theta). \tag{91} \end{equation}\]

Standard results follow directly:

• Gaussian data: \(\hat\mu_{\mathrm{MLE}} = \bar{x}\) and \(\hat\sigma^2_{\mathrm{MLE}} = \frac{1}{n}\sum(x_i-\bar{x})^2\) (note the \(1/n\): the MLE of the variance is biased, which is why \(1/(n-1)\) is used in Eq. 1.79).
• Poisson data: \(\hat\lambda_{\mathrm{MLE}} = \bar{x}\).
• Exponential decay with additive Gaussian error: no closed form; numerical optimization is required.

MLEs are asymptotically unbiased, consistent, efficient, and normally distributed:

\[\begin{equation} \hat{\theta}_{\mathrm{MLE}} \xrightarrow{\ d\ } \mathcal{N}\!\left(\theta,\ I(\theta)^{-1}\right), \qquad I(\theta) = -\mathbb{E}\!\left[\frac{\partial^{2}\log L}{\partial\theta^{2}}\right], \tag{92} \end{equation}\]

with \(I(\theta)\) the Fisher information. Equation (92) is what supplies standard errors for fitted parameters, and the Poisson likelihood is what makes the MLEM algorithm of Chapter 6 the statistically correct reconstruction method for emission tomography.

1.16.3 Error propagation

When a quantity \(f\) is computed from measured variables \(x_1,\ldots,x_n\) with uncertainties \(\sigma_i\), a first-order Taylor expansion (Section 1.2) gives

\[\begin{equation} \sigma_f^{2} = \sum_{i=1}^{n}\left(\frac{\partial f}{\partial x_i}\right)^{2}\sigma_i^{2} + 2\sum_{i<j}\frac{\partial f}{\partial x_i}\frac{\partial f}{\partial x_j}\,\mathrm{Cov}(x_i,x_j). \tag{93} \end{equation}\]

For uncorrelated variables the cross terms vanish. Frequently used special cases:

Function	Propagated uncertainty
\(f = x \pm y\)	\(\sigma_f^2 = \sigma_x^2 + \sigma_y^2\)
\(f = xy\) or \(x/y\)	\((\sigma_f/f)^2 \approx (\sigma_x/x)^2 + (\sigma_y/y)^2\)
\(f = x^{n}\)	\(\sigma_f/\lvert f\rvert = \lvert n\rvert\,\sigma_x/\lvert x\rvert\)
\(f = e^{ax}\)	\(\sigma_f = \lvert a\rvert f\,\sigma_x\)

Do not drop the covariance term. Fitted parameters are almost always correlated. In the exponential fit below, the amplitude and rate constant have correlation \(\approx 0.6\), and a Monte Carlo simulation that samples them independently gives an AUC uncertainty roughly 50% larger than the correct value. Ignoring the covariance is one of the most common errors in reported biomedical uncertainty.

Applications: the uncertainty of a PET standardized uptake value combines uncertainties in the region-of-interest activity, the injected dose calibration, and the patient weight; the uncertainty of a tumour volume from a segmented image combines boundary uncertainty with voxel-size calibration (Chapter 7); and the uncertainty of a radiomic feature under a one-voxel boundary perturbation is exactly Eq. (93) applied to the feature definition (Chapter 8).

1.16.4 Confidence intervals

A \((1-\alpha)\times 100\%\) confidence interval is a procedure that, on repeated application, contains the true parameter that fraction of the time. Standard forms:

\[\begin{equation} \bar{x} \pm z_{\alpha/2}\frac{\sigma}{\sqrt{n}} \quad (\sigma\text{ known}), \qquad \bar{x} \pm t_{\alpha/2,\,n-1}\frac{s}{\sqrt{n}} \quad (\sigma\text{ unknown}), \qquad \hat{p} \pm z_{\alpha/2}\sqrt{\frac{\hat p(1-\hat p)}{n}}. \tag{94} \end{equation}\]

The standard error \(\sigma/\sqrt{n}\) carries the practical message: halving the margin of error requires quadrupling the sample size. The same \(n^{-1/2}\) scaling that governs Poisson imaging noise governs clinical trial size, and for the same reason.

1.16.5 Hypothesis testing

The procedure is: state \(H_0\) (usually “no effect”) and \(H_1\); choose a significance level \(\alpha\); compute a test statistic; compute the \(p\)-value, the probability under \(H_0\) of a statistic at least as extreme as the one observed; and reject \(H_0\) if \(p<\alpha\).

Common tests in biomedical physics: the two-sample \(t\)-test (tumour versus normal SUV), the paired \(t\)-test (before versus after intervention in the same subjects), the chi-squared test (contingency tables of diagnostic accuracy), ANOVA (several groups), and the Kolmogorov–Smirnov test (comparing whole intensity distributions rather than just their means).

Two error types govern study design: Type I (false positive, rate \(\alpha\)) and Type II (false negative, rate \(\beta\)), with power \(= 1-\beta\). Trials are typically designed for power \(\ge 0.80\).

What a \(p\)-value is not. It is not the probability that \(H_0\) is true, and it is not a measure of effect size. A very small \(p\)-value from a large study can accompany a clinically negligible difference. Report effect sizes with confidence intervals; Chapter 8 develops this into the distinction between discrimination, calibration, and clinical utility.

set.seed(42)
true_A <- 100; true_k <- 0.15; sd_noise <- 6
t_obs  <- seq(0, 24, by = 2)
C_obs  <- true_A*exp(-true_k*t_obs) + rnorm(length(t_obs), 0, sd_noise)

## Nonlinear least squares = maximum likelihood under additive Gaussian noise
fit <- nls(C_obs ~ A*exp(-kk*t_obs), start = list(A = 90, kk = 0.1))
cf  <- coef(fit); V <- vcov(fit)
round(summary(fit)$coefficients, 5)

##     Estimate Std. Error  t value Pr(>|t|)
## A  101.93397    5.72783 17.79626        0
## kk   0.13622    0.01254 10.86272        0

cat(sprintf("correlation between A and k: %.3f\n", cov2cor(V)[1, 2]))

## correlation between A and k: 0.611

round(confint.default(fit), 4)

##      2.5 %   97.5 %
## A  90.7076 113.1603
## kk  0.1116   0.1608

## Delta method (Eq. 1.93) for AUC = A/k
grad     <- c(1/cf["kk"], -cf["A"]/cf["kk"]^2)
sd_delta <- sqrt(as.numeric(t(grad) %*% V %*% grad))

## Monte Carlo with the FULL covariance, and (wrongly) ignoring it
set.seed(1); nmc <- 20000
draws_corr <- sweep(matrix(rnorm(2*nmc), ncol = 2) %*% chol(V), 2, cf, "+")
auc_corr   <- draws_corr[, 1]/draws_corr[, 2]
auc_indep  <- rnorm(nmc, cf["A"], sqrt(V[1,1])) / rnorm(nmc, cf["kk"], sqrt(V[2,2]))

cat(sprintf("\nAUC point estimate      : %.1f\n", cf["A"]/cf["kk"]))

## 
## AUC point estimate      : 748.3

cat(sprintf("delta-method SD         : %.2f\n", sd_delta))

## delta-method SD         : 54.52

cat(sprintf("Monte Carlo SD (correct): %.2f\n", sd(auc_corr)))

## Monte Carlo SD (correct): 56.29

cat(sprintf("Monte Carlo SD (ignoring covariance, WRONG): %.2f\n", sd(auc_indep)))

## Monte Carlo SD (ignoring covariance, WRONG): 82.94

cat(sprintf("95%% CI for AUC          : [%.1f, %.1f]\n",
            quantile(auc_corr, 0.025), quantile(auc_corr, 0.975)))

## 95% CI for AUC          : [653.3, 874.0]

tf <- seq(0, 24, length.out = 300)
e1 <- ggplot() +
  geom_point(data = data.frame(t = t_obs, C = C_obs), aes(t, C),
             colour = bpad_pal[1], size = 2.2) +
  geom_line(data = data.frame(t = tf, C = cf["A"]*exp(-cf["kk"]*tf)),
            aes(t, C), colour = bpad_pal[2], linewidth = 0.9) +
  geom_line(data = data.frame(t = tf, C = true_A*exp(-true_k*tf)),
            aes(t, C), colour = "grey45", linetype = "dashed") +
  labs(title = sprintf("fit: k = %.4f (true %.2f)", cf["kk"], true_k),
       x = "time (h)", y = "concentration")

e2 <- ggplot(data.frame(A = draws_corr[, 1], k = draws_corr[, 2]), aes(A, k)) +
  geom_point(alpha = 0.06, colour = bpad_pal[1], size = 0.5) +
  geom_point(data = data.frame(A = cf["A"], k = cf["kk"]), aes(A, k),
             colour = bpad_pal[2], size = 2.6) +
  labs(title = "joint parameter uncertainty", x = "amplitude A", y = "rate k")

e3 <- ggplot(data.frame(auc = auc_corr), aes(auc)) +
  geom_histogram(bins = 60, fill = bpad_pal[3], colour = "white", linewidth = 0.1) +
  geom_vline(xintercept = quantile(auc_corr, c(0.025, 0.975)),
             linetype = "dotted", colour = bpad_pal[2], linewidth = 0.8) +
  labs(title = "propagated uncertainty in AUC = A/k", x = "AUC", y = "count")

bpad_grid(e1, e2, e3, ncol = 3)

Figure 30: Figure 1.30: Estimation and uncertainty for an exponential decay. Left: noisy observations with the fitted curve. Centre: the joint uncertainty of the two fitted parameters is an inclined ellipse, showing that amplitude and rate are correlated. Right: the Monte Carlo distribution of the derived quantity AUC = A/k, obtained by sampling the correlated parameters.

Section 1.16 summary.

• Estimators are judged by bias, variance, and their sum in the MSE: the bias–variance trade-off that justifies regularization.
• Maximum likelihood chooses the parameter that makes the observed data most probable; Fisher information supplies its standard error.
• Uncertainty propagates through a first-order Taylor expansion; the covariance term is essential and is routinely and wrongly omitted.
• Confidence intervals shrink as \(n^{-1/2}\): halving the margin of error quadruples the required sample size.
• A \(p\)-value is not an effect size and is not the probability that the null hypothesis is true.

Checkpoint 1.15. In Figure 1.30, the delta method and the correlated Monte Carlo agree, but the uncorrelated Monte Carlo does not. Explain why, and state what would have to be true of the experiment for the independent calculation to be correct.

1.17 Linear Algebra: Systems, Eigenvalues, and the SVD

Linear algebra is how large problems are made computable. Almost every natural process admits a first-order (linear) approximation, and once a problem is written as \(A\mathbf{x} = \mathbf{b}\) or \(A\mathbf{v} = \lambda\mathbf{v}\), an enormous body of theory and highly optimized software becomes available. Further worked examples are in the DSPA linear algebra chapter.

1.17.1 Systems of linear equations

The system

\[a + b + 2c = 6, \qquad 3a - 2b + c = 2, \qquad 2a + b - c = 3\]

is written compactly as \(A\mathbf{x} = \mathbf{b}\):

\[\begin{bmatrix} 1 & 1 & 2\\ 3 & -2 & 1\\ 2 & 1 & -1\end{bmatrix} \begin{bmatrix} a\\ b\\ c\end{bmatrix} = \begin{bmatrix} 6\\ 2\\ 3\end{bmatrix}.\]

Multiplying by the inverse gives \(\mathbf{x} = A^{-1}\mathbf{b}\), where \(A^{-1}A = I\) and \(I\) is the identity matrix, the multiplicative unit of matrix algebra.

A <- matrix(c(1, 1, 2,
              3, -2, 1,
              2, 1, -1), nrow = 3, byrow = TRUE)
b <- c(6, 2, 3)

solve(A)          %*% b        # explicit inverse

##      [,1]
## [1,] 1.35
## [2,] 1.75
## [3,] 1.45

solve(A, b)                    # numerically preferred: LU factorization

## [1] 1.35 1.75 1.45

all.equal(as.vector(solve(A) %*% b), as.vector(solve(A, b)))

## [1] TRUE

Never invert a matrix to solve a system. solve(A, b) factorizes \(A\) and back-substitutes; solve(A) %*% b forms the full inverse first. The two agree here, but the second is slower and numerically less stable, and for large or ill-conditioned systems (image reconstruction, regression with many correlated features) the difference becomes decisive. The same principle reappears in Chapter 8 as the reason regression software does not compute \((X^TX)^{-1}\) explicitly.

1.17.2 Special matrices

• Symmetric (\(A = A^{T}\)): real eigenvalues and orthogonal eigenvectors. Every covariance matrix and every diffusion tensor is symmetric, which is why their eigendecompositions are so well behaved.
• Positive definite (\(\mathbf{x}^{T}A\mathbf{x}>0\) for all \(\mathbf{x}\ne 0\)): all eigenvalues positive, positive determinant, and an efficient Cholesky factorization \(A = LL^{T}\) (used in Section 1.16 to draw correlated random samples).
• Orthogonal (\(Q^{T}Q = I\)): rotations and reflections; they preserve lengths and angles, \(\lVert Q\mathbf{x}\rVert = \lVert\mathbf{x}\rVert\). Rigid image registration is exactly a search over orthogonal matrices plus translations.
• Sparse: mostly zeros. CT system matrices and graph Laplacians are sparse, and exploiting that structure is what makes iterative reconstruction feasible.

1.17.3 Eigenvalues and eigenvectors

If \(A\mathbf{v} = \lambda\mathbf{v}\) for \(\mathbf{v}\ne\mathbf{0}\), then \(\mathbf{v}\) is an eigenvector and \(\lambda\) the corresponding eigenvalue. Equivalently \((\lambda I - A)\mathbf{v} = \mathbf{0}\), so the eigenvalues are the roots of the characteristic polynomial \(\det(\lambda I - A) = 0\).

Eigenvectors identify the directions along which a linear transformation acts by pure stretching, compression, or reflection; eigenvalues give the amount. Key identities:

\[\begin{equation} \mathrm{tr}(A) = \sum_i\lambda_i, \qquad \det(A) = \prod_i\lambda_i, \qquad \mathrm{rank}(A) = \#\{\lambda_i \ne 0\}. \tag{95} \end{equation}\]

For symmetric \(A\) the eigenvectors are mutually orthogonal, so \(A = Q\Lambda Q^{T}\) with \(Q\) orthogonal.

Applications throughout the book.

• Diffusion tensor imaging (Chapter 4). The diffusion tensor \(D\) is a symmetric \(3\times3\) matrix. Its eigenvalues \(\lambda_1\ge\lambda_2\ge\lambda_3\) describe the shape of the diffusion ellipsoid and its principal eigenvector \(\mathbf{v}_1\) gives the local fibre orientation used in tractography. The rotation-invariant summaries are

\[\begin{equation} \mathrm{MD} = \frac{\lambda_1+\lambda_2+\lambda_3}{3} = \frac{\mathrm{tr}(D)}{3}, \qquad \mathrm{FA} = \sqrt{\frac{3}{2}}\,\frac{\sqrt{\sum_i(\lambda_i - \mathrm{MD})^2}}{\sqrt{\sum_i\lambda_i^2}}, \tag{96} \end{equation}\]

with \(\mathrm{FA} = 0\) for isotropic diffusion and \(\mathrm{FA}\to 1\) for perfectly linear diffusion. * Compartmental kinetics (Section 1.13). The macro rate constants \(\alpha\) and \(\beta\) are eigenvalues of the system matrix. * Stability analysis. Eigenvalues of the Jacobian determine whether an equilibrium is stable (negative real parts) or unstable. * Vibrational modes (Chapter 2). Normal modes are eigenvectors of the dynamical matrix; eigenvalues give squared vibrational frequencies, the Raman-active transitions. * PCA (Section 1.18). Principal components are eigenvectors of the covariance matrix.

## Diffusion tensor: symmetric 3 x 3, units mm^2/s
Dt <- matrix(c(1.5e-3, 0.3e-3, 0.1e-3,
               0.3e-3, 1.2e-3, 0.2e-3,
               0.1e-3, 0.2e-3, 0.8e-3), nrow = 3, byrow = TRUE)
ed  <- eigen(Dt)
lam <- ed$values

MD <- mean(lam)
FA <- sqrt(1.5*sum((lam - MD)^2)/sum(lam^2))

cat("eigenvalues (mm^2/s):", signif(lam, 4), "\n")

## eigenvalues (mm^2/s): 0.001725 0.001058 0.0007163

cat("principal diffusion direction:", round(ed$vectors[, 1], 4), "\n")

## principal diffusion direction: 0.8142 0.5431 0.2054

cat(sprintf("mean diffusivity MD = %.4g mm^2/s;  fractional anisotropy FA = %.4f\n",
            MD, FA))

## mean diffusivity MD = 0.001167 mm^2/s;  fractional anisotropy FA = 0.4140

cat(sprintf("trace check: tr(D) = %.4g equals sum of eigenvalues = %.4g\n",
            sum(diag(Dt)), sum(lam)))

## trace check: tr(D) = 0.0035 equals sum of eigenvalues = 0.0035

## Two limiting cases for calibration
fa_of <- function(l) { m <- mean(l); sqrt(1.5*sum((l - m)^2)/sum(l^2)) }
cat(sprintf("isotropic tensor (1,1,1)e-3   -> FA = %.3f\n", fa_of(rep(1e-3, 3))))

## isotropic tensor (1,1,1)e-3   -> FA = 0.000

cat(sprintf("prolate tensor (1.7,0.2,0.2)e-3 -> FA = %.3f\n",
            fa_of(c(1.7e-3, 0.2e-3, 0.2e-3))))

## prolate tensor (1.7,0.2,0.2)e-3 -> FA = 0.870

1.17.4 The singular value decomposition

Theorem 1.8 (SVD). Every \(m\times n\) matrix \(A\), square or rectangular, singular or not, factorizes as

\[\begin{equation} A = U\Sigma V^{T}, \tag{97} \end{equation}\]

with \(U\) (\(m\times m\)) orthogonal whose columns are the left singular vectors (eigenvectors of \(AA^{T}\)), \(V\) (\(n\times n\)) orthogonal whose columns are the right singular vectors (eigenvectors of \(A^{T}A\)), and \(\Sigma\) diagonal with non-negative singular values \(\sigma_1\ge\sigma_2\ge\cdots\ge 0\) (square roots of the eigenvalues of \(A^{T}A\)).

Two properties make the SVD indispensable:

• \(\mathrm{rank}(A)\) equals the number of non-zero singular values, and the ratio \(\sigma_1/\sigma_{\min}\) is the condition number, which quantifies how much an inverse problem amplifies noise.
• Eckart–Young theorem. Truncating to the largest \(r\) singular values gives the best possible rank-\(r\) approximation of \(A\) in the least-squares sense. Since real images are approximately low-rank while noise is spread across all ranks, truncation denoises.

set.seed(11)
N <- 64; true_rank <- 5
Um <- matrix(rnorm(N*true_rank), N); Vm <- matrix(rnorm(N*true_rank), N)
X  <- Um %*% diag(seq(5, 1, length.out = true_rank)) %*% t(Vm)
X  <- (X - min(X))/(max(X) - min(X))          # a smooth, low-rank "image"
Xn <- X + matrix(rnorm(N*N, 0, 0.12), N)      # add noise

s     <- svd(Xn)
recon <- function(r) s$u[, 1:r, drop = FALSE] %*% diag(s$d[1:r], r, r) %*%
                     t(s$v[, 1:r, drop = FALSE])
rmse  <- sapply(1:25, function(r) sqrt(mean((recon(r) - X)^2)))
best  <- which.min(rmse)

cat(sprintf("RMSE of noisy image        : %.4f\n", sqrt(mean((Xn - X)^2))))

## RMSE of noisy image        : 0.1182

cat(sprintf("best truncation rank       : %d  (true rank %d)\n", best, true_rank))

## best truncation rank       : 4  (true rank 5)

cat(sprintf("RMSE after truncation      : %.4f  (%.0f%% reduction)\n",
            min(rmse), 100*(1 - min(rmse)/sqrt(mean((Xn - X)^2)))))

## RMSE after truncation      : 0.0492  (58% reduction)

cat("first eight singular values:", round(s$d[1:8], 3), "\n")

## first eight singular values: 32.929 3.66 2.849 2.311 2.016 1.763 1.739 1.625

v1 <- ggplot(data.frame(k = 1:30, sv = s$d[1:30]), aes(k, sv)) +
  geom_col(fill = bpad_pal[1], width = 0.7) +
  labs(title = "singular value spectrum", x = "index k", y = "sigma_k")
v2 <- ggplot(data.frame(r = 1:25, rmse = rmse), aes(r, rmse)) +
  geom_line(colour = bpad_pal[2], linewidth = 0.9) + geom_point(size = 1.3) +
  geom_vline(xintercept = best, linetype = "dashed", colour = "grey40") +
  labs(title = sprintf("reconstruction error (minimum at rank %d)", best),
       x = "retained rank r", y = "RMSE vs. clean image")
bpad_grid(v1, v2, ncol = 2)

Figure 31: Figure 1.31: Low-rank denoising by SVD. Left: singular values, showing a small number of large components (signal) above a broad noise floor. Right: reconstruction error against retained rank, with a clear minimum. Retaining too few ranks discards signal; retaining too many readmits noise.

op <- par(mfrow = c(1, 3), mar = c(0.6, 0.6, 2.2, 0.6))
show_img(X,           "clean (rank 5)")
show_img(Xn,          "noisy")
show_img(recon(best), sprintf("SVD truncated to rank %d", best))

Figure 32: Figure 1.31: Low-rank denoising by SVD. Left: singular values, showing a small number of large components (signal) above a broad noise floor. Right: reconstruction error against retained rank, with a clear minimum. Retaining too few ranks discards signal; retaining too many readmits noise.

par(op)

Note that the optimal truncation rank can be below the true rank when the weakest components sit under the noise floor. This is not a failure of the method: it is the bias–variance trade-off of Section 1.16 appearing in matrix form. Discarding a weak true component costs a little bias but saves more variance.

1.17.5 PCA as eigendecomposition

For a centred data matrix \(X\) (\(n\) observations in rows, \(p\) features in columns), the sample covariance is

\[\begin{equation} \Sigma = \frac{1}{n-1}X^{T}X. \tag{98} \end{equation}\]

The principal components are the eigenvectors of \(\Sigma\), equivalently the right singular vectors \(V\) of \(X\); the eigenvalue \(\lambda_k\) is the variance captured by component \(k\). Projecting onto the top \(k\) components, \(Z = XV_k\), minimizes reconstruction error among all \(k\)-dimensional linear projections, again by Eckart–Young.

set.seed(3)
Sig <- matrix(c(2, 1.5, 1.5, 2), 2)
Xd  <- matrix(rnorm(400), ncol = 2) %*% chol(Sig)
Xc  <- scale(Xd, center = TRUE, scale = FALSE)

S <- cov(Xc); e <- eigen(S)
cat("sample covariance:\n"); print(round(S, 3))

## sample covariance:

##       [,1]  [,2]
## [1,] 1.930 1.366
## [2,] 1.366 1.895

cat("eigenvalues:", round(e$values, 3), "\n")

## eigenvalues: 3.278 0.547

cat("eigenvectors (columns):\n"); print(round(e$vectors, 3))

## eigenvectors (columns):

##        [,1]   [,2]
## [1,] -0.712  0.703
## [2,] -0.703 -0.712

cat(sprintf("variance explained by PC1: %.1f%%\n",
            100*e$values[1]/sum(e$values)))

## variance explained by PC1: 85.7%

## Equivalence of PCA via eigen(cov(X)) and via svd(X)
sv <- svd(Xc)
cat("SVD check, d^2/(n-1):", round(sv$d^2/(nrow(Xc) - 1), 3), "\n")

## SVD check, d^2/(n-1): 3.278 0.547

ax <- data.frame(
  x0 = 0, y0 = 0,
  x1 = e$vectors[1, ]*sqrt(e$values)*2,
  y1 = e$vectors[2, ]*sqrt(e$values)*2,
  pc = c("PC1", "PC2"))
ggplot(data.frame(x = Xc[, 1], y = Xc[, 2]), aes(x, y)) +
  geom_point(alpha = 0.4, colour = bpad_pal[7], size = 1.4) +
  geom_segment(data = ax, aes(x = x0, y = y0, xend = x1, yend = y1, colour = pc),
               inherit.aes = FALSE, linewidth = 1.1,
               arrow = grid::arrow(length = grid::unit(0.2, "cm"))) +
  scale_colour_manual(values = bpad_pal[1:2]) +
  coord_fixed() +
  labs(title = "Principal axes of a correlated bivariate sample",
       x = "feature 1 (centred)", y = "feature 2 (centred)")

Figure 33: Figure 1.32: PCA on correlated bivariate data. The principal axes are the eigenvectors of the covariance matrix; their lengths are proportional to the square roots of the eigenvalues. The first component captures most of the variance, so the two-dimensional cloud is nearly one-dimensional.

1.17.6 Image registration as a linear transformation

Registration aligns a moving image \(M\) to a fixed image \(F\) by finding a spatial transformation \(T\) that minimizes a dissimilarity measure, usually with a regularization penalty:

\[\begin{equation} \hat{T} = \arg\min_{T}\ \mathcal{L}\big(F,\ M\circ T\big) + \lambda\,\mathcal{R}(T). \tag{99} \end{equation}\]

The transformation families form a hierarchy of increasing flexibility:

1. Rigid (Euclidean): \(T(\mathbf{x}) = R\mathbf{x} + \mathbf{t}\) with \(R^{T}R = I\). Preserves all distances; 6 parameters in 3-D. Used for intra-subject alignment of the same anatomy.
2. Affine: \(T(\mathbf{x}) = A\mathbf{x}+\mathbf{t}\) with \(A\) invertible; adds scaling and shear, 12 parameters in 3-D. Used for atlas and cross-subject alignment.
3. Projective: \(T(\mathbf{x}) = (A\mathbf{x}+\mathbf{t})/(\mathbf{v}^{T}\mathbf{x}+1)\); accounts for perspective, relevant to endoscopic and microscope mosaicking.
4. Deformable: a dense non-parametric displacement field, developed in Chapter 7.

In homogeneous coordinates a 2-D affine transformation is a single \(3\times3\) matrix

\[\begin{equation} T = \begin{bmatrix} a_{11} & a_{12} & t_x\\ a_{21} & a_{22} & t_y\\ 0 & 0 & 1\end{bmatrix}, \tag{100} \end{equation}\]

so that composing transformations is matrix multiplication and inverting one is matrix inversion. That is the entire practical reason homogeneous coordinates are used.

Clinically, registration underpins multimodal fusion (PET with CT, MRI with ultrasound), image-guided intervention, longitudinal tumour tracking, and motion correction in dynamic imaging.

bilinear <- function(img, xs, ys) {
  n <- nrow(img); m <- ncol(img)
  x0 <- floor(xs); y0 <- floor(ys); dx <- xs - x0; dy <- ys - y0
  gv <- function(a, b) { a <- pmin(pmax(a, 1), n); b <- pmin(pmax(b, 1), m)
                         img[cbind(a, b)] }
  ok  <- xs >= 1 & xs <= n & ys >= 1 & ys <= m
  out <- gv(x0, y0)*(1-dx)*(1-dy) + gv(x0+1, y0)*dx*(1-dy) +
         gv(x0, y0+1)*(1-dx)*dy   + gv(x0+1, y0+1)*dx*dy
  out[!ok] <- 0
  out
}

## Resample "img" onto the fixed grid under p_moving = c + t + s R(theta) (p_fixed - c)
warp <- function(img, theta, s, tx, ty) {
  n <- nrow(img); m <- ncol(img); cx <- (n+1)/2; cy <- (m+1)/2
  g  <- expand.grid(i = 1:n, j = 1:m)
  dx <- g$i - cx; dy <- g$j - cy
  xs <- s*( cos(theta)*dx - sin(theta)*dy) + cx + tx
  ys <- s*( sin(theta)*dx + cos(theta)*dy) + cy + ty
  matrix(bilinear(img, xs, ys), n, m)
}

N  <- 64
gg <- expand.grid(i = 1:N, j = 1:N)
xr <- (gg$i - N/2)/N*2; yr <- (gg$j - N/2)/N*2
fixed <- matrix(as.numeric((xr/0.6)^2 + (yr/0.3)^2 <= 1), N, N)

theta0 <- pi/8; s0 <- 1.1; t0 <- c(5, -3)          # the transformation to recover
moving <- warp(warp(fixed, -theta0, 1/s0, 0, 0), 0, 1, -t0[1], -t0[2])

obj <- function(p) mean((fixed - warp(moving, p[1], p[2], p[3], p[4]))^2)
opt <- optim(c(0, 1, 0, 0), obj, method = "Nelder-Mead",
             control = list(maxit = 800, reltol = 1e-10))
registered <- warp(moving, opt$par[1], opt$par[2], opt$par[3], opt$par[4])

knitr::kable(data.frame(
  parameter = c("rotation (rad)", "scale", "translation x", "translation y"),
  true      = round(c(theta0, s0, t0), 4),
  recovered = round(opt$par, 4)),
  caption = "Registration parameters recovered by minimizing mean squared difference.")

Table 8: Registration parameters recovered by minimizing mean squared difference.

parameter	true	recovered
rotation (rad)	0.3927	0.3937
scale	1.1000	1.1006
translation x	5.0000	5.0109
translation y	-3.0000	-3.0017

cat(sprintf("MSE before registration: %.5f;  after: %.5f\n",
            mean((fixed - moving)^2), opt$value))

## MSE before registration: 0.10436;  after: 0.00160

op <- par(mfrow = c(1, 4), mar = c(0.6, 0.6, 2.2, 0.6))
show_img(fixed,                    "fixed F")
show_img(moving,                   "moving M")
show_img(abs(fixed - moving),      "residual before")
show_img(abs(fixed - registered),  "residual after")

Figure 34: Figure 1.33: Affine registration. A known rotation, scaling, and translation are applied to a synthetic object, then recovered by numerically minimizing the mean squared difference. The residual image after registration is nearly empty apart from interpolation error at the boundary.

par(op)

Section 1.17 summary.

• Write the problem as \(A\mathbf{x}=\mathbf{b}\) and solve by factorization, never by forming \(A^{-1}\).
• Eigenvectors are the directions a linear map merely stretches; eigenvalues give the amount, and they sum to the trace and multiply to the determinant.
• Symmetric matrices (covariances, diffusion tensors) have real eigenvalues and orthogonal eigenvectors; FA and MD are rotation-invariant functions of the eigenvalues.
• The SVD exists for every matrix and gives the optimal low-rank approximation, which is why truncation denoises and compresses.
• PCA is the eigendecomposition of the covariance matrix, equivalently the SVD of the centred data.
• Registration is optimization over a family of linear transformations, composed and inverted by matrix multiplication in homogeneous coordinates.

Checkpoint 1.16. In Figure 1.31 the optimal truncation rank was below the true rank of the underlying image. Explain this in terms of bias and variance, and predict what would happen to the optimal rank if the noise level were halved.

1.18 Dimensionality Reduction

Biomedical datasets are wide: a radiomic feature set has hundreds of columns, a gene-expression array tens of thousands, an fMRI volume \(10^5\) voxels. Dimensionality reduction constructs a small number of new variables that retain the structure that matters, making visualization possible, reducing overfitting, and stabilizing downstream models. Extended treatments appear in the DSPA dimensionality reduction chapter, and interactive demonstrations are available through the SOCR Statistical Data Analyzer.

1.18.1 Principal component analysis

PCA finds an orthogonal basis in which the data have maximum variance along successive, mutually uncorrelated directions. Writing the \(k\)-th component as

\[\begin{equation} \mathrm{pc}_k = \mathbf{a}_k^{T}X = \sum_{i=1}^{N}a_{i,k}X_i, \tag{101} \end{equation}\]

the weights \(\mathbf{a}_k\) are exactly the eigenvectors of the covariance matrix \(\Sigma\) (Section 1.17.5), and the variance along component \(k\) is the eigenvalue \(\lambda_k\).

PCA is linear, unsupervised, and variance-driven. Those three words also define its limitations: it cannot represent curved structure, it does not know what you are trying to predict, and a high-variance direction is not necessarily an informative one (a scanner-intensity drift can dominate PC1 while carrying no biological signal at all, a failure mode taken up as batch effect and harmonization in Chapter 8).

Biomedical uses: reducing voxel time series in fMRI; finding dominant patterns in gene expression; identifying modes of anatomical shape variation in morphometrics; compressing correlated radiomic features; and quality control, where PCA of calibration data exposes scanner drift.

1.18.2 Independent component analysis

ICA seeks components that are statistically independent rather than merely uncorrelated, which requires attending to higher-order moments (skewness, kurtosis) instead of variance alone. The generative model is

\[\begin{equation} X_i = a_{i,1}s_1 + \cdots + a_{i,n}s_n, \tag{102} \end{equation}\]

with the sources \(s_j\) recovered by maximizing their mutual independence (equivalently, minimizing mutual information).

The distinction from PCA is best seen in its canonical applications. In fMRI, ICA separates spatially independent functional networks without prior knowledge of their extent; in EEG and MEG, it separates artifact sources (eye blinks, cardiac interference, line noise) from neural activity. Neither task is a variance-maximization problem, and PCA does not solve either well.

1.18.3 t-SNE and UMAP

t-distributed stochastic neighbour embedding (t-SNE) is a nonlinear method for embedding high-dimensional data in two or three dimensions for visualization. It proceeds in two steps: construct a probability distribution over pairs of high-dimensional points such that similar points have high pairing probability; then find a low-dimensional configuration whose corresponding distribution minimizes the Kullback–Leibler divergence from the first.

Uniform manifold approximation and projection (UMAP) (McInnes and Healy, 2018) builds a weighted graph of high-dimensional neighbourhoods and optimizes a low-dimensional layout that preserves its topology. It tends to retain more global structure than t-SNE and is substantially faster, which has made it the default for single-cell RNA sequencing and large imaging cohorts.

How to read, and not read, a t-SNE or UMAP plot. These are visualization tools, not measurement tools. Distances between clusters are not meaningful, cluster sizes are not meaningful, and the layout depends on hyperparameters (perplexity, number of neighbours) and on the random seed. Apparent separation can be an artifact of the embedding; genuine separation must be confirmed with a quantitative test in the original feature space. Reporting a t-SNE plot as evidence of group difference is a recurring error in the imaging-AI literature that Chapter 8 addresses directly.

Comparing the four methods.

Method	Linear?	Objective	Preserves	Best used for
PCA	yes	maximize variance	global variance, distances	compression, denoising, feature reduction before modelling
ICA	yes (mixing)	maximize independence	statistically independent sources	source separation, artifact removal
t-SNE	no	minimize KL divergence	local neighbourhoods	visualizing cluster structure
UMAP	no	preserve fuzzy topology	local and some global structure	visualizing large cohorts, faster than t-SNE

Section 1.18 summary.

• PCA finds orthogonal directions of maximum variance and is the eigendecomposition of the covariance matrix.
• ICA finds statistically independent sources using higher-order moments, which is what makes it the right tool for artifact and network separation.
• t-SNE and UMAP are nonlinear embeddings for visualization: local structure is meaningful, global distances and cluster sizes are not.
• Dimensionality reduction is the first step of nearly every biomedical machine-learning pipeline, and its failure modes propagate into everything downstream.

1.19 Choosing the Right Tool

The most common failure in applying mathematics to a biomedical problem is not an algebraic error; it is selecting a method whose assumptions the data do not satisfy. The guide below matches problem structure to mathematical structure.

Decision guide.

• Continuous change in one variable \(\rightarrow\) differential or integral calculus ((1)). Several variables varying together \(\rightarrow\) partial derivatives and PDEs ((12)).
• Time evolution of interacting quantities \(\rightarrow\) ODE systems ((13)). Linear systems reduce to an eigenvalue problem ((17)); nonlinear systems require phase-plane or numerical analysis.
• Geometric or vector-valued data \(\rightarrow\) vector and matrix operations ((3)). Rotation-invariant structure \(\rightarrow\) eigen-spectrum or SVD ((17)).
• Periodic or bandlimited signals \(\rightarrow\) Fourier methods ((8), (9)). Repeated sampling \(\rightarrow\) Nyquist analysis ((10)).
• A blurring instrument \(\rightarrow\) convolution and the PSF ((7)). Recovering the object \(\rightarrow\) regularized deconvolution.
• Oscillations, impedance, wave interference \(\rightarrow\) complex numbers and phasors ((6)).
• Counting discrete random events \(\rightarrow\) Poisson ((15)). Many small additive errors \(\rightarrow\) Gaussian. Magnitude of complex data \(\rightarrow\) Rician.
• Combining prior information with new evidence \(\rightarrow\) Bayes’ theorem ((14)).
• Estimating a parameter from noisy data \(\rightarrow\) maximum likelihood ((16)), with error propagation and a confidence interval always reported alongside.
• Many variables, few underlying factors \(\rightarrow\) PCA for variance, ICA for independent sources, t-SNE or UMAP for visualization only ((18)).

The table below summarizes the operating regime and the failure mode of each major technique. The failure column is the more useful one: knowing when a method breaks is what distinguishes competent from mechanical application.

Technique	Natural setting	Fails when
Differentiation	smooth analytic signals	data are noisy (amplifies noise; smooth first)
Numerical integration	any sampled curve	the sampling is too coarse or the integrand is singular
Taylor expansion	local behaviour near a point	far from the expansion point, or the function is not smooth
Fourier transform	stationary, periodic, or bandlimited signals	the spectrum changes over the record (use time–frequency methods)
Nyquist criterion	bandlimited signals	the signal is not bandlimited (aliasing is then unavoidable)
Convolution model	linear, shift-invariant systems	the PSF varies across the field of view, or the system saturates
Exponential ODE model	first-order rate processes	multiple compartments or saturable kinetics are present
Gaussian statistics	many additive small errors	counts are low (use Poisson) or data are magnitudes (use Rician)
Least squares	Gaussian, homoscedastic errors	errors are heteroscedastic or heavy-tailed
Linearization by logarithm	multiplicative noise	noise is additive (it then distorts the error structure)
PCA	linear, variance-dominated structure	signal is low-variance, or structure is nonlinear
t-SNE / UMAP	visualizing local neighbourhood structure	inter-cluster distances or cluster sizes are interpreted
SVD truncation	approximately low-rank data	signal and noise share a similar singular-value scale

1.20 Conclusions and Discussion

This chapter has assembled the mathematical and statistical apparatus of biomedical physics: single-variable calculus, series approximation, tensors, complex analysis, linear systems, Fourier methods, sampling theory, vector calculus and PDEs, ordinary differential equations, probability, noise modelling, statistical estimation, linear algebra, and dimensionality reduction. The physics, instrumentation, and data analysis of every subsequent chapter are written in this language.

Three themes deserve restatement, because they recur in every chapter that follows.

One equation, many faces. Exponential decay describes radioactive decay (Chapter 6), drug clearance (Section 1.13), fluorescence lifetime (Chapter 2), Beer–Lambert attenuation (Chapters 2 and 5), and \(T_1\)/\(T_2\) relaxation (Chapter 4). The Fourier transform links an ECG to its heart-rate-variability spectrum, MRI \(k\)-space to image space, an FTIR interferogram to a molecular spectrum, and an OCT spectral interferogram to a depth profile. The diffusion equation governs Brownian motion, drug distribution, thermal transport, and the apparent diffusion coefficient that detects stroke within minutes. Recognizing these shared roots is what allows a student to move between modalities instead of memorizing each in isolation.

Every method has a regime. Differentiating a smooth function is exact; differentiating a noisy signal amplifies noise. A wider Fourier analysis window improves frequency resolution and degrades time resolution. PCA preserves variance but not local structure; t-SNE preserves local structure but distorts global distance. Truncating an SVD removes noise and, past the optimum, removes signal. Gaussian statistics are correct until the counts get low, and then Poisson statistics are. The scientific and clinical skill is knowing which regime you are in, which is why Section 1.19 tabulates failure modes alongside natural settings.

Uncertainty is part of the answer. A parameter without an uncertainty is not a measurement, it is a number. The chain from noise model (Section 1.15) through maximum likelihood (Section 1.16) to error propagation and confidence intervals is what converts a fitted curve into a quantitative result, and Chapter 8 extends the same chain from single parameters to whole predictive models.

Key points

• The derivative is an instantaneous rate, the integral an accumulated total; the fundamental theorem makes them inverses. Critical points locate extrema, inflection points locate extrema of the rate.
• Taylor expansion converts smooth functions into polynomials with bounded error, generating finite-difference schemes whose accuracy is limited by round-off as well as truncation.
• Scalars, vectors, matrices, and higher-order arrays are the data structures of biomedical imaging; a dynamic volume is a fourth-order array whose slices are the images and time courses actually analyzed.
• Complex numbers carry amplitude and phase; phasors turn the calculus of oscillations into algebra, and underpin bioimpedance, interferometry, and the complex MR signal.
• An LTI system is characterized by its impulse response; the imaging equation \(g = f*h+n\) states that a measurement is the object convolved with the instrument’s PSF plus noise.
• Fourier series decompose periodic signals; the Fourier transform generalizes to non-periodic ones; the 2-D transform reconstructs MRI and CT images. Truncation produces Gibbs ringing.
• The Nyquist criterion \(f_s > 2f_{\max}\) applies identically in time and space, and aliasing, once introduced, is irreversible.
• First-order linear ODEs have exponential solutions; logistic and Gompertz models add saturation; compartmental models turn dosing into a linear-algebra problem in which macro rate constants are eigenvalues and \(\mathrm{AUC} = \mathrm{Dose}/CL\).
• Probability follows from three axioms; Bayes’ theorem converts test characteristics into patient-specific probabilities, and predictive value depends on prevalence.
• Gaussian noise is signal-independent, Poisson noise satisfies \(\sigma = \sqrt{\lambda}\), and Rician noise biases magnitude images upward at low SNR. \(\mathrm{SNR}\propto\sqrt{N}\) sets the dose–quality trade-off.
• Maximum likelihood, error propagation (including the covariance term), and confidence intervals are the working tools of quantitative biomedical research.
• Eigendecomposition and the SVD underlie diffusion tensor metrics, low-rank denoising, PCA, and registration; the SVD gives the optimal low-rank approximation.
• PCA, ICA, t-SNE, and UMAP reduce dimensionality for compression, source separation, and visualization respectively, and each has a characteristic failure mode.

Connections to other BPAD1 chapters

• Chapter 2 (Optical Methods in Medical Diagnosis). This chapter supplies the exponential-decay solution behind the Beer–Lambert law and fluorescence lifetime ((13)), the Fourier transform used in FTIR and Fourier-domain OCT ((9)), the definite integral that produces the Stefan–Boltzmann law from the Planck function ((1)), the critical-point analysis that yields Wien’s displacement law (1.1.7), the phasor description of low-coherence interferometry ((6)), and the reciprocal width relation that ties OCT axial resolution to optical bandwidth (1.9.4).
• Chapter 3 (Photoacoustic Methods). The wave and diffusion equations ((12)) underlie acoustic propagation and photoacoustic generation; the convolution/PSF model ((7)) sets axial and lateral resolution and motivates deconvolution; the Nyquist criterion ((10)) governs Doppler velocity aliasing; exponential attenuation ((13)) describes depth-dependent signal loss; and least-squares linear algebra ((17)) performs multispectral spectral unmixing.
• Chapter 4 (MRI). The Bloch equations are ODEs ((13)) with exponential \(T_1\) and \(T_2\) solutions; the received signal is a complex phasor ((6)); the image is the 2-D inverse Fourier transform of sampled \(k\)-space ((9)); field of view, wrap-around, and truncation ringing follow from sampling theory and Gibbs ((8), (10)); magnitude images carry Rician noise ((15)); and diffusion tensor metrics are eigenvalue computations ((17)).
• Chapter 5 (X-ray Imaging and CT). Exponential attenuation and half-value layer are first-order ODE solutions ((13)); Poisson statistics and the \(\sqrt{N}\) law set the dose–noise trade-off ((15)); the central slice theorem and the ramp filter are Fourier statements ((9)); the Radon transform is a family of line integrals ((1)); and iterative reconstruction is large-scale linear algebra with regularization ((17)).
• Chapter 6 (PET/SPECT Imaging). The decay law and effective half-life are exponential ODEs ((13)); count statistics, noise-equivalent count rate, and detectability are Poisson ((15)); MLEM reconstruction is maximum likelihood under a Poisson model ((16)); tracer kinetic and Patlak models are compartmental ODE systems ((13)); attenuation correction applies the exponential attenuation law; and SUV uncertainty is error propagation ((16)).
• Chapter 7 (General Medical Image Processing). Images are tensors ((3)); filtering is convolution ((7)); frequency-domain processing is the 2-D Fourier transform ((9)); resampling and interpolation are governed by sampling theory ((10)); registration is optimization over linear transformations ((17)); segmentation uses Gaussian mixtures ((14)); and quantification requires uncertainty propagation ((16)).
• Chapter 8 (Data Modeling, AI, and Machine Learning). Design matrices and network activations are tensors ((3)); training is gradient-based optimization ((1), (12)); losses are negative log-likelihoods chosen to match the noise model ((15), (16)); regularization is the bias–variance trade-off ((16)); convolutional layers are learned filters ((7)); PCA and SVD compress radiomic features ((17), (18)); Bayes’ theorem converts model output into post-test probability ((14)); sampling and resolution determine which features are even learnable ((10)); and longitudinal biomarkers are fitted with the growth ODEs of Section 1.13.

Glossary

Term	Meaning
Aliasing	Irreversible folding of frequencies above \(f_s/2\) into the observable band
Antiderivative	A function whose derivative is the integrand; unique up to a constant
Bias–variance trade-off	Decomposition \(\mathrm{MSE} = \mathrm{Bias}^2 + \mathrm{Var}\); the justification for regularization
Central limit theorem	Sums of many independent finite-variance terms tend to a normal distribution
Condition number	\(\sigma_1/\sigma_{\min}\); how strongly an inverse problem amplifies noise
Convolution	Weighted superposition \((x*h)(t)\); the output of an LTI system
Convolution theorem	Convolution in one domain is multiplication in the other
Critical point	Where \(f' = 0\); candidate location of an extremum
Curl	\(\nabla\times\mathbf{F}\); local rotation of a vector field
Divergence	\(\nabla\cdot\mathbf{F}\); local expansion of a vector field
Divergence theorem	Interior accumulation of divergence equals boundary flux
Eigenvalue / eigenvector	\(A\mathbf{v} = \lambda\mathbf{v}\); direction merely stretched, and the stretch factor
Euler’s formula	\(e^{i\theta} = \cos\theta + i\sin\theta\)
Fisher information	Curvature of the log-likelihood; inverse of the asymptotic estimator variance
Fourier series	Expansion of a periodic function in harmonics of its fundamental
Fourier transform	Decomposition of a signal into a continuous frequency spectrum
Fractional anisotropy (FA)	Rotation-invariant measure of diffusion-tensor eigenvalue disparity
Gibbs phenomenon	About 9% overshoot of a truncated Fourier series near a discontinuity
Gradient	\(\nabla f\); vector of partial derivatives, pointing along steepest ascent
Hessian	Matrix of second partial derivatives; classifies multivariate critical points
Impulse response	System output for a delta-function input; determines all other outputs
Inflection point	Where \(f'' = 0\) with a sign change; the rate attains an extremum
\(k\)-space	The spatial-frequency (2-D Fourier) domain, sampled directly by MRI
Laplacian	\(\nabla^2 f\); departure of a value from its local neighbourhood average
Likelihood	Probability of the observed data viewed as a function of the parameter
L’Hôpital’s rule	Resolves \(0/0\) and \(\infty/\infty\) limits via the ratio of derivatives
Maximum likelihood estimate	Parameter value maximizing the likelihood
Mean value theorem	An average rate is attained as an instantaneous rate somewhere in the interval
Nyquist frequency	\(f_s/2\); the highest unambiguously representable frequency
Orthogonality	Zero inner product; what makes Fourier coefficients independent projections
Parseval’s theorem	Total power is equal in the time and frequency domains
Phasor	Complex number encoding the amplitude and phase of a sinusoid
Point spread function (PSF)	Two-dimensional impulse response of an imaging system
Poisson noise	Counting noise with \(\sigma = \sqrt{\lambda}\), hence \(\mathrm{SNR} = \sqrt{\lambda}\)
Positive predictive value	\(P(\text{disease}\mid\text{positive test})\); depends on prevalence
Rician noise	Distribution of the magnitude of complex Gaussian data; biased high at low SNR
Riemann sum	Finite approximation whose limit defines the definite integral
Sifting property	\(\int f(t)\delta(t-t_0)dt = f(t_0)\)
Singular value decomposition	\(A = U\Sigma V^{T}\); exists for every matrix; optimal low-rank truncation
SNR / CNR	Signal-to-noise and contrast-to-noise ratios
Taylor series	Polynomial expansion of a smooth function with an explicit remainder
Tensor	Multi-dimensional array; in physics, one that transforms covariantly
Time constant	\(\tau = 1/k\); time for exponential decay to reach \(1/e\) of its initial value
Wiener–Khinchin theorem	The power spectrum is the Fourier transform of the autocovariance

Review Questions

1. State the \(\epsilon\)–\(\delta\) definition of a limit and explain why \(f(a)\) may be undefined even when \(\lim_{x\to a}f(x)\) exists.
2. Use L’Hôpital’s rule to evaluate \(\lim_{x\to 0}\sin(x)/x\), and identify where this limit appears in sampling theory.
3. State the fundamental theorem of calculus. Why does the definite integral of a derivative equal the difference of boundary values?
4. For \(f(x) = x^3 - 5x\), find all critical points, classify each, and locate the inflection point. Why is the inflection point not a critical point?
5. Why is \(\int_0^{\infty}e^{-x}dx\) an improper integral, and what is its value?
6. Explain the difference between a scalar, a vector, a matrix, and a higher-order array, giving one biomedical example of each.
7. Why is the central finite difference more accurate than the forward difference, and why does making \(h\) arbitrarily small nevertheless degrade a numerical derivative?
8. Write a system of three linear equations as \(A\mathbf{x} = \mathbf{b}\) and explain why solve(A, b) is preferable to solve(A) %*% b.
9. What is the relationship among the trace, the determinant, and the eigenvalues of a matrix?
10. Why is a free-falling projectile’s velocity linear in time while its displacement is quadratic?
11. Compare the growth of \(\log x\), \(x\), \(x^2\), and \(e^x\) at large \(x\), and give a biological process for each.
12. State the Maclaurin series of \(\sin x\) and \(e^x\) and explain their role in numerical methods.
13. State Euler’s formula and explain its physical meaning for an oscillating signal.
14. What is a phasor, and why does it simplify the analysis of alternating-current and wave problems?
15. Define the impulse response of an LTI system, state the convolution theorem, and explain why the theorem is computationally valuable.
16. State Parseval’s theorem and explain what it licenses in physiological band-power analysis.
17. Why does the Gibbs overshoot fail to diminish as more harmonics are added, and what artifact does this produce in MRI?
18. State the Nyquist criterion. Why must an anti-aliasing filter be applied before sampling rather than after?
19. Why are spherical coordinates natural for the eye but cylindrical coordinates natural for a radiotherapy beam?
20. State the diffusion equation and interpret each term. Why does \(\langle x^2\rangle = 2Dt\) imply that organisms larger than a few hundred micrometres require circulation?
21. State the divergence theorem and explain in what sense it generalizes the fundamental theorem of calculus.
22. Show that \(dN/dt = kN\) has solution \(N(t) = N_0e^{kt}\) and derive the doubling time.
23. Contrast the exponential, logistic, and Gompertz growth models by their inflection points, and explain why they are hard to distinguish from early data.
24. In a two-compartment pharmacokinetic model, what are \(\alpha\) and \(\beta\) mathematically, and what does \(\mathrm{AUC} = \mathrm{Dose}/CL\) express?
25. State Bayes’ theorem and explain why the positive predictive value of a fixed test changes with the population screened.
26. Define the four statistical moments and state what each reveals about an emission-imaging intensity histogram.
27. What is the relationship between the mean and the variance of a Poisson distribution, and what practical law follows from it?
28. Distinguish Gaussian, Poisson, and Rician noise, and name the modality in which each dominates.
29. Why are magnitude MR images biased at low SNR, and what is the bias when the true signal is zero?
30. State the maximum-likelihood estimator for the mean of a Gaussian. Is the corresponding MLE of the variance unbiased?
31. Why must the covariance term be retained when propagating uncertainty through a function of fitted parameters?
32. Define the SVD and explain why it is called the most general matrix factorization.
33. Why does truncating an SVD denoise an image, and why can the optimal rank be lower than the true rank?
34. What is the difference between PCA and ICA, and which is appropriate for removing eye-blink artifact from EEG?
35. Compare t-SNE and UMAP, and state two features of such a plot that must not be interpreted.
36. Define the Hessian and explain how it generalizes the second-derivative test.
37. Why is the imaging equation \(g = f*h+n\) a useful organizing statement across four different modalities?
38. Explain why the SNR of a photon-counting measurement scales as the square root of the counts, and what that implies for dose reduction.

Problems

Problem 1. Critical points of a cubic

Find all critical points of \(f(x) = x^3 - 6x^2 + 9x + 1\), classify each, and locate any inflection point. Sketch the graph and mark all three points.

Problem 2. Definite integral by parts

Evaluate \(\int_0^1 x e^{-x}dx\) two ways: (a) by integration by parts, and (b) with integrate() in R. Confirm agreement to three decimal places.

Problem 3. Riemann sums and convergence order

For \(f(x) = e^{-x^2}\) on \([0,2]\), compute the left, midpoint, and trapezoid approximations for \(n = 10, 20, 40, 80\). Tabulate the absolute errors against integrate() and estimate the convergence order of each rule from the ratio of successive errors.

Problem 4. Matrix algebra in R

For \(A = \begin{bmatrix} 1&2&3\\ 4&5&6\\ 7&8&10\end{bmatrix}\) and \(\mathbf{b} = (1,2,3)^{T}\): (a) compute \(A^{-1}\), \(A^{T}\), \(\det(A)\), and \(\mathrm{tr}(A)\); (b) solve \(A\mathbf{x} = \mathbf{b}\); (c) verify by substitution; (d) confirm that \(\mathrm{tr}(A)\) equals the sum of the eigenvalues and \(\det(A)\) their product.

Problem 5. Euler’s identity and De Moivre

Compute \((1+i)^4\) (a) by repeated multiplication and (b) by writing \(1+i = \sqrt{2}e^{i\pi/4}\) and applying De Moivre’s theorem.

Problem 6. Phasor addition

Two current phasors at the same frequency are \(\tilde{I}_1 = 5e^{i\pi/6}\) A and \(\tilde{I}_2 = 3e^{-i\pi/4}\) A. (a) Find the resultant phasor. (b) Give its amplitude and phase. (c) Write the time-domain expression. (d) Explain why the amplitude is less than \(5+3\).

Problem 7. Convolution of two boxcars

Compute \(y = x * h\) where \(x(t) = h(t) = 1\) for \(0\le t\le 1\) and zero otherwise. Sketch \(y(t)\), state its support, and explain how the support followed from the convolution property alone.

Problem 8. Fourier transform of a Gaussian

Using \(\hat{f}(\omega) = \int_{-\infty}^{\infty}f(t)e^{-i2\pi\omega t}dt\), show that the Fourier transform of \(f(t) = e^{-\pi t^2}\) is \(\hat{f}(\omega) = e^{-\pi\omega^2}\). (Hint: complete the square.) State the consequence for the relation between pulse duration and bandwidth.

Problem 9. Nyquist sampling of a pediatric ECG

A pediatric ECG contains QRS energy up to 250 Hz and is digitized at 500 Hz. (a) What is the Nyquist frequency? (b) Is the QRS faithfully reconstructed? (c) At what frequency would a 300 Hz interference component appear? (d) Propose a sampling rate and pre-sampling filter with an adequate margin.

Problem 10. Diffusion in one dimension

An initial Gaussian profile \(C(x,0) = A\exp\!\big(-x^2/(4\sigma_0^2)\big)\) evolves under the diffusion equation into \(C(x,t) = \frac{A\sigma_0}{\sqrt{\sigma_0^2+Dt}}\exp\!\left(-\frac{x^2}{4(\sigma_0^2+Dt)}\right)\). (a) Verify that the total mass \(\int C\,dx\) is conserved. (b) Compute \(\langle x^2\rangle\) as a function of time. (c) With \(\sigma_0 = 1\) µm, the width doubles in 7.1 ms; estimate \(D\) and compare it with the literature value for fluorescein in water.

Problem 11. Tumour growth and doubling time

A tumour of \(N_0 = 10^6\) cells grows exponentially with \(k = 0.08\) day\(^{-1}\). (a) Write the growth equation. (b) How many cells after 30 days? (c) What is the doubling time? (d) With carrying capacity \(K = 10^{10}\), what does the logistic model predict at \(t = 30\) days, and why is it so close to the exponential prediction?

Problem 12. Two-compartment pharmacokinetics

A drug follows the two-compartment model with \(k_{10} = 0.2\), \(k_{12} = 0.4\), \(k_{21} = 0.3\) hr\(^{-1}\), after a 100 mg IV bolus into a central compartment of volume \(V_1 = 10\) L. (a) Write the ODE system. (b) Find \(\alpha\) and \(\beta\) as eigenvalues of the system matrix. (c) Determine the coefficients \(A\) and \(B\) in \(C_1(t) = Ae^{-\alpha t}+Be^{-\beta t}\). (d) Compute \(\mathrm{AUC}_{0\to\infty}\) and \(\mathrm{AUC}_{0\to 24}\). (e) Compute the clearance and verify it equals \(k_{10}V_1\).

Problem 13. Bayes’ theorem and diagnostic testing

A screening test has sensitivity 95% and specificity 90%. (a) Compute the PPV at 1% prevalence. (b) Compute the NPV and distinguish it from \(P(D\mid -)\). (c) Recompute the PPV at 10% prevalence. (d) What sensitivity and specificity would be needed to reach PPV \(= 0.5\) at 1% prevalence, holding one fixed?

Problem 14. Maximum likelihood for a Poisson rate

Forty patients are imaged and the counts in a region of interest have sample mean \(\bar{k} = 128.4\). Assuming a Poisson model, (a) write the log-likelihood, (b) derive \(\hat\lambda_{\mathrm{MLE}}\), (c) compute its standard error from the Fisher information, and (d) give a 95% confidence interval.

Problem 15. Error propagation for the volume of a sphere

A sphere’s radius is measured as \(r = 2.0 \pm 0.1\) cm. (a) Use Gaussian error propagation to estimate \(\sigma_V\) for \(V = \frac{4}{3}\pi r^3\). (b) Compare with a Monte Carlo simulation of \(10^5\) draws from \(r\sim\mathcal{N}(2.0, 0.1^2)\). (c) Express the result as a relative uncertainty and explain the factor of 3.

Problem 16. Confidence interval for a proportion

In a trial, 73 of 200 patients respond. (a) Compute a 95% Wald interval. (b) Compute the Wilson interval. (c) Explain which is preferable and why. (d) How many patients would be needed to halve the interval width?

Problem 17. Eigenvalues, eigenvectors, and PCA

Generate \(n = 200\) points from a bivariate normal with \(\boldsymbol{\mu} = (0,0)\) and \(\Sigma = \begin{bmatrix} 2.0 & 1.5\\ 1.5 & 2.0\end{bmatrix}\). (a) Compute the sample covariance. (b) Find its eigenvalues and eigenvectors. (c) Plot the data with the principal axes overlaid. (d) Report the variance explained by PC1. (e) Verify that the SVD of the centred data gives the same result.

Problem 18. SVD denoising

Construct a \(64\times64\) rank-5 matrix, add Gaussian noise with \(\sigma = 0.12\), and denoise by SVD truncation. (a) Plot the singular-value spectrum. (b) Plot reconstruction RMSE against retained rank. (c) Report the optimal rank and explain why it may differ from 5. (d) Explain why the same procedure applied to a 1-D signal fails, and what would be done instead.

Problem 19. Diffusion tensor metrics

A voxel has diffusion tensor eigenvalues \((1.7, 0.2, 0.2)\times10^{-3}\) mm\(^2\)/s. (a) Compute MD and FA. (b) Repeat for \((0.7, 0.7, 0.7)\times10^{-3}\). (c) Explain why FA is rotation-invariant and why that matters clinically. (d) An acute stroke reduces all three eigenvalues by 40%. What happens to MD and to FA?

Problem 20. Matching the tool to the problem

For each task, name the most appropriate tool from this chapter and the underlying equation or principle.

1. Compute the volume of a tumour from a stack of segmented cross-sections.
2. Estimate a drug’s clearance from serial plasma concentrations.
3. Decide whether a perfusion change in an fMRI region is statistically significant.
4. Reduce a 1000-dimensional radiomic dataset to two dimensions for a figure.
5. Find the oscillation period of a damped harmonic oscillator from its equation of motion.
6. Determine whether two image regions have the same intensity distribution.
7. Predict how much longer a PET scan must run to halve the image noise.
8. Explain why an ultrasound Doppler display shows flow reversal in a high-velocity jet.

Solutions Key

Solution 1

\(f'(x) = 3x^2 - 12x + 9 = 3(x-1)(x-3)\), so the critical points are \(x = 1\) and \(x = 3\). Since \(f''(x) = 6x - 12\), we have \(f''(1) = -6 < 0\) (local maximum, \(f(1) = 5\)) and \(f''(3) = +6 > 0\) (local minimum, \(f(3) = 1\)). There is an inflection point at \(x = 2\), where \(f''\) vanishes and changes sign, with \(f(2) = 3\). It is not a critical point because \(f'(2) = -3 \ne 0\): the curvature reverses there, but the tangent is not horizontal.

Solution 2

With \(u = x\) and \(dv = e^{-x}dx\), so \(du = dx\) and \(v = -e^{-x}\):

\[\int_0^1 x e^{-x}dx = \left[-xe^{-x}\right]_0^1 + \int_0^1 e^{-x}dx = -e^{-1} + \left(1 - e^{-1}\right) = 1 - \frac{2}{e} \approx 0.2642.\]

In R, integrate(function(x) x*exp(-x), 0, 1)$value returns \(0.2642411\), agreeing to seven decimals.

Solution 3

Errors fall by a factor of about 2 per doubling of \(n\) for the left rule (\(O(h)\), first order) and by about 4 for the midpoint and trapezoid rules (\(O(h^2)\), second order). The convergence order is estimated as \(\log_2(E_n/E_{2n})\). The midpoint rule has roughly half the error of the trapezoid rule at equal \(n\), with opposite sign, which is what Simpson’s rule exploits by combining them.

Solution 4

\(\det(A) = -3\) and \(\mathrm{tr}(A) = 16\).

\[A^{-1} = \begin{bmatrix} -\tfrac{2}{3} & -\tfrac{4}{3} & 1\\[2pt] -\tfrac{2}{3} & \tfrac{11}{3} & -2\\[2pt] 1 & -2 & 1\end{bmatrix}, \qquad A^{T} = \begin{bmatrix} 1&4&7\\ 2&5&8\\ 3&6&10\end{bmatrix}.\]

Then \(\mathbf{x} = A^{-1}\mathbf{b} = \left(-\tfrac{1}{3},\ \tfrac{2}{3},\ 0\right)^{T}\), and substituting back gives \(A\mathbf{x} = (1,2,3)^{T}\). The eigenvalues sum to \(16\) and multiply to \(-3\), confirming Eq. (95).

Solution 5

1. \((1+i)^2 = 2i\), so \((1+i)^4 = (2i)^2 = -4\).
2. \(\lvert 1+i\rvert = \sqrt{2}\) and \(\arg(1+i) = \pi/4\), so \((1+i)^4 = (\sqrt{2})^4 e^{i\pi} = 4e^{i\pi} = -4\). Both give \(-4\).

Solution 6

1. \(\tilde{I}_1 = 5(\cos\tfrac{\pi}{6} + i\sin\tfrac{\pi}{6}) \approx 4.3301 + 2.5000i\) and \(\tilde{I}_2 \approx 2.1213 - 2.1213i\), so \(\tilde{I} \approx 6.4514 + 0.3787i\).
2. \(\lvert\tilde{I}\rvert \approx 6.463\) A and \(\arg\tilde{I}\approx 0.0586\) rad (\(3.36^\circ\)).
3. \(I(t)\approx 6.463\cos(\omega t + 0.0586)\) A.
4. The two phasors are \(75^\circ\) apart, so their imaginary parts largely cancel; only co-phased components add arithmetically.

Solution 7

\(y(t) = t\) for \(0\le t\le 1\), \(y(t) = 2-t\) for \(1\le t\le 2\), and zero elsewhere: a triangle peaking at \(t = 1\) with height 1. The support follows from the general rule that convolving supports of length \(L_1\) and \(L_2\) gives support of length \(L_1+L_2\), here \(1+1 = 2\). The area is \(\int y = 1 = (\int x)(\int h)\), consistent with \(\hat{y}(0) = \hat{x}(0)\hat{h}(0)\).

Solution 8

\[\hat{f}(\omega) = \int_{-\infty}^{\infty}e^{-\pi t^2}e^{-i2\pi\omega t}dt = \int_{-\infty}^{\infty}\exp\!\left[-\pi\left(t^2 + 2i\omega t\right)\right]dt.\]

Completing the square, \(t^2 + 2i\omega t = (t+i\omega)^2 + \omega^2\), so

\[\hat{f}(\omega) = e^{-\pi\omega^2}\int_{-\infty}^{\infty}e^{-\pi(t+i\omega)^2}dt = e^{-\pi\omega^2},\]

using \(\int e^{-\pi u^2}du = 1\) after a contour shift. The Gaussian is its own transform. Consequently, a Gaussian pulse of temporal width \(\sigma_t\) has spectral width \(\sigma_\omega \propto 1/\sigma_t\): short pulses are necessarily broadband, which is why high axial resolution in ultrasound and OCT requires wide bandwidth.

Solution 9

1. \(f_{\mathrm{Nyq}} = 500/2 = 250\) Hz.
2. Only marginally. The Nyquist frequency exactly equals the highest signal frequency, so any imperfect analog roll-off or out-of-band noise aliases into the diagnostic band. Equation (50) requires strict inequality.
3. \(\lvert 300 - 500\rvert = 200\) Hz, squarely inside the diagnostic band and indistinguishable from genuine 200 Hz content.
4. Sample at 1000 Hz or more (four times the highest signal frequency) with an analog low-pass filter whose passband ends near 250 Hz and whose stopband begins by about 500 Hz, giving a 250 Hz guard band.

Solution 10

1. \(\int_{-\infty}^{\infty}C(x,t)dx = \frac{A\sigma_0}{\sqrt{\sigma_0^2+Dt}}\cdot\sqrt{4\pi(\sigma_0^2+Dt)} = 2A\sigma_0\sqrt{\pi}\), independent of \(t\): mass is conserved.
2. Since \(\exp\!\big(-x^2/(4w^2)\big)\) is a Gaussian of variance \(2w^2\) with \(w^2 = \sigma_0^2 + Dt\), we get \(\langle x^2\rangle = 2(\sigma_0^2 + Dt)\), growing linearly in \(t\).
3. Doubling the width means \(w = 2\sigma_0\), hence \(Dt = 3\sigma_0^2\) and \(D = 3(10^{-6}\,\text{m})^2/(7.1\times10^{-3}\,\text{s}) = 4.2\times10^{-10}\) m\(^2\)/s, in good agreement with the literature value of about \(4.25\times10^{-10}\) m\(^2\)/s for fluorescein in water near room temperature.

Solution 11

1. \(N(t) = 10^6 e^{0.08t}\).
2. \(N(30) = 10^6 e^{2.4} \approx 1.10\times10^7\) cells.
3. \(t_d = \ln 2/0.08 \approx 8.66\) days.
4. \(N(30) = 10^{10}/\big(1 + (10^4-1)e^{-2.4}\big)\approx 1.101\times10^{7}\), only 0.1% below the exponential value. At \(10^7\) cells the tumour occupies \(10^{-3}\) of its carrying capacity, so the logistic brake \((1 - N/K)\approx 0.999\) is essentially inactive. Distinguishing the models requires observation much closer to \(K\).

Solution 12

1. \(\dot{C}_1 = -(k_{10}+k_{12})C_1 + k_{21}C_2 = -0.6C_1 + 0.3C_2\) and \(\dot{C}_2 = 0.4C_1 - 0.3C_2\), with \(C_1(0) = 100/10 = 10\) mg/L and \(C_2(0) = 0\).
2. The system matrix \(\begin{bmatrix}-0.6 & 0.3\\ 0.4 & -0.3\end{bmatrix}\) has characteristic equation \(\lambda^2 + 0.9\lambda + 0.06 = 0\), giving \(\lambda = -0.8275\) and \(-0.0725\), so \(\alpha = 0.8275\) hr\(^{-1}\) (\(t_{1/2} = 0.84\) h) and \(\beta = 0.0725\) hr\(^{-1}\) (\(t_{1/2} = 9.56\) h).
3. \(A = C_0(\alpha - k_{21})/(\alpha-\beta) = 6.987\) and \(B = C_0(k_{21}-\beta)/(\alpha-\beta) = 3.013\) mg/L, and \(A + B = 10 = C_0\) as required.
4. \(\mathrm{AUC}_{0\to\infty} = A/\alpha + B/\beta = 8.443 + 41.557 = 50.0\) mg·h/L; \(\mathrm{AUC}_{0\to24} = 42.71\) mg·h/L, so 15% of total exposure occurs after 24 h.
5. \(CL = \mathrm{Dose}/\mathrm{AUC}_{0\to\infty} = 100/50 = 2.0\) L/h \(= k_{10}V_1 = 0.2\times10\). ✓

Solution 13

1. \(\mathrm{PPV} = \dfrac{0.95\times0.01}{0.95\times0.01 + 0.10\times0.99} = \dfrac{0.0095}{0.1085} = 0.0876\), about 8.8%.
2. \(\mathrm{NPV} = P(\bar{D}\mid-) = \dfrac{0.90\times0.99}{0.90\times0.99 + 0.05\times0.01} = 0.99944\). The complementary quantity \(P(D\mid-) = 1-\mathrm{NPV} = 0.00056\) (0.056%) is the false omission rate; the two are often confused.
3. At 10% prevalence, \(\mathrm{PPV} = 0.095/0.185 = 0.514\).
4. Holding sensitivity at 0.95, PPV \(= 0.5\) at 1% prevalence requires \(0.95\times0.01 = (1-\mathrm{spec})\times0.99\), so specificity \(= 1 - 0.0096 = 0.9904\). Raising specificity from 90% to 99% is what buys a usable screening test; raising sensitivity alone cannot achieve it, since even perfect sensitivity gives PPV \(= 0.01/(0.01+0.099) = 0.092\).

Solution 14

1. \(\ell(\lambda) = \sum_{i=1}^{40}\left[k_i\log\lambda - \lambda - \log k_i!\right] = 40\bar{k}\log\lambda - 40\lambda - \sum\log k_i!\).
2. \(\partial\ell/\partial\lambda = 40\bar{k}/\lambda - 40 = 0 \Rightarrow \hat\lambda_{\mathrm{MLE}} = \bar{k} = 128.4\).
3. The Fisher information is \(I(\lambda) = n/\lambda\), so \(\mathrm{SE}(\hat\lambda) = \sqrt{\lambda/n} = \sqrt{128.4/40} = 1.79\).
4. \(128.4 \pm 1.96\times1.79 = (124.9,\ 131.9)\).

Solution 15

1. \(\partial V/\partial r = 4\pi r^2\), so \(\sigma_V = 4\pi(2.0)^2(0.1) = 1.6\pi \approx 5.03\) cm\(^3\).
2. Monte Carlo with \(10^5\) draws gives \(\sigma_V \approx 5.06\) cm\(^3\); the small excess reflects the third-order nonlinearity that the first-order expansion omits.
3. \(V = 33.5\) cm\(^3\), so the relative uncertainty is \(5.03/33.5 = 15\%\), exactly three times the 5% relative uncertainty in \(r\). This is the power rule \(\sigma_f/f = n\,\sigma_x/x\) with \(n = 3\): volume errors are amplified threefold, which is why segmentation-boundary uncertainty dominates tumour-volume uncertainty in Chapter 7.

Solution 16

1. \(\hat{p} = 73/200 = 0.365\); \(\mathrm{SE} = \sqrt{0.365\times0.635/200} = 0.0340\); the Wald interval is \(0.365 \pm 0.0667 = (0.298,\ 0.432)\).
2. The Wilson interval is \((0.301,\ 0.434)\), obtained from \(\left(\hat{p} + \tfrac{z^2}{2n}\right)\big/\left(1+\tfrac{z^2}{n}\right) \pm \tfrac{z}{1+z^2/n}\sqrt{\tfrac{\hat p(1-\hat p)}{n} + \tfrac{z^2}{4n^2}}\).
3. Wilson is preferable: it has better coverage at small \(n\) or extreme \(\hat{p}\), never extends outside \([0,1]\), and is what prop.test() reports. The two nearly coincide here because \(n\) is large and \(\hat{p}\) is near 0.5.
4. Width scales as \(n^{-1/2}\), so halving it requires quadrupling \(n\) to 800.

Solution 17

1. For \(n=200\) the sample covariance is close to \(\Sigma\).
2. Theoretically \(\lambda_1 = 2.0+1.5 = 3.5\) and \(\lambda_2 = 2.0-1.5 = 0.5\), with eigenvectors \((1,1)/\sqrt{2}\) and \((1,-1)/\sqrt{2}\).
3. The scatter is elongated along the first principal axis at \(45^\circ\).
4. PC1 explains \(3.5/(3.5+0.5) = 87.5\%\) of the variance.
5. For centred \(X\), \(\mathrm{svd}(X)\) gives \(d_k^2/(n-1) = \lambda_k\) and right singular vectors equal to the eigenvectors of \(\Sigma\) (up to sign).

Solution 18

1. The spectrum shows a few large singular values above a broad noise floor.
2. RMSE against the clean image falls, reaches a minimum, then rises as noise-dominated components are readmitted.
3. The optimal rank is often 4 rather than 5. The fifth component’s singular value sits near the noise floor, so retaining it adds more variance than the bias it removes. This is the bias–variance trade-off of Section 1.16 expressed in matrix form; halving the noise would lower the floor and typically restore the optimum to 5.
4. A length-\(N\) vector, viewed as an \(N\times1\) matrix, has exactly one non-zero singular value equal to its norm, so truncation does nothing. The standard remedy is singular spectrum analysis: embed the signal in a Hankel (trajectory) matrix of lagged copies, truncate the SVD of that matrix, and average the anti-diagonals to recover a denoised series.

Solution 19

1. \(\mathrm{MD} = (1.7+0.2+0.2)/3 = 0.7\times10^{-3}\) mm\(^2\)/s; with \(\sum\lambda_i^2 = (2.89+0.04+0.04)\times10^{-6}\) and \(\sum(\lambda_i-\mathrm{MD})^2 = (1.0+0.25+0.25)\times10^{-6}\), \(\mathrm{FA} = \sqrt{1.5\times1.5/2.97} = 0.870\).
2. For \((0.7,0.7,0.7)\), MD \(= 0.7\times10^{-3}\) and FA \(= 0\) exactly, since all eigenvalues equal the mean.
3. FA and MD are functions of the eigenvalues alone, which are invariant under rotation of the coordinate frame. Clinically this means the measured anisotropy does not depend on how the patient was positioned in the scanner, so values are comparable across sessions and sites.
4. Scaling all eigenvalues by 0.6 scales MD by 0.6 (to \(0.42\times10^{-3}\)), but FA is a ratio of like-degree quantities and is therefore unchanged. This is exactly the observation in acute stroke: diffusivity falls sharply while anisotropy is initially preserved.

Solution 20

1. Integral calculus, the disk method \(V = \int \pi r(z)^2 dz\) (Section 1.1.3), applied slice by slice.
2. First-order ODE plus regression: fit \(C(t) = C_0e^{-k_{el}t}\) and use \(CL = k_{el}V_d\), or model-independently \(CL = \mathrm{Dose}/\mathrm{AUC}\) (Eq. 1.73).
3. Hypothesis testing, a \(t\)-test or a general linear model against a null distribution of no activation, with multiple-comparison control (Section 1.16.5).
4. PCA if a linear summary suffices, t-SNE or UMAP for visualization only, remembering that inter-cluster distances in the latter are not interpretable (Section 1.18).
5. Complex roots of the characteristic equation: for \(m\ddot{x}+c\dot{x}+kx = 0\) the roots of \(mr^2+cr+k = 0\) are complex in the underdamped case, and the damped period is \(2\pi/\omega_d\) with \(\omega_d = \sqrt{k/m - c^2/(4m^2)}\) (Sections 1.6, 1.13.6).
6. Two-sample \(t\)-test for means, or the Kolmogorov–Smirnov test to compare entire distributions, which is the correct choice when the histograms differ in shape rather than location (Section 1.16.5).
7. Poisson statistics and the square-root law: halving the noise requires quadrupling the counts, so the scan must run four times as long at constant activity, subject to correction for tracer decay over the extended acquisition (Section 1.15).
8. Nyquist criterion: the pulse repetition frequency sets a Nyquist velocity, and flow exceeding it aliases to a negative apparent velocity (Section 1.10).

References and Further Reading

1. Kreyszig, E. (2011). Advanced Engineering Mathematics, 10th ed. Wiley.
2. Strang, G. (2016). Introduction to Linear Algebra, 5th ed. Wellesley–Cambridge Press.
3. Boyce, W. E., DiPrima, R. C., and Meade, D. B. (2017). Elementary Differential Equations and Boundary Value Problems, 11th ed. Wiley.
4. Oppenheim, A. V. and Schafer, R. W. (2009). Discrete-Time Signal Processing, 3rd ed. Pearson.
5. Bracewell, R. N. (2000). The Fourier Transform and Its Applications, 3rd ed. McGraw-Hill.
6. Kak, A. C. and Slaney, M. (2001). Principles of Computerized Tomographic Imaging. SIAM.
7. Nishimura, D. G. (2010). Principles of Magnetic Resonance Imaging. Stanford University.
8. Bushberg, J. T., Seibert, J. A., Leidholdt, E. M., and Boone, J. M. (2011). The Essential Physics of Medical Imaging, 3rd ed. Lippincott Williams & Wilkins.
9. Barrett, H. H. and Myers, K. J. (2004). Foundations of Image Science. Wiley.
10. Gudbjartsson, H. and Patz, S. (1995). The Rician distribution of noisy MRI data. Magnetic Resonance in Medicine, 34(6), 910–914.
11. Taylor, J. R. (1997). An Introduction to Error Analysis, 2nd ed. University Science Books.
12. Casella, G. and Berger, R. L. (2002). Statistical Inference, 2nd ed. Duxbury.
13. Hastie, T., Tibshirani, R., and Friedman, J. (2009). The Elements of Statistical Learning, 2nd ed. Springer.
14. Rowland, M. and Tozer, T. N. (2010). Clinical Pharmacokinetics and Pharmacodynamics, 4th ed. Lippincott Williams & Wilkins.
15. McInnes, L., Healy, J., and Melville, J. (2018). UMAP: Uniform Manifold Approximation and Projection for dimension reduction. arXiv:1802.03426.
16. Dinov, I. D. (2021). Data Science: Time Complexity, Inferential Uncertainty, and Spacekime Analytics. De Gruyter.
17. Dinov, I. D. and Velev, M. V. (2021). Data Science: Time Complexity, Inferential Uncertainty, and Spacekime Analytics. https://tciu.predictive.space.
18. SOCR Probability and Statistics EBook.
19. DSPA2 Chapter 4: Dimensionality Reduction.
20. DSPA Chapter 4: Linear Algebra and Matrix Computing.
21. DSPA2 Chapter 13: Function Optimization.
22. SOCR Probability Distributome Project.
23. SOCR Statistical Data Analyzer (SDA).

SOCR/BPAD Resource Visitor number