Jim Chen's Blog
0%

Topology-Preserving RDP

Posted on Edited on

The Douglas-Peucker (RDP) algorithm is excellent for simplifying polylines by removing points within a tolerance epsilon. However, standard RDP only guarantees geometric proximity — it doesn’t prevent topological violations like self-intersections, collapsed holes, or accidental crossings between nearby features.

In this post, we’ll explore what topology-preserving simplification means and how to implement it using RhinoCommon.

Why Standard RDP Can Break Topology

While RDP with epsilon controls point deviation, it can still create invalid geometries:

Common topology violations:

  • Self-intersections: A simplified polyline crosses itself
  • Collapsed holes: Interior vertices of a polygon are removed, collapsing a hole
  • Feature crossings: Two close parallel lines collapse into crossing lines after simplification
  • Lost connectivity: Road intersections disappear or new false intersections appear

These issues arise because RDP makes local decisions (remove this point based on distance) without considering global constraints (does this create an intersection?).

What “Topology-Preserving” Actually Means

In spatial geometry, topology refers to relationships between features:

Relationship Example
Adjacency Two polygons sharing a border
Containment A hole inside a polygon
Non-intersection Polylines that never cross
Connectivity Road segments meeting at intersections

A topology-preserving simplification ensures:

✔ Polygon rings stay valid (no self-intersections)
✔ Holes remain inside boundary shells
✔ Shared boundaries aren’t broken or duplicated
✔ Lines that didn’t intersect before still don’t intersect after
✔ Connectivity is not accidentally removed or added

This is stronger than just controlling point deviation with epsilon. Standard RDP only guarantees geometric proximity, not topological validity.

How Topology-Preserving RDP Works

Most topology-preserving algorithms extend RDP with additional checks:

1. Self-Intersection Checks

Before removing a point, test whether the new segment would intersect any existing segments in the polyline.

1
2
3
4
Original:  A─B─C─D─E
           └───┘ 
Remove C:  A───┬─D─E  ← Does A-D intersect B-C? If yes, keep C.
               B

2. Constraint Propagation

When simplifying adjacent features (e.g., shared polygon borders), ensure both sides simplify consistently to maintain adjacency.

3. Global Validity Guarantees

After simplification, verify the output is a valid geometry (e.g., closed polygons without crossing boundaries).

Implementation in RhinoCommon

Let’s implement a topology-preserving RDP algorithm for polylines in C#/RhinoCommon.

Core Algorithm: Standard RDP

First, here’s the standard RDP as a baseline:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
using Rhino.Geometry;
using System;
using System.Collections.Generic;
using System.Linq;

public static class PolylineSimplification
{
    /// <summary>
    /// Standard Douglas-Peucker simplification
    /// </summary>
    public static List<Point3d> SimplifyRDP(List<Point3d> points, double epsilon)
    {
        if (points.Count < 3)
            return new List<Point3d>(points);

        // Find the point with maximum distance from line segment
        double maxDist = 0;
        int maxIndex = 0;
        
        Point3d start = points[0];
        Point3d end = points[points.Count - 1];
        Line baseLine = new Line(start, end);

        for (int i = 1; i < points.Count - 1; i++)
        {
            double dist = baseLine.DistanceTo(points[i], false);
            if (dist > maxDist)
            {
                maxDist = dist;
                maxIndex = i;
            }
        }

        // If max distance is greater than epsilon, recursively simplify
        if (maxDist > epsilon)
        {
            // Recursive call for left and right segments
            List<Point3d> left = SimplifyRDP(points.GetRange(0, maxIndex + 1), epsilon);
            List<Point3d> right = SimplifyRDP(points.GetRange(maxIndex, points.Count - maxIndex), epsilon);

            // Combine results (remove duplicate middle point)
            left.RemoveAt(left.Count - 1);
            left.AddRange(right);
            return left;
        }
        else
        {
            // All points can be removed except endpoints
            return new List<Point3d> { start, end };
        }
    }
}

Topology-Preserving Extension

Now let’s add self-intersection checks:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
public static class TopologyPreservingSimplification
{
    /// <summary>
    /// Topology-preserving Douglas-Peucker simplification
    /// Prevents self-intersections in the simplified polyline
    /// </summary>
    public static List<Point3d> SimplifyPreserveTopology(List<Point3d> points, double epsilon)
    {
        if (points.Count < 3)
            return new List<Point3d>(points);

        return SimplifyRecursive(points, epsilon, 0, points.Count - 1);
    }

    private static List<Point3d> SimplifyRecursive(List<Point3d> points, double epsilon, 
        int startIdx, int endIdx)
    {
        if (endIdx - startIdx < 2)
            return new List<Point3d> { points[startIdx], points[endIdx] };

        // Find point with maximum distance
        double maxDist = 0;
        int maxIndex = startIdx;
        
        Point3d start = points[startIdx];
        Point3d end = points[endIdx];
        Line baseLine = new Line(start, end);

        for (int i = startIdx + 1; i < endIdx; i++)
        {
            double dist = baseLine.DistanceTo(points[i], false);
            if (dist > maxDist)
            {
                maxDist = dist;
                maxIndex = i;
            }
        }

        // Check if we can safely simplify (remove all intermediate points)
        if (maxDist <= epsilon)
        {
            // TOPOLOGY CHECK: Would creating segment start->end cause self-intersection?
            Line newSegment = new Line(start, end);
            
            if (WouldCauseSelfIntersection(points, startIdx, endIdx, newSegment))
            {
                // Cannot simplify - keep the furthest point and recurse
                List<Point3d> left = SimplifyRecursive(points, epsilon, startIdx, maxIndex);
                List<Point3d> right = SimplifyRecursive(points, epsilon, maxIndex, endIdx);
                
                left.RemoveAt(left.Count - 1);
                left.AddRange(right);
                return left;
            }
            else
            {
                // Safe to simplify
                return new List<Point3d> { start, end };
            }
        }
        else
        {
            // Recursively simplify both sides
            List<Point3d> left = SimplifyRecursive(points, epsilon, startIdx, maxIndex);
            List<Point3d> right = SimplifyRecursive(points, epsilon, maxIndex, endIdx);
            
            left.RemoveAt(left.Count - 1);
            left.AddRange(right);
            return left;
        }
    }

    /// <summary>
    /// Check if a new segment would cause self-intersection with the polyline
    /// </summary>
    private static bool WouldCauseSelfIntersection(List<Point3d> points, 
        int startIdx, int endIdx, Line newSegment)
    {
        // Check all existing segments in the polyline
        for (int i = 0; i < points.Count - 1; i++)
        {
            // Skip adjacent segments (they naturally share endpoints)
            if (i == startIdx - 1 || i == startIdx || i == endIdx - 1 || i == endIdx)
                continue;

            // Skip if we're checking the segment being replaced
            if (i >= startIdx && i < endIdx)
                continue;

            Line existingSegment = new Line(points[i], points[i + 1]);
            
            // Check for intersection
            double a, b;
            if (Rhino.Geometry.Intersect.Intersection.LineLine(
                newSegment, existingSegment, out a, out b, 0.001, false))
            {
                // Check if intersection is at interior of both segments (not just endpoints)
                if (a > 0.001 && a < 0.999 && b > 0.001 && b < 0.999)
                {
                    return true; // Self-intersection detected
                }
            }
        }
        
        return false; // No self-intersection
    }

    /// <summary>
    /// Simplify a polyline curve with topology preservation
    /// </summary>
    public static Polyline SimplifyPolylineCurve(Polyline polyline, double epsilon)
    {
        bool isClosed = polyline.IsClosed;
        List<Point3d> points = polyline.ToList();
        
        // For closed polylines, temporarily open them
        if (isClosed && points.Count > 0)
        {
            points.RemoveAt(points.Count - 1); // Remove duplicate closing point
        }
        
        List<Point3d> simplified = SimplifyPreserveTopology(points, epsilon);
        
        // Restore closed state if needed
        if (isClosed && simplified.Count > 0)
        {
            simplified.Add(simplified[0]); // Re-close the polyline
        }
        
        return new Polyline(simplified);
    }
}

Handling Closed Polylines

For closed polylines (polygons), we need additional care:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
public static class PolygonSimplification
{
    /// <summary>
    /// Simplify a closed polygon with topology preservation
    /// </summary>
    public static Polyline SimplifyPolygon(Polyline polygon, double epsilon)
    {
        if (!polygon.IsClosed || polygon.Count < 4)
            return polygon;

        List<Point3d> points = polygon.ToList();
        points.RemoveAt(points.Count - 1); // Remove duplicate closing point
        
        // Find the point furthest from centroid as a stable anchor
        Point3d centroid = AveragePoint(points);
        int anchorIdx = FindFurthestPoint(points, centroid);
        
        // Rotate list so anchor is first
        List<Point3d> rotated = RotateList(points, anchorIdx);
        
        // Simplify as open polyline
        List<Point3d> simplified = TopologyPreservingSimplification
            .SimplifyPreserveTopology(rotated, epsilon);
        
        // Check for self-intersections in the final closed polygon
        simplified.Add(simplified[0]); // Close it
        
        if (HasSelfIntersection(simplified))
        {
            // Fall back to original if simplification created invalid geometry
            return polygon;
        }
        
        return new Polyline(simplified);
    }

    private static Point3d AveragePoint(List<Point3d> points)
    {
        double x = 0, y = 0, z = 0;
        foreach (var pt in points)
        {
            x += pt.X;
            y += pt.Y;
            z += pt.Z;
        }
        int count = points.Count;
        return new Point3d(x / count, y / count, z / count);
    }

    private static int FindFurthestPoint(List<Point3d> points, Point3d reference)
    {
        double maxDist = 0;
        int maxIdx = 0;
        for (int i = 0; i < points.Count; i++)
        {
            double dist = reference.DistanceTo(points[i]);
            if (dist > maxDist)
            {
                maxDist = dist;
                maxIdx = i;
            }
        }
        return maxIdx;
    }

    private static List<Point3d> RotateList(List<Point3d> list, int startIdx)
    {
        var result = new List<Point3d>();
        for (int i = 0; i < list.Count; i++)
        {
            result.Add(list[(startIdx + i) % list.Count]);
        }
        return result;
    }

    private static bool HasSelfIntersection(List<Point3d> points)
    {
        for (int i = 0; i < points.Count - 1; i++)
        {
            Line seg1 = new Line(points[i], points[i + 1]);
            
            for (int j = i + 2; j < points.Count - 1; j++)
            {
                // Skip adjacent segments
                if (j == i + 1 || (i == 0 && j == points.Count - 2))
                    continue;
                
                Line seg2 = new Line(points[j], points[j + 1]);
                
                double a, b;
                if (Rhino.Geometry.Intersect.Intersection.LineLine(
                    seg1, seg2, out a, out b, 0.001, false))
                {
                    if (a > 0.001 && a < 0.999 && b > 0.001 && b < 0.999)
                        return true;
                }
            }
        }
        return false;
    }
}

Usage Example

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
// Example: Simplify a polyline with topology preservation
public void SimplifyPolylineExample()
{
    // Create a test polyline (potentially with close segments that might intersect)
    var points = new List<Point3d>
    {
        new Point3d(0, 0, 0),
        new Point3d(1, 0.1, 0),
        new Point3d(2, -0.1, 0),
        new Point3d(3, 0.2, 0),
        new Point3d(4, 0, 0),
        new Point3d(3.5, 2, 0),
        new Point3d(3, 4, 0),
        new Point3d(2, 3.9, 0),
        new Point3d(1, 4.1, 0),
        new Point3d(0, 4, 0)
    };
    
    Polyline polyline = new Polyline(points);
    
    double epsilon = 0.5;
    
    // Standard RDP (might create self-intersections)
    List<Point3d> standardSimplified = PolylineSimplification.SimplifyRDP(points, epsilon);
    
    // Topology-preserving RDP (prevents self-intersections)
    Polyline topoSimplified = TopologyPreservingSimplification
        .SimplifyPolylineCurve(polyline, epsilon);
    
    Rhino.RhinoApp.WriteLine($"Original points: {points.Count}");
    Rhino.RhinoApp.WriteLine($"Standard RDP: {standardSimplified.Count}");
    Rhino.RhinoApp.WriteLine($"Topology-preserving: {topoSimplified.Count}");
}

Performance Considerations

Topology-preserving simplification is slower than standard RDP because:

  • Each simplification decision requires intersection checks: $O(n)$ per check
  • Overall complexity: $O(n^2 \log n)$ vs $O(n \log n)$ for standard RDP

Optimization strategies:

  1. Spatial indexing: Use R-trees to speed up intersection queries
  2. Coarse filtering: Check bounding boxes before detailed intersection tests
  3. Progressive simplification: Apply aggressive epsilon first, then refine

Real-World Applications

Topology-preserving simplification is crucial in:

Domain Use Case
GIS Simplifying administrative boundaries while preserving adjacency
Cartography Map generalization at different zoom levels
CAD/BIM Reducing mesh complexity while maintaining valid geometry
Game Development LOD (Level of Detail) generation for terrain

References & Further Reading

  • PostGIS: ST_SimplifyPreserveTopology - PostgreSQL/PostGIS implementation
  • GEOS Library: C++ library used by PostGIS for geometric operations
  • Visvalingam-Whyatt: Alternative algorithm that can also be made topology-aware (see Visvalingam-Whyatt post)

Conclusion

While standard RDP is fast and simple, topology-preserving RDP is essential when geometric validity matters. By adding self-intersection checks, we ensure that simplified geometries remain valid and maintain their topological relationships.

The performance cost is worth it for applications where invalid geometry could cause downstream failures — such as spatial databases, GIS analysis, or manufacturing CAD models.

Key takeaway: Always consider whether your simplification needs to preserve topology. If yes, implement intersection checks. If performance is critical, invest in spatial indexing structures.