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(N2) 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 p1, p2, ... and so on up to pn. 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 p1 ... pn given a tolerance t"

  1. Find the point whose perpendicular distance from the line segment connecting  p1 ... pn is the greatest.  Call the point pb and this maximal distance  db
  2. If db < t,
    • then throw away all of the points p2...pn-1
    • else
      1. Simplify the array of n points p1 ... pb given tolerance t
      2. Simplify the array of points pb..pn 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;
     vector = point;

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

var v1p:   vector;
    slash: real;

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)
           else begin
                slash := v1p.x * v12.y - v1p.y * v12.x;
                pt_to_seg_dist := abs (slash / sqrt (m12))


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

  procedure keep (i: integer);

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

  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
                di := pt_to_seg_dist (p [first], vfl, p [i]);
                if di < db
                   then begin
                        b := i;
                        db := di
                i := i + 1
           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;  { simplify_part }

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

Log in or register to write something here or to contact authors.