An algorithm for determining the minimum distance between two non-convex polygons

minimum_polygon_distance.md

Adapted from https://www.quora.com/How-do-you-compute-the-distance-between-two-non-convex-polygons-in-linear-time/answer/Tom-Dreyfus

Calculating the Minimum Distance Between Two Non-Convex Polygons

See Amato, Nancy M (1992) for details of the difference between separation (sigma: σ) and closest visible vertex (CVV).

Refer to P and Q as the two polygons with n and m vertices, respectively.
For the purposes of this discussion, a key insight is that it is enough to find the closest edge to each vertex in order to compute the minimum separation between P and Q.
This means iterating over all vertices, and finding a nearest neighbour. Thus, a time complexity in O((m + n) * log(m * n)) should be expected.

It should be noted that the minimum distance involves at least one vertex:

• If the nearest edges of each polygon are parallel, you can slide along the edges so that one of the points is a vertex of one of the polygons
• if not, it is obvious that one of the points is a vertex of one of the polygons.

Step 1

In order to locate vertices of Q within P, we compute a triangulation of P such that edges of P are edges of the triangulation. This is known as a constrained triangulation:

• First, a Delaunay triangulation is built from the vertices of P
• Edges of P are inserted one by one by removing triangles from the triangulation
• Building this triangulation has an output sensitive complexity: constraining an edge of the polygon has a cost proportional to the number of triangles which must be removed.

The triangulation should include a special vertex called the infinite vertex, which easily allows the convex hull of the vertices of the triangulation to be found. This will help to determine if a located point q is inside or outside P.

Step 2: Tagging Triangles of the Triangulation as Interior or Exterior.

This step is necessary to determine whether a vertex is inside or outside P.

• First find a vertex of P belonging to an “infinite” edge (i.e., the edge includes the “infinite” vertex)

• From this vertex, find an “infinite” triangle, tag it as outside P.

• Iterate over the triangles which share a common edge with this triangle, such that triangles contain either an edge of P or two vertices belonging to P

  • Also tag these triangles as outside P; they are the triangles comprising the non-convex vertices of P. This operation should be O(n) since triangles have only three neighbour triangles, and always have a vertex of P as vertex

Step 3: Locating and Computing Distances

• For each vertex q of Q, check whether it is contained in one of the triangles tagged as outside:
  • If q is not contained within any of the triangles the algorithm terminates: the distance is 0, since P and Q intersect
• If q is contained within one of the triangles:
  • Compute the distance between q and each edge and each vertex of the triangle, storing the minimum distance.

The minimum distance between vertices of Q and P has now been calculated.

Repeat the operation, reversing the role of the polygons:

• The triangulation of Q is computed and the closest vertex p of P to Q is calculated.
• The minimum Polygon distance is the minimum distance between this operation and the minimum distance found in Step 3.

Analysis: Triangulation-Based Polygon Distance Algorithm

Summary

We evaluated a triangulation-based approach for computing minimum distance between two non-convex polygons, as an alternative to the existing separable_geometry_distance_fast() projection-based algorithm.

Conclusion: The triangulation approach is theoretically sound but unlikely to outperform the existing implementation for one-off polygon distance queries. The analysis is documented here for future reference.

The Proposed Algorithm

Adapted from Amato, Nancy M (1992) and a Quora answer by Tom Dreyfus.

Key Insight

The minimum distance between two polygons must involve at least one vertex:

• If nearest edges are parallel, slide along until you hit a vertex
• If not parallel, the contact point is necessarily at a vertex

Algorithm Outline

1. Build a Constrained Delaunay Triangulation (CDT) of polygon P
2. Tag triangles as interior/exterior via flood-fill from the outer face
3. For each vertex q of polygon Q:
   • Locate q in the CDT
   • If q is in an interior triangle: distance = 0 (intersection)
   • If q is in an exterior triangle: find nearest constraint edge
4. Repeat with P and Q swapped
5. Return minimum distance found

Issues Identified

Problem with Original Step 2 Description

The original algorithm's description of interior/exterior tagging was incomplete. The correct approach:

// Flood-fill from outer face, stopping at constraint edges
fn tag_faces(cdt: &ConstrainedDelaunayTriangulation) -> HashSet<FaceHandle> {
    let mut exterior = HashSet::new();
    let mut queue = VecDeque::new();

    queue.push_back(cdt.outer_face());
    exterior.insert(cdt.outer_face());

    while let Some(face) = queue.pop_front() {
        for edge in face.adjacent_edges() {
            if !edge.is_constraint_edge() {
                let adjacent = edge.rev().face();
                if exterior.insert(adjacent) {
                    queue.push_back(adjacent);
                }
            }
        }
    }
    exterior  // All faces NOT in this set are interior
}

Spade provides the necessary API: outer_face(), is_constraint_edge(), adjacent_edges(), rev().face().

Problem with Step 3: Nearest Constraint Edge

For concave polygons, the nearest constraint edge to a point q in an exterior triangle may not be an edge of that triangle:

    +----+   +----+
    |    |   |    |
    |    | q |    |     <- q is IN the notch (exterior to P, inside convex hull)
    |    |   |    |
    |    +---+    |
    |             |
    +-------------+
         P (U-shape opening upward)

The CDT triangle containing q might have vertices from the top corners and some vertex below. Depending on q's exact position, the nearest point on P could be on the left inner wall, right inner wall, or bottom of the notch - and the nearest constraint edge might not be an edge of q's current triangle.

Solution: Walk-and-prune through exterior triangles:

fn nearest_constraint_distance(q: Coord, start_face: FaceHandle, cdt: &CDT) -> f64 {
    let mut best = f64::INFINITY;
    let mut queue = BinaryHeap::new();  // min-heap by lower_bound
    queue.push((Reverse(0.0), start_face));

    while let Some((Reverse(lower_bound), face)) = queue.pop() {
        if lower_bound >= best { break; }

        for edge in face.adjacent_edges() {
            if edge.is_constraint_edge() {
                best = best.min(distance(q, edge));
            } else {
                let edge_dist = distance(q, edge);
                if edge_dist < best {
                    queue.push((Reverse(edge_dist), edge.rev().face()));
                }
            }
        }
    }
    best
}

Complexity Comparison

Triangulation + Walk-and-Prune

Phase	Complexity
Build CDT of P	O(n log n)
Tag faces	O(n)
Per query: locate	O(log n)
Per query: walk	O(1) best, O(n) worst
Total for P-Q distance	O(n log n) + O(m * walk)

Walk cost depends on polygon shape and query point position. Worst case for complex concave polygons with query points far from constraint edges.

Existing: separable_geometry_distance_fast()

Phase	Complexity
Sort vertices by 1D projection	O((n+m) log(n+m))
Iterate with pruning	O((n+m) * k) where k = nearby pairs
Total	O((n+m) log(n+m)) + O((n+m) * k)

Best case k is small (well-separated geometries). Worst case k approaches n*m.

When Each Approach Wins

Triangulation might win when:

• Many queries against the same polygon P (amortise CDT construction)
• Queries are "near miss" cases (close to P, short walks)
• Need point-in-polygon as side effect

Projection wins when:

• One-off distance queries (no setup amortisation)
• Polygons are well-separated (small k, good pruning)
• Cache-friendly sorted array traversal

Why We Chose Not to Implement

1. separable_geometry_distance_fast() is only used when bounding boxes don't overlap, meaning polygons are guaranteed to be separated - exactly where projection excels

2. The walk-and-prune worst case (O(n) faces visited) is unpredictable and depends on polygon shape

3. CDT construction has non-trivial overhead for one-off queries

4. The existing implementation is already well-optimised with R-tree fallback for overlapping bboxes

References

• Amato, Nancy M (1992)
• Spade documentation
• Existing implementation: geo/src/algorithm/line_measures/metric_spaces/euclidean/distance.rs