Farthest-point sampling

Random picks follow the crowd; farthest picks follow the shape

Many methods need a smaller set of points from a large cloud: a model that accepts a fixed number of points, or a set of centers to group nearby points around. A random sample keeps the cloud's density. If one part of a scan holds ten times as many points, it gets about ten times as many picks, and a thin part may get none.

Farthest-point sampling (FPS) spreads its picks over the space the shape fills, regardless of how crowded each part is. The quantity it watches is each point's gap: its distance to the nearest point picked so far. A point with a large gap is poorly covered.

Keep one gap per point, and pick the largest each round

1. Pick a starting point. Set every point's gap to its distance from that start.
2. Pick the point with the largest gap.
3. Measure every point's distance to the new pick. Where it is smaller than the point's gap, it becomes the new gap.
4. Repeat steps 2 and 3 until you have as many points as you need.

Step 3 is cheap: the nearest pick of each point either stays the same or becomes the new pick, so one comparison per point updates every gap.

fps.py

import numpy as np

def farthest_point_sampling(pts, k, start=0):
    chosen = [start]
    gap = np.sqrt(((pts - pts[start]) ** 2).sum(axis=1))   # distance to the chosen set
    for _ in range(k - 1):
        # argmax() finds the index of the largest gap, and int() makes it a plain
        # Python number. append adds that index to the end of the list chosen.
        i = int(gap.argmax())                               # the worst-covered point
        chosen.append(i)
        d = np.sqrt(((pts - pts[i]) ** 2).sum(axis=1))
        gap = np.minimum(gap, d)                            # the new point may be closer
    return chosen

pts = np.array([[0, 0], [1, 0], [2, 0], [6, 0], [6, 3]], dtype=float)
print(farthest_point_sampling(pts, 4))

[0, 4, 3, 2]

Indexes 0, 4, 3, 2 are A, E, D, C. pts - pts[i] subtracts one point from every row at once, and broadcasting explains how. np.minimum takes the smaller of two numbers at each position, which is step 3 for every point at once. The name _ is a loop variable that the loop never uses.

Each pick measures one distance per point, so K picks from N points cost K × N distances. Picking 1,024 points from 100,000 takes about 10⁸ distances: work that a GPU does quickly, since every point's update is the same simple step.