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^{2}) 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_{1},
p_{2},
... 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_{1}
... p_{n} given a tolerance t"

Find the point whose perpendicular distance from the line segment connecting
p_{1} ... p_{n} is the
greatest. Call the point p_{b}_{ }and
this maximal distance d_{b}_{}

If d_{b }< t,

then throw away all of the points p_{2}...p_{n1}

else

Simplify the array of n points p_{1}
... p_{b} given tolerance t

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 }