Ground Plane Removal in Point Cloud Data

Ground points are usually the majority of an outdoor scan and almost never the thing you are looking for. Removing them typically deletes half to three-quarters of the data and, more importantly, disconnects everything else into separate objects — which is why it is the first step in nearly every pipeline.

Why this is the first step in almost everything

Two effects, and the second is the one that matters. The obvious effect is volume: fewer points means every subsequent stage is faster in direct proportion.

The structural effect is connectivity. With the ground present, every object in the scene is joined to every other object through it, so Euclidean clustering returns one enormous component containing the entire scan. Remove the ground and a car, a tree and a pedestrian become three connected components that a clustering algorithm separates with no training data, no model, and no parameters beyond a distance threshold. That is a large amount of segmentation obtained for the cost of one plane fit.

How RANSAC fits a plane

Least squares cannot fit the ground, because it minimises error against every point and the building, the trees and the cars all pull on it. The fit ends up somewhere between the ground and the roofline, matching nothing.

RANSAC, published by Martin Fischler and Robert Bolles in 1981, inverts the approach. Instead of fitting all the data, it hypothesises a model from the minimum number of points that determine one, then asks how many other points agree. A plane needs three points, so:

1. Sample three points at random.
2. Compute the plane through them: the normal is the cross product of two edge vectors, normalised, and the offset follows from any of the three points.
3. Count how many of all the points lie within a distance threshold of that plane. These are the inliers.
4. Keep the hypothesis with the most inliers.
5. Optionally refit by least squares to the winning inlier set, which improves the estimate without reintroducing the outliers.

The distance threshold is the parameter that matters and it should be set from physics, not by trial: roughly two to three times the scanner’s range noise plus the actual roughness of the surface. Asphalt and a ploughed field need very different values.

The property that makes RANSAC the right tool here is that its cost does not grow with the number of outliers, only with their proportion. A scan containing a million points of building above two million points of ground is no harder than a scan of a hundred points in the same ratio, because the sampling probability depends on the fraction alone. That is exactly the wrong way round for least squares, where every additional non-ground point makes the fit worse. It is also why the crop in the script below matters more than it looks: raising the ground fraction from 30% to 70% cuts the iteration requirement by a factor of fifteen.

Worked: six points by hand

the point set
  p1 = (0.0, 0.0, 0.02)     near-ground
  p2 = (1.0, 0.0, -0.01)    near-ground
  p3 = (0.0, 1.0,  0.00)    near-ground
  p4 = (1.0, 1.0,  0.03)    near-ground
  p5 = (0.5, 0.5,  1.80)    a person's head
  p6 = (0.4, 0.6,  1.20)    the same person's torso

sample p1, p2, p3 and build the plane

  v1 = p2 - p1 = ( 1.0,  0.0, -0.03)
  v2 = p3 - p1 = ( 0.0,  1.0, -0.02)

  n  = v1 x v2
     = (0*(-0.02) - (-0.03)*1,
        (-0.03)*0 - 1*(-0.02),
        1*1 - 0*0)
     = (0.03, 0.02, 1.00)

  |n| = sqrt(0.0009 + 0.0004 + 1.0) = 1.00065
  n_hat = (0.02998, 0.01999, 0.99935)

  d = -(n_hat . p1) = -(0.99935 x 0.02) = -0.019987

  plane: 0.02998x + 0.01999y + 0.99935z - 0.019987 = 0

distance of each point from the plane

  p4: 0.02998(1) + 0.01999(1) + 0.99935(0.03) - 0.019987
    = 0.02998 + 0.01999 + 0.029981 - 0.019987
    = 0.0600 m

  p5: 0.02998(0.5) + 0.01999(0.5) + 0.99935(1.80) - 0.019987
    = 0.014990 + 0.009995 + 1.798830 - 0.019987
    = 1.8038 m

  p6: 0.02998(0.4) + 0.01999(0.6) + 0.99935(1.20) - 0.019987
    = 0.011992 + 0.011994 + 1.199220 - 0.019987
    = 1.2032 m

with a threshold of 0.10 m
  inliers: p1, p2, p3 (exactly on it) and p4 at 0.060 m
  outliers: p5 at 1.80 m, p6 at 1.20 m
  -> 4 of 6 inliers, and the two person points survive

with a threshold of 0.05 m
  p4 at 0.060 m is now excluded: a real ground point
  rejected because the threshold is tighter than the
  surface roughness in this sample

That last case is the failure people meet first. A threshold below the combined range noise and surface roughness rejects genuine ground, leaving a thin carpet of leftover points that the next clustering stage turns into thousands of spurious objects.

How many iterations you actually need

RANSAC’s iteration count is not a guess. To draw at least one all-inlier sample with probability p, given an inlier fraction w and a sample size s:

N = log(1 - p) / log(1 - w^s)     with s = 3 for a plane

p = 0.999
  w = 0.70 (70 % ground):
    w^3 = 0.343
    N = log(0.001) / log(0.657) = -6.908 / -0.420 = 17

  w = 0.50:
    w^3 = 0.125
    N = -6.908 / log(0.875) = -6.908 / -0.1335 = 52

  w = 0.30:
    w^3 = 0.027
    N = -6.908 / log(0.973) = -6.908 / -0.02737 = 253

  w = 0.10:
    w^3 = 0.001
    N = -6.908 / log(0.999) = -6.908 / -0.001 = 6,905

On a scan where the ground is most of the data, seventeen iterations is enough. A default of 1,000 is therefore generous rather than careful, and raising it further does nothing — if a fit is failing at 1,000 iterations on a ground-dominated scan, the model is wrong, not the sampling. The one case where a large count genuinely helps is a scan cropped so tightly that ground is a small minority of the points, which is why cropping before fitting can quietly break a pipeline that worked on the full cloud.

The script

This runs against Open3D and NumPy and writes two files. The guards after the fit are the part that distinguishes it from a snippet: a plane fit with no verification will happily hand you a wall.

import numpy as np
import open3d as o3d

pcd = o3d.io.read_point_cloud("scan.ply")
print("input points:", len(pcd.points))

# 1. downsample. the fit does not need full density, and
#    5 cm voxels make the RANSAC loop far cheaper.
work = pcd.voxel_down_sample(voxel_size=0.05)

# 2. crop to a height band around the expected ground.
#    this raises the inlier fraction, which is the term
#    the iteration count is most sensitive to.
z = np.asarray(work.points)[:, 2]
band = np.where((z > np.percentile(z, 1) - 0.5) &
                (z < np.percentile(z, 1) + 2.0))[0]
work = work.select_by_index(band)

# 3. fit. threshold = 2-3x range noise + surface roughness.
plane_model, inliers = work.segment_plane(
    distance_threshold=0.10,
    ransac_n=3,
    num_iterations=1000,
)
a, b, c, d = plane_model
print("plane: %.4f x + %.4f y + %.4f z + %.4f = 0" % (a, b, c, d))

# 4. GUARD ONE: is the plane roughly horizontal?
tilt_deg = np.degrees(np.arccos(min(1.0, abs(c))))
print("tilt from horizontal: %.2f deg" % tilt_deg)
if tilt_deg > 15.0:
    raise SystemExit("fitted plane is not horizontal - "
                     "probably a wall. check the crop.")

# 5. GUARD TWO: did it claim a plausible share of points?
frac = len(inliers) / len(work.points)
print("inlier fraction: %.3f" % frac)
if frac < 0.20:
    raise SystemExit("ground is under 20 % of the cropped "
                     "cloud - the model is probably wrong.")

# 6. apply the plane to the FULL cloud, not the downsampled
#    one, so no resolution is lost in the output.
pts = np.asarray(pcd.points)
dist = pts @ np.array([a, b, c]) + d
ground_mask = np.abs(dist) <= 0.10

ground = pcd.select_by_index(np.where(ground_mask)[0])
rest   = pcd.select_by_index(np.where(~ground_mask)[0])

print("ground: %d  non-ground: %d" % (len(ground.points),
                                      len(rest.points)))
o3d.io.write_point_cloud("ground.ply", ground)
o3d.io.write_point_cloud("non_ground.ply", rest)

Step 6 is the detail most snippets get wrong. Fitting on a downsampled cloud is correct and fast, but selecting the output from that same downsampled cloud silently throws away everything the voxel filter removed. Fit on the small cloud, apply to the large one.

When one plane is the wrong model

A single plane assumes the ground is flat within the distance threshold across the whole scan. Over 20 m of a car park that is true. Over 500 m of a road with a 4% gradient it is not: the ground rises 20 m end to end, and no single plane with a 10 cm threshold can contain it.

• Tile it. Split into 20 m × 20 m tiles and fit a plane per tile. Simple, parallel, and effective on gently varying terrain. Fit on a slightly larger window than the tile you keep, so the planes agree across the seams.
• Radial bins, for a vehicle-mounted sensor: divide by azimuth and range into wedge-shaped segments and fit a line per segment. This matches the sensor’s own geometry and handles the varying density that comes with it.
• Progressive morphological filtering, described by Keqi Zhang and colleagues in IEEE Transactions on Geoscience and Remote Sensing in 2003, opens the surface with a window that grows over successive passes, removing objects by size while preserving terrain relief. It is the classical choice in airborne survey processing and it handles hills that no plane model can.
• Cloth simulation filtering inverts the cloud and drops a simulated cloth onto it; where the cloth settles is the terrain. It has few parameters, handles undulation naturally, and is widely used in forestry and topographic work.

Whichever model you use, verify the output before trusting it. The cheapest check is that the non-ground cloud, clustered by Euclidean distance, produces a plausible number of components of plausible size — thousands of tiny fragments means the threshold was too tight, and one giant component means the ground is still connected. That check belongs in the pipeline alongside the ones on the quality gate page.