2 Mathematical Preliminaries
The daily life of industrial vision engineering looks like tuning parameters, configuring lighting, and writing pipelines, but underneath it the same cast of mathematical characters is always on stage: “straightening up” a part relies on matrix multiplication, turning a string of edge points into a datum line relies on least squares, deciding whether a gray value counts as anomalous relies on probability distributions, and pushing localization accuracy from whole pixels to subpixels relies on interpolation. In university curricula this material is scattered across three courses — linear algebra, probability and statistics, and numerical analysis. This chapter gathers the small slice that genuinely keeps reappearing on the production line, selected by one criterion: which later chapter will need it. We aim not for completeness, but for every formula to have a definite destination. The table below is the learning map for this chapter:
| Section | Content | Later chapters it supports |
|---|---|---|
| Section 2.1 | Vectors, transformation matrices, homogeneous coordinates, eigenvectors | Chapter 10, Chapter 5, Chapter 14 |
| Section 2.2 | Normal equations, weighted least squares, numerical stability | Chapter 14, Chapter 5 |
| Section 2.3 | Expectation and variance, the Gaussian distribution, histograms, the central limit theorem | Chapter 6, Chapter 7, Chapter 26 |
| Section 2.4 | Bilinear interpolation, three-point parabolic interpolation | Chapter 10, Chapter 20, Chapter 14 |
A note up front: this is a theory chapter, with no standalone companion code project; every conclusion derived here will be cashed in repeatedly in the real experiments of later chapters — for example, Chapter 6 uses Gaussian-noise experiments with \(\sigma=18\) to confirm that “averaging cancels independent errors”, and Chapter 14 uses an inlier-RMSE comparison of 0.723 px versus 0.257 px to demonstrate the squared penalty’s sensitivity to outliers. Readers can treat this chapter as a “memo with proofs”: whenever a formula in a later chapter leaves you unsure of its origin, flip back here.
2.1 Vectors, Matrices, and Geometric Transformations
Throughout the book, bold lowercase letters denote vectors, always column vectors: a continuous coordinate point in the image plane is written \(\mathbf{x} = (x, y)^\top\), consistent with the Notation chapter; in a discrete image \(f[n,m]\), \(n\) is the row index (downward) and \(m\) is the column index (rightward). Matrices are denoted by uppercase letters. A \(2\times 2\) matrix \(A\) acting on a point (\(\mathbf{x}' = A\mathbf{x}\)) is a linear transformation: rotation, scaling, shear — each is some specific \(A\).
The three 2D transformations industrial vision uses most form a family ordered from strict to permissive. A rigid transformation contains only rotation and translation:
\[ \mathbf{x}' = R\,\mathbf{x} + \mathbf{t}, \qquad R = \begin{pmatrix} \cos\theta & -\sin\theta \\ \sin\theta & \cos\theta \end{pmatrix}, \]
with 3 degrees of freedom (the angle \(\theta\) plus the two components of the translation \(\mathbf{t}\)), and preserves all lengths and angles — exactly the mathematical statement of “the part shifted and rotated on the stage, but the part itself did not change”, and the standard model for position correction. A similarity transformation adds a global scale factor \(s\) on top of the rigid one (that is, \(\mathbf{x}' = sR\,\mathbf{x} + \mathbf{t}\), 4 degrees of freedom): lengths are no longer preserved, but angles and shape still are — suitable when the working distance changes and the image grows or shrinks. An affine transformation opens this up further to an arbitrary invertible \(2\times 2\) matrix plus translation (6 degrees of freedom), allowing independent scaling along the two axes and shear; straight lines remain straight and parallel lines remain parallel, but angles are no longer preserved. Chapter 10 will put this family to work for image resampling and position correction.
Note that all three expressions above are written as “matrix multiplication plus translation” — the translation cannot be folded into the \(2\times 2\) matrix, because a linear transformation must map the origin to the origin. Homogeneous coordinates solve this with a simple dimension lift: append a 1 to every point, writing \(\tilde{\mathbf{x}} = (x, y, 1)^\top\), and the entire affine transformation becomes a single \(3\times 3\) matrix multiplication:
\[ \begin{pmatrix} x' \\ y' \\ 1 \end{pmatrix} = \begin{pmatrix} a_{11} & a_{12} & t_x \\ a_{21} & a_{22} & t_y \\ 0 & 0 & 1 \end{pmatrix} \begin{pmatrix} x \\ y \\ 1 \end{pmatrix}. \]
The engineering value of homogeneous coordinates lies in chaining: image coordinates → distortion correction → position correction → physical coordinates — each step is a \(3\times 3\) matrix, and the whole chain can be pre-multiplied into one matrix, so at run time each point costs a single multiplication. If the last row is not restricted to \((0,0,1)\), we obtain the perspective transformation (homography) — Chapter 5 cannot do without it when describing how a planar calibration target is imaged.
Finally, the geometric intuition of the eigenvector, which will take center stage in the closed-form line-fitting solution of Chapter 14. For a set of points \(\{\mathbf{p}_i\}_{i=1}^{N}\), write the centroid as \(\bar{\mathbf{p}} = \frac{1}{N}\sum_i \mathbf{p}_i\) and define the scatter matrix
\[ S = \sum_{i=1}^{N} (\mathbf{p}_i - \bar{\mathbf{p}})(\mathbf{p}_i - \bar{\mathbf{p}})^\top, \]
a \(2\times 2\) symmetric matrix. Its two eigenvectors are mutually perpendicular, and the one whose eigenvalue is larger points in the direction along which the point set is most spread out — feed in a string of edge points distributed along some edge, and the principal eigenvector hands you the edge’s direction. This is no coincidence: one can prove that the line minimizing the sum of squared perpendicular distances to the points passes exactly through the centroid, with its direction given by the principal eigenvector of \(S\). Compute one centroid, solve one \(2\times 2\) eigenproblem, and line fitting has a closed-form solution — this is the entire content of the phrase “the solution is closed-form” in Chapter 14.
2.2 Least Squares and Robustness
Least squares is the engine of measurement-class algorithms, and the simplest case is worth deriving in full. Given \(N\) points \((x_i, y_i)\), we want to fit the line \(y = kx + b\), defining “fits well” as minimizing the sum of squared vertical residuals:
\[ E(k, b) = \sum_{i=1}^{N} (k x_i + b - y_i)^2 . \]
\(E\) is a smooth convex function of \(k, b\), so at the minimum the partial derivatives vanish:
\[ \frac{\partial E}{\partial k} = 2\sum_i x_i (k x_i + b - y_i) = 0, \qquad \frac{\partial E}{\partial b} = 2\sum_i (k x_i + b - y_i) = 0 . \]
Rearranged into a linear system in \((k,b)\), these are the normal equations:
\[ \begin{pmatrix} \sum_i x_i^2 & \sum_i x_i \\ \sum_i x_i & N \end{pmatrix} \begin{pmatrix} k \\ b \end{pmatrix} = \begin{pmatrix} \sum_i x_i y_i \\ \sum_i y_i \end{pmatrix}. \]
Two rows, two columns — one solve yields \(k, b\); and the second row throws in a property worth remembering: \(b = \bar{y} - k\bar{x}\), that is, the least-squares line always passes through the centroid of the point set. The general case is completely isomorphic: write the model as \(A\boldsymbol{\theta} \approx \mathbf{y}\) (each row of \(A\) is one observation, \(\boldsymbol{\theta}\) the parameters to be estimated), and the normal equations are \(A^\top A\,\boldsymbol{\theta} = A^\top \mathbf{y}\). Whether it is circle fitting, homography estimation for a calibration target, or solving the hand–eye relation, most problems eventually land in this form.
The condition number measures how much a solution amplifies perturbations of the input, and the condition number of \(A^\top A\) is the square of that of \(A\). If you fit directly in pixel coordinates (easily in the thousands), \(\sum x_i^2\) and \(N\) differ by six orders of magnitude, the matrix is nearly ill-conditioned, and floating-point errors are amplified sharply; circle fitting contains an \(x^2+y^2\) term, making things even worse. The engineering rule: first subtract the centroid to center the data (and, if necessary, divide by a scale to normalize), solve in small coordinates, then transform the result back — a few lines of code buy back several significant digits.
Two common extensions. First, weighted least squares: if the observations differ in trustworthiness, give each term a weight \(w_i \ge 0\) and minimize \(\sum_i w_i r_i^2\); the normal equations become \(A^\top W A\,\boldsymbol{\theta} = A^\top W \mathbf{y}\) (\(W\) the diagonal weight matrix) — the derivation is word-for-word the same as above, only with an extra \(w_i\) inside every sum. Second, replace the “vertical residual” with the perpendicular distance from point to line: the optimum is then no longer given by the normal equations, but is exactly the closed form of the previous section — through the centroid, direction along the principal eigenvector of the scatter matrix; it favors no coordinate axis and is the right way to do geometric fitting.
The Achilles’ heel of least squares hides in that same square: a single spurious point with a 12 px residual contributes nearly 600 times as much to the objective as a normal point with a 0.5 px residual, and the fit can be hijacked by a tiny minority of outliers — the controlled experiment in Chapter 14 puts striking numbers on this (inlier RMSE 0.723 px versus 0.257 px). The cure is to replace the squared penalty with one more forgiving of large residuals (such as Huber weighting (Huber 1964), solved iteratively via weighted least squares) or with random sample consensus, RANSAC (Fischler and Bolles 1981); both are left for that chapter to develop alongside the experiments. Here it suffices to remember the root cause: squaring amplifies the voice of large residuals.
2.3 Probability and Statistics Essentials
Noise cannot be predicted pixel by pixel, but it can be characterized statistically. Treat a single gray-value observation as a random variable \(X\); its two most important summaries are the expectation \(\mu = \mathrm{E}[X]\) (what repeated observations average toward) and the variance \(\sigma^2 = \mathrm{E}[(X-\mu)^2]\) (how spread out values are around the expectation). The square root of the variance, \(\sigma\), is the standard deviation; it carries the same units as the gray values, and in practice serves directly as the “noise amplitude”. Together, the linearity of expectation and the additivity of variance for independent observations immediately yield a law used every day on the production line: averaging \(N\) independent, identically distributed observations gives a mean whose variance is \(1/N\) of a single observation’s, i.e. the standard deviation shrinks to \(\sigma/\sqrt{N}\) — multi-frame averaging, the neighborhood smoothing of Chapter 6, and the projection averaging along search lines in Chapter 14 are all incarnations of this \(\sqrt{N}\) law.
The most important distribution in machine vision is the Gaussian distribution \(\mathcal{N}(\mu, \sigma^2)\), with probability density
\[ p(x) = \frac{1}{\sqrt{2\pi}\,\sigma}\exp\!\left(-\frac{(x-\mu)^2}{2\sigma^2}\right). \]
Its values are highly concentrated: the probabilities of falling within \(\mu \pm \sigma\), \(\mu \pm 2\sigma\), and \(\mu \pm 3\sigma\) are about 68.3%, 95.4%, and 99.7% respectively — the last is the famous 3σ rule. It gives threshold setting a direct vocabulary: if the gray values of a good region follow \(\mathcal{N}(\mu, \sigma^2)\), then flagging “outside \(\mu \pm 3\sigma\)” as anomalous yields a per-pixel false-alarm rate of only about 0.3%. But beware the effect of scale: in a five-megapixel image, 0.3% means over ten thousand pixels cross the line out of sheer luck, which is why defect-detection practice often loosens the gate to \(4\sigma\)–\(6\sigma\), or follows the 3σ test with a connected-area filter — Chapter 26 will return to this accounting.
Estimating \(\sigma\) itself also needs protection against outliers: the sample standard deviation contains a squared term, and a few bad points can inflate it. The robust alternative is the median absolute deviation (MAD): \(\hat{\sigma} = 1.4826 \times \operatorname{median}_i\big(|x_i - \operatorname{median}(x)|\big)\), where the factor 1.4826 makes it agree with the standard deviation on Gaussian data. The median naturally ignores a small number of extreme values, which is why MAD is a standard fixture of threshold statistics on the line.
Where does the distribution of a real image come from? Just count: tally how many pixels take each gray level and you get the histogram; divide by the total pixel count and it becomes the empirical estimate of the gray-value distribution — expectation and variance can both be computed straight from it. An image of a dark part on a bright background shows a bimodal histogram, and the valley between the two peaks is a natural segmentation threshold; Otsu’s method in Chapter 7 treats the histogram precisely as a probability distribution and searches for the threshold maximizing the between-class variance — the expectation and variance of this section are its entire raw material.
Finally, an answer to a question that has hung over two chapters: why model sensor noise as Gaussian? The answer is the central limit theorem: the sum of a large number of mutually independent, individually tiny random perturbations tends toward a Gaussian distribution, regardless of what each perturbation’s own distribution is. Pixel readout perturbations are exactly this situation — thermal fluctuations of dark current, readout-circuit noise, and electronic noise from gain amplification stack layer upon layer, each small and independent, so the total effect is approximately Gaussian. This is the justification for Chapter 6 running all of its controlled experiments with Gaussian noise at \(\sigma=18\); it also explains that chapter’s counterexample — salt-and-pepper noise comes from single, large-amplitude events such as dead pixels and transmission errors, violating the “many tiny independent perturbations” premise, so it is not Gaussian, and methods built for Gaussian noise should not be used against it.
2.4 Interpolation
An image has values only at integer grid points, yet algorithms keep demanding values at non-integer positions: rotate an image, and a target pixel mapped back into the source almost never lands on a grid point (Chapter 10); take a gray profile along a search line in an arbitrary direction, and the sample points fall between grid points too (Chapter 20). Interpolation is responsible for answering “what lies between the grid points”.
The cheapest answer is nearest-neighbor interpolation: take the value of the closest grid point. It creates no new gray values and is the fastest, but positionally it amounts to rounding by up to 0.5 pixel; edges come out jagged after geometric transformations, so it is essentially never used for measurement. The standard answer is bilinear interpolation: let the sample point fall inside the unit square bounded by rows \(n\), \(n+1\) and columns \(m\), \(m+1\), with fractional offsets \(\beta, \alpha \in [0,1)\) in the row and column directions; then
\[ f(n+\beta,\, m+\alpha) = (1-\alpha)(1-\beta)\, f[n,m] + \alpha(1-\beta)\, f[n,m+1] + (1-\alpha)\beta\, f[n+1,m] + \alpha\beta\, f[n+1,m+1]. \]
The four weights are exactly the four diagonally opposite areas into which the sample point divides the unit square: the closer the point is to a grid point, the larger that grid point’s say; the four weights always sum to 1, so a flat region interpolates to a flat region. Bilinear output is continuous and cheap to compute, making it the default choice for industrial geometric transformations; the higher-order bicubic interpolation trades a \(4\times 4\) neighborhood for a smoother result, and Chapter 10 will compare them on real measurements.
The other kind of interpolation serves not resampling but finding the precise location of an extremum — the mathematical core of subpixel localization. A discrete extremum of a one-dimensional profile is accurate only to the whole pixel: place the extremum at local coordinate \(u=0\) with value \(g_0\), and its left and right neighbors at \(u=\mp 1\) with values \(g_{-1}, g_{+1}\). The true continuous extremum need not sit at \(u=0\), but three points suffice to determine a parabola \(g(u) = au^2 + bu + c\), whose vertex approximates the true extremum. Substituting the three samples:
\[ g(0) = c = g_0, \qquad g(\pm 1) = a \pm b + c = g_{\pm 1}, \]
and adding and subtracting the two equations gives immediately
\[ a = \frac{g_{+1} + g_{-1} - 2g_0}{2}, \qquad b = \frac{g_{+1} - g_{-1}}{2}. \]
The parabola’s vertex lies where the derivative vanishes, at \(u^\ast = -\,b/(2a)\); substituting \(a, b\):
\[ \delta = u^\ast = \frac{g_{-1} - g_{+1}}{2\,(g_{-1} - 2g_0 + g_{+1})}. \]
This is precisely the subpixel formula of step 4 in the edge-localization procedure of Chapter 14. Two sanity checks deepen one’s trust in it. First, if \(g_{-1} = g_{+1}\) (left–right symmetric), then \(\delta = 0\) and the extremum sits exactly on the whole pixel, as intuition demands. Second, as long as \(g_0\) is a strict discrete maximum, the denominator satisfies \(g_{-1} - 2g_0 + g_{+1} < 0\) without exception and \(|\delta| \le 1/2\) — the vertex never escapes past the midpoint to a neighboring pixel, so the formula is numerically safe. Three samples and one division push localization accuracy from the 1-pixel scale to the 0.1-pixel scale — that is the entire cost behind “subpixel accuracy is almost free”.
Industry Case: Three Small Things That Make Measurements Right
A vision engineer at a precision structural-parts plant once summed up three “small math” lessons at a postmortem meeting. First: an early roundness-inspection program fitted circles directly in full-frame pixel coordinates — coordinate values in the thousands, and with the \(x^2+y^2\) term the normal equations’ condition number was enormous; measuring the same part ten times in a row, the fitted center drifted by 0.3 px. Subtracting the centroid before fitting and adding it back afterwards dropped the drift to 0.02 px on the spot. Second: when cleaning edge points with “mean ±3σ rejection”, the σ was computed with the bad points included — the bad points first inflated σ, then sailed through the loosened gate en masse; only after switching to a robust σ estimated by MAD, or to iterative re-estimation after rejection, did the gate actually hold. Third: a line-width measurement station’s repeatability stubbornly hovered around 0.1 px off target; lighting and vibration were investigated to no avail, until it emerged that the profile sampling used nearest-neighbor interpolation — switching to bilinear brought the repeatability within spec immediately. Not one of the three involved any deep algorithm; all of it is in the few pages of this chapter.
2.5 Summary
- Homogeneous coordinates unify “matrix multiplication plus translation” into a single \(3\times 3\) multiplication, and a chained sequence of transformations can be pre-multiplied into one matrix; rigid (3 degrees of freedom), similarity (4), and affine (6) form the 2D transformation family from strict to permissive, preserving lengths, shape, and parallelism respectively.
- The normal equations come from “set the partial derivatives of the objective to zero”: the general form is \(A^\top A\,\boldsymbol{\theta} = A^\top \mathbf{y}\), and the weighted version just multiplies each sum by \(w_i\); the closed-form solution for fitting a line by perpendicular distance is “through the centroid, direction along the principal eigenvector of the scatter matrix”.
- Numerical stability is part of measurement accuracy: \(A^\top A\) squares the condition number, so with large coordinates always center the data before fitting; the squared penalty amplifies the voice of outliers, and robust fitting (Huber, RANSAC) takes over in Chapter 14.
- The Gaussian distribution plus the 3σ rule is the vocabulary of threshold setting, and the central limit theorem explains why sensor noise is nearly Gaussian; the histogram is the empirical estimate of the distribution, and averaging makes noise converge as \(\sigma/\sqrt{N}\) — smoothing, projection, and multi-frame averaging all share this one law.
- Bilinear interpolation weights the four neighbors by diagonally opposite areas; three-point parabolic interpolation gives the vertex offset \(\delta = (g_{-1}-g_{+1})/(2(g_{-1}-2g_0+g_{+1}))\), with \(|\delta|\le 1/2\) always — it is the source of every subpixel formula in this book.
A systematic treatment of the probability, statistics, and Gaussian-model thread is given in (Bishop 2006); the normal equations, condition number, and numerical stability of least squares are classic numerical-linear-algebra material, for which the authoritative reference is (Golub and Van Loan 2013); a quick reference on geometric transformations, linear algebra, and optimization in a vision context can be found in the appendices of (Szeliski 2022). This chapter’s material is developed more completely, and closer to industrial measurement, in the book by Steger et al. (Steger, Ulrich, and Wiedemann 2018), particularly the parts on geometric transformations and on least-squares fitting with its uncertainty analysis, which make good further reading.