{ "cells": [ { "cell_type": "markdown", "id": "add51d1a", "metadata": {}, "source": [ "# Lesson 22: Model Fitting, Robust Estimation, and RANSAC\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": null, "id": "d33e4310", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import cv2\n", "import matplotlib.pyplot as plt" ] }, { "cell_type": "markdown", "id": "54b4f929", "metadata": {}, "source": [ "## Ordinary vs. total least squares\n", "\n", "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.\n", "\n", "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.)" ] }, { "cell_type": "code", "execution_count": null, "id": "e424ca75", "metadata": {}, "outputs": [], "source": [ "rng_tls = np.random.default_rng(2)\n", "true_m, true_b = 3.0, 1.0 # a steep line, where OLS struggles most\n", "t = rng_tls.uniform(0, 5, 40)\n", "x_noisy = t + rng_tls.normal(0, 0.4, 40)\n", "y_noisy = true_m * t + true_b + rng_tls.normal(0, 0.4, 40)\n", "\n", "A_line = np.vstack([x_noisy, np.ones_like(x_noisy)]).T\n", "slope_ols, intercept_ols = np.linalg.lstsq(A_line, y_noisy, rcond=None)[0]\n", "\n", "points = np.stack([x_noisy, y_noisy], axis=1)\n", "centroid = points.mean(axis=0)\n", "cov = np.cov((points - centroid).T)\n", "eigvals, eigvecs = np.linalg.eigh(cov) # ascending order\n", "direction = eigvecs[:, -1] # largest eigenvalue = direction of max spread\n", "slope_tls = direction[1] / direction[0]\n", "intercept_tls = centroid[1] - slope_tls * centroid[0]\n", "\n", "print(f'true line: y = {true_m:.2f}x + {true_b:.2f}')\n", "print(f'OLS fit: y = {slope_ols:.2f}x + {intercept_ols:.2f} <- biased, since x is noisy too')\n", "print(f'TLS fit: y = {slope_tls:.2f}x + {intercept_tls:.2f} <- closer to the true slope')\n", "\n", "xs_tls = np.array([x_noisy.min(), x_noisy.max()])\n", "plt.scatter(x_noisy, y_noisy, s=15, alpha=0.6, label='noisy points (both x and y)')\n", "plt.plot(xs_tls, true_m * xs_tls + true_b, '--', color='gray', label='true line')\n", "plt.plot(xs_tls, slope_ols * xs_tls + intercept_ols, color='tab:orange', label='OLS (vertical error)')\n", "plt.plot(xs_tls, slope_tls * xs_tls + intercept_tls, color='tab:green', label='TLS (perpendicular error)')\n", "plt.legend(fontsize=8)\n", "plt.title('OLS vs. TLS when both x and y are noisy')\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "94803ef4", "metadata": {}, "source": [ "TLS noticeably recovers the true slope better here, since the noise is genuinely symmetric in $x$ and $y$." ] }, { "cell_type": "markdown", "id": "3099b493", "metadata": {}, "source": [ "## Why least squares breaks under outliers\n", "\n", "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." ] }, { "cell_type": "code", "execution_count": null, "id": "b2e60e5f", "metadata": {}, "outputs": [], "source": [ "rng = np.random.default_rng(0)\n", "true_slope, true_intercept = 2.0, 5.0\n", "\n", "x_inliers = rng.uniform(0, 10, 40)\n", "y_inliers = true_slope * x_inliers + true_intercept + rng.normal(0, 0.5, 40)\n", "\n", "x_outliers = rng.uniform(0, 10, 15)\n", "y_outliers = rng.uniform(-20, 40, 15) # unrelated to the line at all\n", "\n", "x = np.concatenate([x_inliers, x_outliers])\n", "y = np.concatenate([y_inliers, y_outliers])\n", "\n", "A = np.vstack([x, np.ones_like(x)]).T\n", "slope_ols, intercept_ols = np.linalg.lstsq(A, y, rcond=None)[0]\n", "\n", "print(f'true line: y = {true_slope:.2f}x + {true_intercept:.2f}')\n", "print(f'ordinary least squares: y = {slope_ols:.2f}x + {intercept_ols:.2f} <- dragged off by outliers')\n", "\n", "xs = np.array([0, 10])\n", "plt.scatter(x_inliers, y_inliers, c='tab:blue', label='inliers')\n", "plt.scatter(x_outliers, y_outliers, c='tab:red', marker='x', label='outliers')\n", "plt.plot(xs, true_slope * xs + true_intercept, '--', color='gray', label='true line')\n", "plt.plot(xs, slope_ols * xs + intercept_ols, color='tab:orange', label='OLS fit')\n", "plt.legend(fontsize=8)\n", "plt.title('Ordinary least squares, corrupted by 15 outliers among 40 inliers')\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "c22a3203", "metadata": {}, "source": [ "## RANSAC: fit from the inside out\n", "\n", "**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)." ] }, { "cell_type": "code", "execution_count": null, "id": "0d337cee", "metadata": {}, "outputs": [], "source": [ "def ransac_line(x, y, threshold=2.0, n_iterations=200, rng=None):\n", " rng = rng or np.random.default_rng()\n", " best_inliers, best_count = None, -1\n", "\n", " for _ in range(n_iterations):\n", " i, j = rng.choice(len(x), 2, replace=False) # minimal sample: 2 points define a line\n", " if x[i] == x[j]:\n", " continue\n", " m = (y[j] - y[i]) / (x[j] - x[i])\n", " b = y[i] - m * x[i]\n", " distance = np.abs(m * x - y + b) / np.sqrt(m**2 + 1)\n", " inliers = distance < threshold\n", " if inliers.sum() > best_count:\n", " best_count, best_inliers = inliers.sum(), inliers\n", "\n", " A = np.vstack([x[best_inliers], np.ones(best_inliers.sum())]).T\n", " slope, intercept = np.linalg.lstsq(A, y[best_inliers], rcond=None)[0] # final refit, inliers only\n", " return slope, intercept, best_inliers\n", "\n", "slope_ransac, intercept_ransac, inlier_mask = ransac_line(x, y, rng=np.random.default_rng(1))\n", "\n", "print(f'true line: y = {true_slope:.2f}x + {true_intercept:.2f}')\n", "print(f'RANSAC fit: y = {slope_ransac:.2f}x + {intercept_ransac:.2f}')\n", "print(f'inliers found: {inlier_mask.sum()} / {len(x)} (planted {40} true inliers)')\n", "\n", "plt.scatter(x[inlier_mask], y[inlier_mask], c='tab:blue', label='found inliers')\n", "plt.scatter(x[~inlier_mask], y[~inlier_mask], c='tab:red', marker='x', label='rejected as outliers')\n", "plt.plot(xs, slope_ransac * xs + intercept_ransac, color='tab:green', label='RANSAC fit')\n", "plt.legend(fontsize=8)\n", "plt.title('RANSAC recovers the true line despite 27% outlier contamination')\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "684587d2", "metadata": {}, "source": [ "## How many iterations does RANSAC need?\n", "\n", "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", "\n", "$$N = \\frac{\\log(1-p)}{\\log(1-w^n)}$$\n", "\n", "More outliers, or a larger minimal sample size, cause this required number of samples to increase." ] }, { "cell_type": "code", "execution_count": null, "id": "18cbaca6", "metadata": {}, "outputs": [], "source": [ "def required_iterations(inlier_fraction, sample_size, confidence=0.99):\n", " return np.log(1 - confidence) / np.log(1 - inlier_fraction**sample_size)\n", "\n", "print(f'{\"inlier %\":>10} {\"line (n=2)\":>12} {\"homography (n=4)\":>18} {\"fundamental matrix (n=8)\":>18}')\n", "for w in [0.9, 0.7, 0.5, 0.3]:\n", " n_line = int(np.ceil(required_iterations(w, 2)))\n", " n_homog = int(np.ceil(required_iterations(w, 4)))\n", " n_fund = int(np.ceil(required_iterations(w, 8)))\n", " print(f'{100*w:>9.0f}% {n_line:>12} {n_homog:>18} {n_fund:>18}')" ] }, { "cell_type": "markdown", "id": "fe15c5cc", "metadata": {}, "source": [ "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." ] }, { "cell_type": "markdown", "id": "52afd09b", "metadata": {}, "source": [ "## A gentler alternative: robust loss functions\n", "\n", "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." ] }, { "cell_type": "code", "id": "c5fe888a", "source": "def huber_reweighted_fit(x, y, delta=2.0, n_iterations=10):\n weights = np.ones_like(x)\n for _ in range(n_iterations):\n sw = np.sqrt(weights)\n A = np.vstack([x, np.ones_like(x)]).T * sw[:, None] # weighted least squares\n slope, intercept = np.linalg.lstsq(A, y * sw, rcond=None)[0]\n residuals = np.abs(y - (slope * x + intercept))\n weights = np.where(residuals > delta, delta / np.maximum(residuals, 1e-6), 1.0)\n return slope, intercept\n\nslope_huber, intercept_huber = huber_reweighted_fit(x, y)\n\nprint(f'true line: y = {true_slope:.2f}x + {true_intercept:.2f}')\nprint(f'OLS fit: y = {slope_ols:.2f}x + {intercept_ols:.2f}')\nprint(f'Huber fit: y = {slope_huber:.2f}x + {intercept_huber:.2f}')\nprint(f'RANSAC fit: y = {slope_ransac:.2f}x + {intercept_ransac:.2f}')\n\nplt.scatter(x_inliers, y_inliers, c='tab:blue', label='inliers')\nplt.scatter(x_outliers, y_outliers, c='tab:red', marker='x', label='outliers')\nplt.plot(xs, true_slope * xs + true_intercept, '--', color='gray', label='true line')\nplt.plot(xs, slope_ols * xs + intercept_ols, color='tab:orange', label='OLS')\nplt.plot(xs, slope_huber * xs + intercept_huber, color='tab:purple', label='Huber-reweighted')\nplt.plot(xs, slope_ransac * xs + intercept_ransac, color='tab:green', label='RANSAC')\nplt.legend(fontsize=8)\nplt.title('Same contaminated data, three different fits')\nplt.show()", "metadata": {}, "execution_count": null, "outputs": [] }, { "cell_type": "markdown", "id": "e89edd12", "source": "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.", "metadata": {} }, { "cell_type": "markdown", "id": "33384aac", "metadata": {}, "source": "### Exercise\n\n1. 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 the `required_iterations` formula to check whether 200 is even theoretically sufficient at that contamination level.\n2. Lower `ransac_line`'s distance `threshold` from `2.0` to `0.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)?\n3. Try `huber_reweighted_fit` with a much smaller `delta` (e.g. `0.5`) and a much larger one (e.g. `5.0`). How does `delta` trade off between OLS-like behavior (barely any downweighting) and RANSAC-like hard rejection (almost binary in/out)?" } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "name": "python", "version": "3.x" } }, "nbformat": 4, "nbformat_minor": 5 }