Visvalingam-Whyatt

Visvalingam-Whyatt is an algorithm similar to Douglas-Peucker algorithm for approximating a list of points (Polyline), but it often produces results that are visually more pleasing because it iteratively removes the “least important” features based on area, rather than distance.

Concept

While RDP removes points based on perpendicular distance to a line segment (ignoring local geometry), Visvalingam-Whyatt (VW) focuses on Effective Area.
The “Effective Area” of a point is defined as the area of the triangle formed by the point and its two immediate neighbors.

1. Calculate the triangular area for every internal point.
2. Find the point with the minimum effective area.
3. Remove that point.
4. Update the effective areas of the adjacent points (since their neighbors have changed).
5. Repeat until the desired number of points is left or the minimum area exceeds a tolerance.

Why Determinant?

The area of a triangle with vertices $A=(x_1, y_1)$, $B=(x_2, y_2)$, and $C=(x_3, y_3)$ can be calculated using the determinant of a matrix:

$$\text{Area} = \frac{1}{2} | x_1(y_2 - y_3) + x_2(y_3 - y_1) + x_3(y_1 - y_2) |$$

Graphically Meaning

Graphically, this is related to the Cross Product of two vectors. If we form vectors $\vec{AB}$ and $\vec{AC}$, their cross product represents the area of the parallelogram formed by them. The triangle area is exactly half of that parallelogram.

$$\text{Area} = \frac{1}{2} | \vec{AB} \times \vec{AC} |$$

In 2D (or projected on XY plane), this magnitude corresponds to the determinant formula above. Points that form a straight line have an area of 0 (collinear). Points that form a sharp spike have a small area (if the base is narrow). This metric ensures we keep points that contribute significantly to the shape’s overall character.

Implementation with Priority Queue

A naive implementation would scan the entire list to find the minimum area point in every iteration, leading to $O(N^2)$ complexity.
To make this efficient, we can use a Priority Queue.

• Initialization: Calculate areas for all points and insert them into a Min-Heap (Priority Queue).
• Loop:
  1. ExtractMin() to get the point with the smallest area ($O(\log N)$).
  2. Mark it as removed.
  3. Recalculate areas for its left and right neighbors.
  4. Update their positions in the Priority Queue ($O(\log N)$).

This reduces the complexity to $O(N \log N)$.

Code

Below is a C# implementation using RhinoCommon. Note that since .NET 6, PriorityQueue<TElement, TPriority> is built-in. For older versions (like Rhino 7/Classic), you might need a custom MinHeap class, but logic remains the same. Here I use SortedSet as a makeshift priority queue for simplicity.

using System;
using System.Collections.Generic;
using Rhino.Geometry;

public Polyline SimplifyVisvalingamWhyatt(Polyline pl, int targetCount)
{
    if (pl.Count <= targetCount || pl.Count < 3) return pl;

    int count = pl.Count;
    var areas = new double[count];
    var left = new int[count];
    var right = new int[count];
    
    // 1. Initialize Linked List structure and Areas
    for (int i = 0; i < count; i++)
    {
        left[i] = i - 1;
        right[i] = i + 1;
        // Endpoints have infinite area (never removed)
        areas[i] = (i == 0 || i == count - 1) ? double.MaxValue : CalculateTriangleArea(pl[i - 1], pl[i], pl[i + 1]);
    }

    // 2. Build Priority Queue
    // SortedSet stores tuples of (Area, Index). It sorts by Area, then Index.
    var minHeap = new SortedSet<(double area, int index)>();
    
    for(int i = 1; i < count - 1; i++)
    {
        minHeap.Add((areas[i], i));
    }

    int removedCount = 0;
    int toRemove = count - targetCount;

    while (removedCount < toRemove && minHeap.Count > 0)
    {
        // 3. Extract Min
        var minItem = minHeap.Min;
        minHeap.Remove(minItem);
        int idx = minItem.index;

        // "Remove" the point by updating neighbors
        int prev = left[idx];
        int next = right[idx];

        if (prev < 0 || next >= count) continue; // Boundary check

        // Update Links
        right[prev] = next;
        left[next] = prev;

        // 4. Update Neighbors in Heap
        // We must remove the old neighbor entries from heap before updating their keys
        if (prev > 0)
        {
            minHeap.Remove((areas[prev], prev));
            areas[prev] = CalculateTriangleArea(pl[left[prev]], pl[prev], pl[next]);
            minHeap.Add((areas[prev], prev));
        }

        if (next < count - 1)
        {
            minHeap.Remove((areas[next], next));
            areas[next] = CalculateTriangleArea(pl[prev], pl[next], pl[right[next]]);
            minHeap.Add((areas[next], next));
        }

        removedCount++;
    }

    // Reconstruct Polyline
    Polyline simplified = new Polyline();
    int current = 0;
    while(current < count && current != -1)
    {
        simplified.Add(pl[current]);
        current = right[current];
        if (current >= count) break;
    }

    return simplified;
}

private double CalculateTriangleArea(Point3d a, Point3d b, Point3d c)
{
    // Cross Product magnitude / 2 (Shoelace Formula equivalent)
    // Area = |(Bx - Ax)(Cy - Ay) - (By - Ay)(Cx - Ax)| / 2
    return 0.5 * Math.Abs((b.X - a.X) * (c.Y - a.Y) - (b.Y - a.Y) * (c.X - a.X));
}