Solved

# Line generalization algorithm

Posted on 2003-11-24

Dear Sir/Madam,

I need to write a program in delphi to do line generalization. Douglas Peucker Algorithm is adopted to generalize the line. There are 2 units for my program and they are 'main' and 'douglas peucker'. The algorithm works quite good but there is some problem for the 'main' unit.

At present, I can draw a line on the canva by using the mouse and press 'Simplify' to generalize the line. However, I wish to draw a line by the following formula instead of using mouse as an input but I failed at last. So I need help to write some codes to incorporate the formula to 'main' unit.

X:=-3+6*i/100 (where i is 1, 2, 3......)

Y:=cos(X)+0.3*sin(4*X)+random(10)/100

The source of 'main' unit and 'douglas peucker' unit are as follows:

1. 'Main' Unit

unit Main;

interface

uses

Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms, Dialogs,

StdCtrls, ExtCtrls, Contnrs, DouglasPeuckers, OpenGL;

type

TForm1 = class(TForm)

pbMain: TPaintBox;

Label5: TLabel;

Label6: TLabel;

lbNumPtsOrig: TLabel;

Label7: TLabel;

lbNumPtsSimple: TLabel;

Button1: TButton;

procedure pbMainPaint(Sender: TObject);

// procedure pbMainMouseDown(Sender: TObject; Button: TMouseButton;

// Shift: TShiftState; X, Y: Integer);

procedure pbMainMouseMove(Sender: TObject; Shift: TShiftState; X,

Y: Integer);

procedure Button1Click(Sender: TObject);

private

{ Private declarations }

procedure AddPointToCurve(X, Y: integer);

procedure CreateSimplifiedPolyline;

public

OrigList: array of TPoint;

SimpleList: array of TPoint;

end;

var

Form1: TForm1;

implementation

{$R *.DFM}

procedure TForm1.AddPointToCurve(X, Y: integer);

var

APoint: TPoint;

//i:integer;

begin

//randomize;

//PosX := X;

//PosY := Y;

APoint.X := X;

APoint.Y := Y;

SetLength(OrigList, Length(OrigList) + 1);

OrigList[Length(OrigList) - 1] := APoint;

end;

procedure TForm1.pbMainPaint(Sender: TObject);

begin

with pbMain.Canvas do begin

// Draw original polyline

if Length(OrigList) > 0 then begin

Pen.Color := clBlack;

Pen.Width := 1;

PolyLine(OrigList);

end;

// Draw simplification

if Length(SimpleList) > 0 then begin

Pen.Color := clRed;

Pen.Width := round(1);

PolyLine(SimpleList);

end;

end;

// Other controls

lbNumPtsOrig.Caption := IntToStr(Length(OrigList));

lbNumPtsSimple.Caption := IntToStr(Length(SimpleList));

end;

procedure TForm1.pbMainMouseMove(Sender: TObject; Shift: TShiftState; X,

Y: Integer);

begin

// Store the new point

AddPointToCurve(X, Y);

pbMain.Invalidate;

end;

//end;

//end;

procedure TForm1.CreateSimplifiedPolyline;

// Create the simple polyline approximation

var

ALength: integer;

s:string;

r:integer;

begin

s:=InputBox('Input Box','Please Input the Tolerance','0');

r:=strtoint(s);

// Create the simple polyline approximation

SetLength(SimpleList, Length(OrigList));

if length(OrigList) > 2 then begin

ALength := PolySimplifyInt2D(r, OrigList, SimpleList);

SetLength(SimpleList, ALength);

end;

end;

procedure TForm1.Button1Click(Sender: TObject);

begin

// the button is pressed to draw the simplified line

CreateSimplifiedPolyline;

pbMain.Invalidate;

end;

end.

2. 'Douglas Peucker' Unit

unit DouglasPeuckers;

interface

uses

OpenGL, Windows;// We use TPoint from the windows unit

type

// Generalized float and int types

TFloat = double;

{ 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;

implementation

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;

end;

function PolySimplifyInt2D(Tol: TFloat; const Orig: array of TPoint;

var Simple: array of TPoint): integer;

var

i, N: integer;

Marker: array of boolean;

Tol2: TFloat;

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;

end;

end;

end.

Many thanks

Regards,

Kenneth