unit polygonalyzation;

interface
uses
  image;

type
  TApproximationMethod = (amDouglasPeucker,amSklanskyGonzalez);

function polygonalization_SklanskyGonzalez( const segmentList : TSegmentList; const apprimationAccuracy:Single; const imageMonitoring : PBitmap32) : TSegmentList;
function polygonalization_DouglasPeucker( const segmentList : TSegmentList; const apprimationAccuracy:Single; const imageMonitoring : PBitmap32) : TSegmentList;


implementation
uses
  Math;

{******************************************************************************}
{polygonalization method from Sklansky and Gonzalez (1980)}
function polygonalization_SklanskyGonzalez( const segmentList : TSegmentList; const apprimationAccuracy:Single; const imageMonitoring : PBitmap32) : TSegmentList;
var
  i : Integer;
  lastCriticalPoint : TFPoint;
  wedge, lastWedge : TWedge;
  vectorWidth : Single;
  vectorWidthCount : Integer;
  segment : TSegment;
  appoximatedSegmentList : TSegmentList;
begin
  SetLength(appoximatedSegmentList,0);
  if Length(segmentList)>0 then begin
    // approximation
    vectorWidth:=0;
    vectorWidthCount:=0;
    lastCriticalPoint:=segmentList[0].p1;
    lastWedge.line1.a:=math.MaxSingle;
    // for each segment
    for i:=0 to Length(segmentList)-1 do begin
      segment:=segmentList[i];
      vectorWidth:=vectorWidth+segment.width;
      Inc(vectorWidthCount);
      // for the first segment, we only calculate the first wedge
      if i>0 then begin
        // for next segment, check if the new point is inside the last wedge
        if isPointInWedge(segment.p2,lastWedge)=true then begin
          // if the point is inside, then it isn't critical
        end else begin
          // if the point is outside, then the last point is critical
          vectorWidth:=vectorWidth / vectorWidthCount;
          addSegment(appoximatedSegmentList,lastCriticalPoint,segment.p1,vectorWidth);
          vectorWidth:=segment.width;
          vectorWidthCount:=1;
          lastCriticalPoint:=segment.p1;
        end;
      end;
      // calculate the wedge at this new point
      wedge:=getWedge(lastCriticalPoint,segment.p2,lastWedge,apprimationAccuracy,imageMonitoring);
      lastWedge:=wedge;
    end;
    // the last point is critical in all case
    vectorWidth:=vectorWidth / vectorWidthCount;
    addSegment(appoximatedSegmentList,lastCriticalPoint,segment.p2,vectorWidth);
  end;
  Result:=appoximatedSegmentList;
end;

{******************************************************************************}
{polygonalization method from Douglas-Peucker polyline simplification algorithm}
{adaptation of the implementation of "http://www.simdesign.nl/Components/DouglasPeucker.html"
 here their licence :
  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.
}
function VecMinFloat3D(const A, B: TFPoint): TFPoint;
// Result = A - B
begin
  Result.X := A.X - B.X;
  Result.Y := A.Y - B.Y;
end;

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

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

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

procedure SimplifyFloat3D(var Tol2: double;
  const Orig: array of TFPoint;
  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: double;    // Maximum value squared
  CU, CW, B: double;
  DV2: double;
  P0, P1, PB, U, W: TFPoint;
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;
        PB.X := P0.X + B * U.X;
        PB.Y := P0.Y + B * U.Y;
        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;
  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]
    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;
end;

function PolySimplifyFloat3D(Tol: double;
    const Orig: array of TFPoint;
    var Simple: array of TFPoint): integer;
var
  i, N: integer;
  Marker: array of boolean;
  Tol2: double;
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 polygonalization_DouglasPeucker( const segmentList : TSegmentList; const apprimationAccuracy:Single; const imageMonitoring : PBitmap32) : TSegmentList;
var
  originalPointArray, simplifiedPointArray : array of TFPoint;
  s, p : integer;
  simplifiedPointArrayLength : integer;
  segmentListLength : integer;
begin
  segmentListLength:=Length(segmentList);
  if segmentListLength>0 then begin
    // convert a list of segment to a list of point
    SetLength(originalPointArray,segmentListLength+1);
    p:=0;
    for s:=0 to segmentListLength-1 do begin
      originalPointArray[p]:=segmentList[s].p1;
      Inc(p);
    end;
    originalPointArray[p]:=segmentList[segmentListLength-1].p2;
    // call implementation of Douglas-Peucker algorithm
    SetLength(simplifiedPointArray, Length(originalPointArray));
    simplifiedPointArrayLength:=PolySimplifyFloat3D(apprimationAccuracy,originalPointArray,simplifiedPointArray);
    SetLength(simplifiedPointArray,simplifiedPointArrayLength);
    // convert a list of point to a list of segment
    SetLength(Result,simplifiedPointArrayLength-1);
    s:=0;
    for p:=0 to simplifiedPointArrayLength-2 do begin
      Result[s].p1:=simplifiedPointArray[p];
      Result[s].p2:=simplifiedPointArray[p+1];
      Inc(s);
    end;
    SetLength(originalPointArray,0);
    SetLength(simplifiedPointArray,0);
  end;
end;

{******************************************************************************}

end.