RANSAC Applications

The previous section derived RANSAC on a line: the number of trials

L

for a minimal sample of size

s

and inlier fraction

\pi

, the threshold

\tau

as a chi-squared gate, and the refit on the consensus set. This section applies the same steps to problems from robotics and vision. It finds the floor plane in a depth camera’s point cloud and checks the settings of the homography estimate in an earlier section. It ends with the way RANSAC and the ICP algorithm of scan registration work together to align two point clouds. The code below repeats the RANSAC loop and the formula for the number of trials from the previous section, so that this section runs on its own.

import numpy as np
from scipy.stats import chi2, norm

def ransac(m, s, fit_minimal, residuals, tau, L, rng, score="inliers", L_adaptive=None):
    # fit_minimal(idx) returns theta or None for a degenerate sample; residuals(theta) returns e for all points.
    # score: "inliers" (count), "msac" (truncated squared error) or a callable cost(e) to minimize.
    # L_adaptive(n_inliers) optionally returns an updated number of trials after each new best.
    best_theta, best_cost, best_inliers, trial, history = None, np.inf, None, 0, []
    while trial < L:
        trial += 1
        theta = fit_minimal(rng.choice(m, s, replace=False))
        if theta is None:
            continue
        e = residuals(theta)
        inliers = np.abs(e) < tau
        if score == "inliers":
            cost = -inliers.sum()
        elif score == "msac":
            cost = np.minimum(e**2, tau**2).sum()
        else:
            cost = score(e)
        if cost < best_cost:
            best_theta, best_cost, best_inliers = theta, cost, inliers
            if L_adaptive is not None:
                L = min(L, L_adaptive(inliers.sum()))
            history.append((trial, inliers.sum(), L))
    return best_theta, best_inliers, trial, history



def trials_needed(pi, s, P):
    # smallest L with (1 - pi^s)^L <= 1 - P
    return int(np.ceil(np.log(1 - P) / np.log1p(-(pi**s))))

Plane fitting in a point cloud

A depth camera on a mobile robot sees the floor, and finding it is often the first step of obstacle detection. The floor is a plane

\mathbf n^\top \mathbf p = d

with unit normal

\mathbf n

, and three points that are not collinear determine it, so

s = 3

. The residual of a point is its distance to the plane, which is one-dimensional, so the threshold is

1.96\,\sigma

at 95%, as it was for the line. The synthetic cloud below has 300 floor points with 1 cm of noise on a slightly tilted plane, 150 points on the faces of two boxes standing on the floor, and 50 points scattered at random.

def make_floor_cloud(seed=3, sigma_floor=0.01):
    rng = np.random.default_rng(seed)
    n_true = np.array([0.02, -0.01, 1.0]) / np.linalg.norm([0.02, -0.01, 1.0])
    xy = rng.uniform(0, 4, (300, 2))
    z = (0.0 - n_true[0] * xy[:, 0] - n_true[1] * xy[:, 1]) / n_true[2] + sigma_floor * rng.standard_normal(300)
    floor = np.column_stack([xy, z])
    boxes = []
    for (cx, cy, w, h) in [(1.0, 1.0, 0.6, 0.5), (2.8, 2.5, 0.8, 0.3)]:
        face = rng.integers(0, 3, 75)
        u, v = rng.uniform(-0.5, 0.5, (2, 75))
        p = np.column_stack([cx + w * u, cy + w * v, h * np.ones(75)])  # top face
        side = face > 0
        p[side, 2] = rng.uniform(0, h, side.sum())
        p[side & (face == 1), 0] = cx + w / 2
        p[side & (face == 2), 1] = cy + w / 2
        boxes.append(p)
    clutter = rng.uniform([0, 0, 0.05], [4, 4, 1.0], (50, 3))
    cloud = np.vstack([floor] + boxes + [clutter])
    is_floor = np.r_[np.ones(300, bool), np.zeros(len(cloud) - 300, bool)]
    return cloud, is_floor, n_true


def plane_through(idx, pts):
    # plane through three points as (n, d) with |n| = 1; None if the points are nearly collinear
    p1, p2, p3 = pts[idx]
    n = np.cross(p2 - p1, p3 - p1)
    norm_n = np.linalg.norm(n)
    if norm_n < 1e-6:
        return None
    n = n / norm_n
    return np.r_[n, n @ p1]


def plane_least_squares(pts):
    # total least squares plane: normal = direction of least variance
    c = pts.mean(axis=0)
    n = np.linalg.svd(pts - c)[2][-1]
    return np.r_[n, n @ c]


cloud, is_floor, n_true = make_floor_cloud()
sigma_floor = 0.01
tau_plane = sigma_floor * np.sqrt(chi2.ppf(0.95, 1))
plane_res = lambda th: cloud @ th[:3] - th[3]
pi_floor = is_floor.mean()
L_plane = trials_needed(pi_floor, 3, 0.99)
plane_r, inl_p, _, _ = ransac(len(cloud), 3, lambda idx: plane_through(idx, cloud), plane_res, tau_plane, L=L_plane,
                              rng=np.random.default_rng(0))
plane_r = plane_least_squares(cloud[inl_p])
inl_p = np.abs(plane_res(plane_r)) < tau_plane  # consensus set of the refit plane
angle = np.degrees(np.arccos(abs(plane_r[:3] @ n_true)))
print(f"{len(cloud)} points, floor fraction {pi_floor:.2f}, L = {L_plane} trials for P = 0.99, tau = {100 * tau_plane:.2f} cm")
print(f"consensus set: {inl_p.sum()} points, {(inl_p & is_floor).sum()} floor and {(inl_p & ~is_floor).sum()} other")
print(f"normal error after refit: {angle:.3f} degrees")
print(f"least squares plane through all points, normal error: "
      f"{np.degrees(np.arccos(abs(plane_least_squares(cloud)[:3] @ n_true))):.1f} degrees")

500 points, floor fraction 0.60, L = 19 trials for P = 0.99, tau = 1.96 cm
consensus set: 280 points, 280 floor and 0 other
normal error after refit: 0.217 degrees
least squares plane through all points, normal error: 1.2 degrees

Three-dimensional scatter of a point cloud. Blue points form a flat floor; red points form two box shapes standing on it and a scatter of clutter above it.

RANSAC plane fit with $s = 3$ . Blue points are the consensus set of the floor plane.

Homographies: four correspondences

The homographies section estimates a planar homography from point correspondences with RANSAC. A homography has

M = 8

parameters, and each correspondence gives two equations, so a minimal sample is

s = 4

correspondences. That section runs 400 trials with a threshold of 0.20 on the reprojection error, and reports that RANSAC accepted 25 of 36 correspondences, an inlier ratio of 0.694. Its synthetic destination points have Gaussian noise with standard deviation 0.08 in each coordinate. The results of the previous section explain both settings.

pi_h = 25 / 36
print(f"trials for pi = {pi_h:.3f}, s = 4: L = {trials_needed(pi_h, 4, 0.99)} at P = 0.99, {trials_needed(pi_h, 4, 0.999)} at P = 0.999")
print(f"95% chi-squared gate for a 2D error with sigma = 0.08: tau = {0.08 * np.sqrt(chi2.ppf(0.95, 2)):.3f}")
print(f"trials at the parameter study's largest outlier ratio, pi = 0.3: L = {trials_needed(0.3, 4, 0.99)} at P = 0.99")

trials for pi = 0.694, s = 4: L = 18 at P = 0.99, 27 at P = 0.999
95% chi-squared gate for a 2D error with sigma = 0.08: tau = 0.196
trials at the parameter study's largest outlier ratio, pi = 0.3: L = 567 at P = 0.99

The threshold of 0.20 is the 95% gate of a two-dimensional error at that noise level, and 18 trials would have reached 99% confidence where the section runs 400. The section’s parameter study uses 500 trials up to an outlier ratio of 0.7, where

\pi = 0.3

requires 567 for

P = 0.99

, so its last point is slightly under-sampled. The same reasoning applies to the essential matrix in the visual SLAM section, where RANSAC keeps 97 of 113 matches between the first two keyframes with a five-point minimal sample.

print(f"five-point essential matrix, pi = 97/113 = {97 / 113:.2f}: L = {trials_needed(97 / 113, 5, 0.999)} at P = 0.999")

five-point essential matrix, pi = 97/113 = 0.86: L = 12 at P = 0.999

RANSAC and ICP

The scan registration section aligns two point clouds with the iterative closest point algorithm (ICP) (Besl & McKay, 1992). ICP alternates two steps: pair every point of one scan with its nearest neighbor in the other, then solve for the rigid motion that best aligns the pairs. This is the same loop as k-means. The nearest-neighbor pairing is a hard assignment like

r_{ik}

, and the rigid fit plays the role of recomputing the cluster means. ICP is hard EM over correspondences, and variants that weight every candidate pair by a probability are EM proper. Like k-means and EM, it climbs to the nearest local optimum and needs a good starting alignment. On matched features, such as descriptor matches between two scans or two images, RANSAC draws minimal samples of correspondences (three for a rigid motion in 3D). It needs no starting alignment and tolerates many wrong matches, but its accuracy is limited by the matches it keeps. Registration pipelines therefore use both: RANSAC on feature matches gives a coarse global alignment, and ICP refines it on all the points. The correspondence estimation section uses the same split for SLAM: descriptor matches are checked against the geometry, and only the matches that pass go to the filter. Key references: (Hartley & Zisserman, 2004; Besl & McKay, 1992)

References

Besl, P., McKay, N. (1992). A method for registration of 3-D shapes. IEEE Trans. Pattern Anal. Mach. Intell.
Hartley, R., Zisserman, A. (2004). Multiple View Geometry in Computer Vision.

Found an error on this page? Report an issue.

Foundations

Learning & Regression

Optimization

Maximum Likelihood

Classification

Dimensionality Reduction

Mixture Models and EM

Robust Estimation

RANSAC Applications

Plane fitting in a point cloud

Homographies: four correspondences

RANSAC and ICP

References

​Plane fitting in a point cloud

​Homographies: four correspondences

​RANSAC and ICP

​References

Plane fitting in a point cloud

Homographies: four correspondences

RANSAC and ICP

References