Lesson 22: Model Fitting, Robust Estimation, and RANSAC¶
Many computer vision problems boil down to fitting a model — a line, a circle, a plane, a homography (Lesson 23), a fundamental matrix (Lesson 26) — to a set of corresponding points. In practice, that data almost always contains outliers: correspondences that are simply wrong. This lesson confronts that directly: least-squares fitting is catastrophically fragile to even a few bad points, and RANSAC is the standard fix.
import numpy as np
import cv2
import matplotlib.pyplot as plt
Ordinary vs. total least squares¶
Fitting $y=mx+b$ by minimizing $\sum(y_i-(mx_i+b))^2$ — ordinary least squares (OLS) — measures error vertically. That's the right choice when $x$ is known exactly and only $y$ is noisy, but not when both coordinates carry comparable noise, as they typically do for 2D image points: the fit gets biased toward whichever direction has less spread. (And it can't even represent a vertical line.) Total least squares (TLS) instead minimizes the perpendicular distance from each point to the line, treating both coordinates symmetrically.
For a 2D line, TLS needs no new machinery: the best-fit line passes through the centroid, and its direction is the eigenvector of the (centered) points' covariance matrix with the largest eigenvalue — exactly Lesson 6's covariance-and-eigenvectors idea, just picking out the axis of maximum spread instead of describing a blob's shape. (With more unknowns — such as with a homography or fundamental matrix — the approach generalizes to the singular value decomposition (SVD), introduced in Lesson 26.)
rng_tls = np.random.default_rng(2)
true_m, true_b = 3.0, 1.0 # a steep line, where OLS struggles most
t = rng_tls.uniform(0, 5, 40)
x_noisy = t + rng_tls.normal(0, 0.4, 40)
y_noisy = true_m * t + true_b + rng_tls.normal(0, 0.4, 40)
A_line = np.vstack([x_noisy, np.ones_like(x_noisy)]).T
slope_ols, intercept_ols = np.linalg.lstsq(A_line, y_noisy, rcond=None)[0]
points = np.stack([x_noisy, y_noisy], axis=1)
centroid = points.mean(axis=0)
cov = np.cov((points - centroid).T)
eigvals, eigvecs = np.linalg.eigh(cov) # ascending order
direction = eigvecs[:, -1] # largest eigenvalue = direction of max spread
slope_tls = direction[1] / direction[0]
intercept_tls = centroid[1] - slope_tls * centroid[0]
print(f'true line: y = {true_m:.2f}x + {true_b:.2f}')
print(f'OLS fit: y = {slope_ols:.2f}x + {intercept_ols:.2f} <- biased, since x is noisy too')
print(f'TLS fit: y = {slope_tls:.2f}x + {intercept_tls:.2f} <- closer to the true slope')
xs_tls = np.array([x_noisy.min(), x_noisy.max()])
plt.scatter(x_noisy, y_noisy, s=15, alpha=0.6, label='noisy points (both x and y)')
plt.plot(xs_tls, true_m * xs_tls + true_b, '--', color='gray', label='true line')
plt.plot(xs_tls, slope_ols * xs_tls + intercept_ols, color='tab:orange', label='OLS (vertical error)')
plt.plot(xs_tls, slope_tls * xs_tls + intercept_tls, color='tab:green', label='TLS (perpendicular error)')
plt.legend(fontsize=8)
plt.title('OLS vs. TLS when both x and y are noisy')
plt.show()
TLS noticeably recovers the true slope better here, since the noise is genuinely symmetric in $x$ and $y$.
Why least squares breaks under outliers¶
Least-squares fitting minimizes the sum of squared residuals. But squaring has a drawback: every outlier (point far from the model) contributes an enormous amount to the total error, causing the whole fit to move away from all the other, correct points.
rng = np.random.default_rng(0)
true_slope, true_intercept = 2.0, 5.0
x_inliers = rng.uniform(0, 10, 40)
y_inliers = true_slope * x_inliers + true_intercept + rng.normal(0, 0.5, 40)
x_outliers = rng.uniform(0, 10, 15)
y_outliers = rng.uniform(-20, 40, 15) # unrelated to the line at all
x = np.concatenate([x_inliers, x_outliers])
y = np.concatenate([y_inliers, y_outliers])
A = np.vstack([x, np.ones_like(x)]).T
slope_ols, intercept_ols = np.linalg.lstsq(A, y, rcond=None)[0]
print(f'true line: y = {true_slope:.2f}x + {true_intercept:.2f}')
print(f'ordinary least squares: y = {slope_ols:.2f}x + {intercept_ols:.2f} <- dragged off by outliers')
xs = np.array([0, 10])
plt.scatter(x_inliers, y_inliers, c='tab:blue', label='inliers')
plt.scatter(x_outliers, y_outliers, c='tab:red', marker='x', label='outliers')
plt.plot(xs, true_slope * xs + true_intercept, '--', color='gray', label='true line')
plt.plot(xs, slope_ols * xs + intercept_ols, color='tab:orange', label='OLS fit')
plt.legend(fontsize=8)
plt.title('Ordinary least squares, corrupted by 15 outliers among 40 inliers')
plt.show()
RANSAC: fit from the inside out¶
RANSAC (RANdom SAmple Consensus, Fischler & Bolles, 1981★) flips the strategy: instead of using all the data and hoping outliers don't matter, it repeatedly picks the smallest possible random subset needed to define a candidate model, counts how many of the remaining points agree with it (the consensus set), and keeps whichever candidate has the most agreement. A final least-squares refit on the winning inlier set gives the polished result (since least squares works well without outliers).
def ransac_line(x, y, threshold=2.0, n_iterations=200, rng=None):
rng = rng or np.random.default_rng()
best_inliers, best_count = None, -1
for _ in range(n_iterations):
i, j = rng.choice(len(x), 2, replace=False) # minimal sample: 2 points define a line
if x[i] == x[j]:
continue
m = (y[j] - y[i]) / (x[j] - x[i])
b = y[i] - m * x[i]
distance = np.abs(m * x - y + b) / np.sqrt(m**2 + 1)
inliers = distance < threshold
if inliers.sum() > best_count:
best_count, best_inliers = inliers.sum(), inliers
A = np.vstack([x[best_inliers], np.ones(best_inliers.sum())]).T
slope, intercept = np.linalg.lstsq(A, y[best_inliers], rcond=None)[0] # final refit, inliers only
return slope, intercept, best_inliers
slope_ransac, intercept_ransac, inlier_mask = ransac_line(x, y, rng=np.random.default_rng(1))
print(f'true line: y = {true_slope:.2f}x + {true_intercept:.2f}')
print(f'RANSAC fit: y = {slope_ransac:.2f}x + {intercept_ransac:.2f}')
print(f'inliers found: {inlier_mask.sum()} / {len(x)} (planted {40} true inliers)')
plt.scatter(x[inlier_mask], y[inlier_mask], c='tab:blue', label='found inliers')
plt.scatter(x[~inlier_mask], y[~inlier_mask], c='tab:red', marker='x', label='rejected as outliers')
plt.plot(xs, slope_ransac * xs + intercept_ransac, color='tab:green', label='RANSAC fit')
plt.legend(fontsize=8)
plt.title('RANSAC recovers the true line despite 27% outlier contamination')
plt.show()
How many iterations does RANSAC need?¶
If a fraction $w$ of the data are inliers, and the model needs a minimal sample of $n$ points, the probability that any one random sample is entirely inliers is $w^n$. To be at least p confident of drawing an all-inlier sample at least once across $N$ independent tries:
$$N = \frac{\log(1-p)}{\log(1-w^n)}$$
More outliers, or a larger minimal sample size, cause this required number of samples to increase.
def required_iterations(inlier_fraction, sample_size, confidence=0.99):
return np.log(1 - confidence) / np.log(1 - inlier_fraction**sample_size)
print(f'{"inlier %":>10} {"line (n=2)":>12} {"homography (n=4)":>18} {"fundamental matrix (n=8)":>18}')
for w in [0.9, 0.7, 0.5, 0.3]:
n_line = int(np.ceil(required_iterations(w, 2)))
n_homog = int(np.ceil(required_iterations(w, 4)))
n_fund = int(np.ceil(required_iterations(w, 8)))
print(f'{100*w:>9.0f}% {n_line:>12} {n_homog:>18} {n_fund:>18}')
Fitting a line only needs 2 points, so even fairly heavy contamination (50% outliers) needs just a few dozen iterations. Fitting a homography needs a minimal sample of 4 points, so the same inlier fraction needs an order of magnitude more iterations — the price of a more complex model. Similarly for fitting a fundamental matrix, which needs a minimal sample of 8 points.
A gentler alternative: robust loss functions¶
RANSAC makes a hard inlier/outlier decision. An alternative family, M-estimators, instead reweights every point's contribution smoothly — e.g. the Huber loss behaves like ordinary squared error for small residuals but switches to linear (much less aggressive) growth beyond a threshold, so a single very-wrong point can no longer dominate the total cost the way it does in Loss = residual^2. RANSAC and M-estimators are complementary in practice: RANSAC is excellent at rejecting gross outliers (completely wrong matches), while an M-estimator refinement afterward can down-weight smaller, more subtle deviations among the remaining inliers.
def huber_reweighted_fit(x, y, delta=2.0, n_iterations=10):
weights = np.ones_like(x)
for _ in range(n_iterations):
sw = np.sqrt(weights)
A = np.vstack([x, np.ones_like(x)]).T * sw[:, None] # weighted least squares
slope, intercept = np.linalg.lstsq(A, y * sw, rcond=None)[0]
residuals = np.abs(y - (slope * x + intercept))
weights = np.where(residuals > delta, delta / np.maximum(residuals, 1e-6), 1.0)
return slope, intercept
slope_huber, intercept_huber = huber_reweighted_fit(x, y)
print(f'true line: y = {true_slope:.2f}x + {true_intercept:.2f}')
print(f'OLS fit: y = {slope_ols:.2f}x + {intercept_ols:.2f}')
print(f'Huber fit: y = {slope_huber:.2f}x + {intercept_huber:.2f}')
print(f'RANSAC fit: y = {slope_ransac:.2f}x + {intercept_ransac:.2f}')
plt.scatter(x_inliers, y_inliers, c='tab:blue', label='inliers')
plt.scatter(x_outliers, y_outliers, c='tab:red', marker='x', label='outliers')
plt.plot(xs, true_slope * xs + true_intercept, '--', color='gray', label='true line')
plt.plot(xs, slope_ols * xs + intercept_ols, color='tab:orange', label='OLS')
plt.plot(xs, slope_huber * xs + intercept_huber, color='tab:purple', label='Huber-reweighted')
plt.plot(xs, slope_ransac * xs + intercept_ransac, color='tab:green', label='RANSAC')
plt.legend(fontsize=8)
plt.title('Same contaminated data, three different fits')
plt.show()
The Huber fit lands between the two: much closer to the true line than OLS, though not quite as exact as RANSAC on data this heavily contaminated, since a handful of gross outliers still pull a little on every iteration rather than being cut out entirely.
Exercise¶
- Increase
ransac_line's outlier count until roughly 70% of the points are outliers. Does 200 iterations remain enough to reliably recover the true line? Use therequired_iterationsformula to check whether 200 is even theoretically sufficient at that contamination level. - Lower
ransac_line's distancethresholdfrom2.0to0.5. Does the number of found inliers change, and why might too-strict a threshold actually start rejecting genuine inliers (hint: think about what noise, not outliers, does to correct points)? - Try
huber_reweighted_fitwith a much smallerdelta(e.g.0.5) and a much larger one (e.g.5.0). How doesdeltatrade off between OLS-like behavior (barely any downweighting) and RANSAC-like hard rejection (almost binary in/out)?