douglaspeuckers.pas

计算机图形学

unit DouglasPeuckers;

{ Implementation of the famous Douglas-Peucker polyline simplification
  algorithm.

  This file contains a 3D floating point implementation, for spatial
  polylines, as well as a 2D integer implementation for use with
  Windows GDI.

  Loosely based on C code from SoftSurfer (www.softsurfer.com)
  http://geometryalgorithms.com/Archive/algorithm_0205/algorithm_0205.htm

  References:
  David Douglas & Thomas Peucker, "Algorithms for the reduction of the number of
  points required to represent a digitized line or its caricature", The Canadian
  Cartographer 10(2), 112-122  (1973)

  Delphi code by Nils Haeck (c) 2003 Simdesign (www.simdesign.nl)
  http://www.simdesign.nl/components/douglaspeucker.html

  ****************************************************************
  The contents of this file are subject to the Mozilla Public
  License Version 1.1 (the "License"); you may not use this file
  except in compliance with the License. You may obtain a copy of
  the License at:
  http://www.mozilla.org/MPL/

  Software distributed under the License is distributed on an
  "AS IS" basis, WITHOUT WARRANTY OF ANY KIND, either express or
  implied. See the License for the specific language governing
  rights and limitations under the License.
}

interface

uses
     Windows; // We use TPoint from the windows unit

type

     // Generalized float and int types
     TFloat = Double;

     // Float point 3D
     TPointFloat3D = packed record
          X: TFloat;
          Y: TFloat;
          Z: TFloat;
     end;

     { PolySimplifyFloat3D:
       Approximates the polyline with 3D float vertices in Orig, with a simplified
       version that will be returned in Simple. The maximum deviation from the
       original line is given in Tol.
       Input:  Tol      = approximation tolerance
               Orig[]   = polyline array of vertex points
       Output: Simple[] = simplified polyline vertices. This array must initially
                          have the same length as Orig
       Return: the number of points in Simple
     }

function PolySimplifyFloat3D(Tol: TFloat; const Orig: array of TPointFloat3D;
     var Simple: array of TPointFloat3D): Integer;

{ PolySimplifyInt2D:
  Approximates the polyline with 2D integer vertices in Orig, with a simplified
  version that will be returned in Simple. The maximum deviation from the
  original line is given in Tol.
  Input:  Tol      = approximation tolerance
          Orig[]   = polyline array of vertex points
  Output: Simple[] = simplified polyline vertices. This array must initially
                     have the same length as Orig
  Return: the number of points in Simple
}

function PolySimplifyInt2D(Tol: TFloat; const Orig: array of TPoint;
     var Simple: array of TPoint): Integer;

{ MinDistPointLine calculates the minimum distance of a point to a line.
  P is the point, the line is between points A and B.
}

function MinDistPointLine(Px, Py, Ax, Ay, Bx, By: Double): Double;

implementation

function VecMinFloat3D(const A, B: TPointFloat3D): TPointFloat3D;
// Result = A - B
begin
     Result.X := A.X - B.X;
     Result.Y := A.Y - B.Y;
     Result.Z := A.Z - B.Z;
end;

function DotProdFloat3D(const A, B: TPointFloat3D): TFloat;
// Dotproduct = A * B
begin
     Result := A.X * B.X + A.Y * B.Y + A.Z * B.Z;
end;

function NormSquaredFloat3D(const A: TPointFloat3D): TFloat;
// Square of the norm |A|
begin
     Result := A.X * A.X + A.Y * A.Y + A.Z * A.Z;
end;

function DistSquaredFloat3D(const A, B: TPointFloat3D): TFloat;
// Square of the distance from A to B
begin
     Result := NormSquaredFloat3D(VecMinFloat3D(A, B));
end;

procedure SimplifyFloat3D(var Tol2: TFloat; const Orig: array of TPointFloat3D;
     var Marker: array of boolean; j, k: integer);
// Simplify polyline in OrigList between j and k. Marker[] will be set to True
// for each point that must be included
var
     i, MaxI: integer; // Index at maximum value
     MaxD2: TFloat; // Maximum value squared
     CU, CW, B: TFloat;
     DV2: TFloat;
     P0, P1, PB, U, W: TPointFloat3D;
begin
     // Is there anything to simplify?
     if k <= j + 1 then exit;

     P0 := Orig[j];
     P1 := Orig[k];

     U := VecMinFloat3D(P1, P0); // Segment vector
     CU := DotProdFloat3d(U, U); // Segment length squared

     MaxD2 := 0;
     MaxI := 0;

     // Loop through points and detect the one furthest away
     for i := j + 1 to k - 1 do
     begin
          W := VecMinFloat3D(Orig[i], P0);
          CW := DotProdFloat3D(W, U);

          // Distance of point Orig[i] from segment
          if CW <= 0 then
          begin
               // Before segment
               DV2 := DistSquaredFloat3D(Orig[i], P0)
          end
          else
          begin
               if CW > CU then
               begin
                    // Past segment
                    DV2 := DistSquaredFloat3D(Orig[i], P1);
               end
               else
               begin
                    // Fraction of the segment
                    try
                         B := CW / CU;
                    except
                         B := 0; // in case CU = 0
                    end; //try

                    PB.X := P0.X + B * U.X;
                    PB.Y := P0.Y + B * U.Y;
                    PB.Z := P0.Z + B * U.Z;
                    DV2 := DistSquaredFloat3D(Orig[i], PB);
               end; //
          end;

          // test with current max distance squared
          if DV2 > MaxD2 then
          begin
               // Orig[i] is a new max vertex
               MaxI := i;
               MaxD2 := DV2;
          end; //if
     end; //for

     // If the furthest point is outside tolerance we must split
     if MaxD2 > Tol2 then
     begin // error is worse than the tolerance

          // split the polyline at the farthest vertex from S
          Marker[MaxI] := True; // mark Orig[maxi] for the simplified polyline

          // recursively simplify the two subpolylines at Orig[maxi]
          SimplifyFloat3D(Tol2, Orig, Marker, j, MaxI); // polyline Orig[j] to Orig[maxi]
          SimplifyFloat3D(Tol2, Orig, Marker, MaxI, k); // polyline Orig[maxi] to Orig[k]
     end; //if
end;

function VecMinInt2D(const A, B: TPoint): TPoint;
// Result = A - B
begin
     Result.X := A.X - B.X;
     Result.Y := A.Y - B.Y;
end;

function DotProdInt2D(const A, B: TPoint): TFloat;
// Dotproduct = A * B
begin
     Result := A.X * B.X + A.Y * B.Y;
end;

function NormSquaredInt2D(const A: TPoint): TFloat;
// Square of the norm |A|
begin
     Result := A.X * A.X + A.Y * A.Y;
end;

function DistSquaredInt2D(const A, B: TPoint): TFloat;
// Square of the distance from A to B
begin
     Result := NormSquaredInt2D(VecMinInt2D(A, B));
end;

procedure SimplifyInt2D(var Tol2: TFloat; const Orig: array of TPoint;
     var Marker: array of Boolean; J, K: Integer);
// Simplify polyline in OrigList between j and k. Marker[] will be set to True
// for each point that must be included
var
     I, MaxI: integer; // Index at maximum value
     MaxD2: TFloat; // Maximum value squared
     CU, CW, B: TFloat;
     DV2: TFloat;
     P0, P1, PB, U, W: TPoint;
begin
     // Is there anything to simplify?
     if K <= J + 1 then Exit;

     P0 := Orig[J];
     P1 := Orig[K];
     U := VecMinInt2D(P1, P0); // Segment vector
     CU := DotProdInt2D(U, U); // Segment length squared

     MaxD2 := 0;
     MaxI := 0;

     // Loop through points and detect the one furthest away
     for I := J + 1 to K - 1 do
     begin
          W := VecMinInt2D(Orig[I], P0);
          CW := DotProdInt2D(W, U);

          // Distance of point Orig[i] from segment
          if CW <= 0 then
          begin
               // Before segment
               DV2 := DistSquaredInt2D(Orig[I], P0)
          end
          else
          begin
               if CW > CU then
               begin
                    // Past segment
                    DV2 := DistSquaredInt2D(Orig[I], P1);
               end
               else
               begin
                    // Fraction of the segment
                    try
                         B := CW / CU;
                    except
                         B := 0; // in case CU = 0
                    end;

                    PB.X := round(P0.X + B * U.X);
                    PB.Y := round(P0.Y + B * U.Y);
                    DV2 := DistSquaredInt2D(Orig[I], PB);
               end;
          end;

          // test with current max distance squared
          if DV2 > MaxD2 then
          begin
               // Orig[i] is a new max vertex
               MaxI := I;
               MaxD2 := DV2;
          end;
     end;

     // If the furthest point is outside tolerance we must split
     if MaxD2 > Tol2 then
     begin // error is worse than the tolerance

          // split the polyline at the farthest vertex from S
          Marker[MaxI] := True; // mark Orig[maxi] for the simplified polyline

          // recursively simplify the two subpolylines at Orig[maxi]
          SimplifyInt2D(Tol2, Orig, Marker, J, MaxI); // polyline Orig[j] to Orig[maxi]
          SimplifyInt2D(Tol2, Orig, Marker, MaxI, K); // polyline Orig[maxi] to Orig[k]
     end; //if
end;

function PolySimplifyFloat3D(Tol: TFloat; const Orig: array of TPointFloat3D;
     var Simple: array of TPointFloat3D): integer;
var
     i, N: integer;
     Tol2: TFloat;
     Marker: array of boolean;
begin
     Result := 0;
     if length(Orig) < 2 then exit;
     Tol2 := sqr(Tol);

     // Create a marker array
     N := Length(Orig);
     SetLength(Marker, N);
     // Include first and last point
     Marker[0] := True;
     Marker[N - 1] := True;
     // Exclude intermediate for now
     for i := 1 to N - 2 do
          Marker[i] := False;

     // Simplify
     SimplifyFloat3D(Tol2, Orig, Marker, 0, N - 1);

     // Copy to resulting list
     for i := 0 to N - 1 do
     begin
          if Marker[i] then
          begin
               Simple[Result] := Orig[i];
               inc(Result);
          end;
     end;
end;

function PolySimplifyInt2D(Tol: TFloat; const Orig: array of TPoint;
     var Simple: array of TPoint): Integer;
var
     I, N: Integer;
     Tol2: TFloat;
     Marker: array of Boolean;
begin
     Result := 0;

     if Length(Orig) < 2 then exit;

     Tol2 := sqr(Tol);
     // Create a marker array
     N := Length(Orig);
     SetLength(Marker, N);
     // Include first and last point
     Marker[0] := True;
     Marker[N - 1] := True;
     // Exclude intermediate for now
     for i := 1 to N - 2 do
          Marker[i] := False;

     // Simplify
     SimplifyInt2D(Tol2, Orig, Marker, 0, N - 1);

     // Copy to resulting list
     for i := 0 to N - 1 do
     begin
          if Marker[i] then
          begin
               Simple[Result] := Orig[i];
               Inc(Result);
          end; //if
     end; //for
end;

{
Tip #1: Minimum distance between a point and a line
Added: 15Nov2002  Author: Nils Haeck  Category: Geometry

Question:
How can I calculate the distance between a point and a line?

Applicable:
You can use this code when you need to detect whether the mouse click of the user is near or on a line segment or not, like in GetHitTestInfo events.

Answer:
We'll use some optimisation theory: the minimum distance between the point and the line
(which can be expressed as a function) will be at the exact location where the derivative of this
function is zero.
Perhaps you remember this from school; the minimum or maximum in a parabola is where its gradient
is zero.

We parametrise the distance of point P to the line between A and B as the distance of point P to a
point Q on the line:

point Q = (1-q)A+qB where 0 <= q <= 1

The distance PQ is:

|PQ| = sqrt( ((1-q)Ax + qBx - Px)^2 + (... Y term) )

Differentiating gives dPQ/dq = 2((Bx-Ax)q + (Ax-Px))(Bx - Ax) + (... Y term).

dPQ/dq must be zero for minimum so:

q = (Px-Ax)(Bx-Ax)+(Py-Ay)(By-Ay) / ((Bx-Ax)^2+(By-Ay)^2)

Code Sample:
Note that this code also takes into account situations where your line is actually not a line (A=B)
and situations where the point P is past any of the two endpoints. In this case,
the distance to the closest endpoint is calculated.
}

// MinDistPointLine calculates the minimum distance of a point to a line.
// P is the point, the line is between points A and B.

function PointToPointDist(Ax, Ay, Bx, By: Double): Double;
begin
     Result := Sqrt(Sqr(Bx - Ax) + Sqr(By - Ay));
end;

function MinDistPointLine(Px, Py, Ax, Ay, Bx, By: Double): Double;
var
     q: Double;
begin
     if (Ax = Bx) and (Ay = By) then
     begin
          // Point to point
          Result := PointToPointDist(Px, Py, Ax, Ay);
     end
     else
     begin
          // Minimum
          q := ((Px - Ax) * (Bx - Ax) + (Py - Ay) * (By - Ay)) / (Sqr(Bx - Ax) + Sqr(By - Ay));

          // Limit q to 0 <= q <= 1
          if q < 0 then q := 0;

          if q > 1 then q := 1;

          // Distance
          Result := PointToPointDist(Px, Py, (1 - q) * Ax + q * Bx, (1 - q) * Ay + q * By);
     end; //if
end;

end.