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"
-
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
-
If db < t,
-
then throw away all of the points p2...pn-1
-
else
-
Simplify the array of n points p1
... pb given tolerance t
-
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;
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 }