Douglas-Peucker Algorithm

A classic algorithm used for line simplification.   This algorithm was introduced to the world by David H. Douglas and Thomas K. Peucker in a 1973 article in the journal Canadian Cartographer.  It is by no means the fastest algorithm, taking O(n log N) time at best and O(N²) at worst.  However, since researchers at the University of Kansas have shown that this algorithm gives the closest approximation to the points a human being would choose when manually simplifying a line, it is probably the "best" in terms of the results obtained.

You start with an array of points.  Let's call them p₁, p₂, ... and so on up to p_n. n being the number of points. You also have a tolerance distance.  The algorithm is recursive, and essentially goes like this:

"Simplify an array of n points p₁ ... p_n given a tolerance t"
 

1. Find the point whose perpendicular distance from the line segment connecting  p₁ ... p_n is the greatest.  Call the point p_b and this maximal distance  d_b
2. If d_b < t,
• then throw away all of the points p₂...p_n-1
• else
1. Simplify the array of n points p₁ ... p_b given tolerance t
2. Simplify the array of points p_b..p_n given tolerance t.

The algorithm can be implemented with recursion or not, although it requires a stack without recursion.

You can use an array of points or a linked list.

---

In the tradition of noding more Pascal,  here is an implementation of the algorithm in that language, written for clarity rather than efficency:
 

type point = record
             x: real;
             y: real;
             end;
 
     vector = point;
 
 
 

function pt_to_seg_dist (var p1: point; var v12: vector; var p: point): real;

var v1p:   vector;
    m12,
    dot,
    slash: real;

begin
m12 := v12.x * v12.x + v12.y * v12.y;
v1p.x := p.x - p1.x;
v1p.y := p.y - p1.y;
dot := v1p.x * v12.x + v1p.y * v12.y;
if (dot <= 0.0)
   then pt_to_seg_dist := sqrt (v1p.x * v1p.x + v1p.y * v1p.y)
   else if dot >= m12
           then begin
                v1p.x := v1p.x + v12.x;
                v1p.y := v1p.y + v12.y;
                pt_to_seg_dist := sqrt (v1p.x * v1p.x + v1p.y * v1p.y)
                end
           else begin
                slash := v1p.x * v12.y - v1p.y * v12.x;
                pt_to_seg_dist := abs (slash / sqrt (m12))
                end
end;
 

 

function simplify (var p:         array [lo..hi: integer] of point;
                       ct:        integer;
                       tolerance: real): integer;
var kept: integer;
 

  procedure keep (i: integer);

  begin
  kept := kept + 1;
  if kept < i
     then begin
          p [kept].x := p [i].x;
          p [kept].y := p [i].y
          end
  end;
 

  procedure simplify_part (first: integer;
                           last:  integer);
  var i:  integer;
      b:  integer;
      di: real;
      db: real;
      vfl: point;

  begin { simplify_part }
  if last > first + 1
     then begin
          vfl.x := p [last].x - p [first].x;
          vfl.y := p [last].y - p [first].y;
          {
           find the intermediate point
           furthest from the segment
           connecting first and last
          }
          b := first+1;
          db := pt_to_seg_dist (p [first], vfl, p [b]);
          i := b + 1;
          while i < last do
                begin
                di := pt_to_seg_dist (p [first], vfl, p [i]);
                if di < db
                   then begin
                        b := i;
                        db := di
                        end
                i := i + 1
                end;
          {
           if the furthest distance beats the tolerance,
           recursively simplify the rest of the array.
          }
          if db >= tolerance
             then begin
                  simplify_part (first, b);
                  keep (b);
                  simplify_part (b, last)
                  end
          end
  end;  { simplify_part }
 

begin { simplify }
kept := 0;
keep (lo);
simplify_part (lo, ct);
keep (ct);
simplify := kept
end;  { simplify }