{*********************************************}
{  TDownSamplingFunction                      }
{  Copyright (c) 2004-2020 by Steema Software }
{  With permission from Mr. Nils Haeck        }
{*********************************************}
unit TeeDownSampling;
{$I TeeDefs.inc}

{ Implementation of the famous Douglas-Peucker 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
  {$IFNDEF LINUX}
  {$IFNDEF FMX}
  Windows,
  {$ENDIF}
  {$ENDIF}

  Classes, SysUtils,

  {$IFDEF FMX}
  FMX.Types,
  {$IFDEF D19}
  FMX.Graphics, // TBitmap
  {$ENDIF}
  {$ELSE}
  Graphics,
  {$ENDIF}
  {$IFDEF FMX}
  FMXTee.Engine, FMXTee.Canvas
  {$ELSE}
  TeEngine
  {$ENDIF};

type
  // Float point 3D
  TPointFloat2D = packed record
    X: TChartValue;
    Y: TChartValue;
  end;

  // Different style for downsampling calculation:

  TDownSampleMethod=(dsDouglas, dsMax, dsMin, dsMinMax, dsAverage,
                     dsMinMaxFirstLast, dsMinMaxFirstLastNull);

  TDownSamplingFunction=class(TTeeFunction)
  private
    FDisplayedPointCount: Integer;
    FDownSampleMethod: TDownSampleMethod;
    FTolerance: TChartValue;

    procedure SetDownSampleMethod(const Value: TDownSampleMethod);
    procedure SetTolerance(const Value:TChartValue);
    function Reduce(Method: TDownSampleMethod; const Tol: TChartValue;
                const lbound, ubound: Integer;
                const s: TChartSeries;
                var RX, RY: TChartValues;
                var colors: array of TColor;
                const color: TColor): Integer;
  protected
    class function GetEditorClass: String; override;
    procedure PrepareForGallery(const Chart:TCustomAxisPanel); override;
  public
    { Reduced size. Reduction factor is then Source.Count/ReducedSize}
    ReducedSize : Integer;

    Constructor Create(AOwner: TComponent); override;
    procedure AddPoints(Source:TChartSeries); override;
  published
    { Different downsample methods}
    property DownSampleMethod : TDownSampleMethod read FDownSampleMethod write SetDownSampleMethod default dsDouglas;

    { Tolerance, expressed in real chart values}
    property Tolerance:TChartValue read FTolerance write SetTolerance;

    { Force the actual number of points displayed }
    property DisplayedPointCount:Integer read FDisplayedPointCount write FDisplayedPointCount default 0;
  end;

implementation

uses Math,
     {$IFDEF FMX}
     FMXTee.Chart, FMXTee.ProConstants, FMXTee.Constants
     {$ELSE}

     {$IFDEF LCL}
     LCLType, Types,
     {$ENDIF}

     Chart, TeeConst, TeeProCo
     {$ENDIF}
     {$IFDEF D6}, DateUtils{$ENDIF};

function VecMin2D(const x1,y1,x2,y2: TChartValue): TPointFloat2D; {$IFDEF D9}inline;{$ENDIF}
// Result = (x1,y1) - (x2,y2)
begin
  Result.X := x1 - x2;
  Result.Y := y1 - y2;
end;

function DotProd2D(const x1,y1,x2,y2:TChartValue): TChartValue; {$IFDEF D9}inline;{$ENDIF}
// Dotproduct = (x1,y1)*(x2,y2)
begin
  Result := x1*x2 + y1*y2;
end;

function DistSquared2D(const x1,y1,x2,y2: TChartValue): TChartValue;
// Square of the distance from (x1,y1) to (x2,y2)
var v : TPointFloat2D;
begin
  v  := VecMin2D(x1,y1,x2,y2);
  result := Sqr(v.x)+Sqr(v.y);
end;

procedure Simplify2D(const Tol2: TChartValue; const X,Y: TChartValues;
  var Marker: array of boolean; j, k: integer);
// Simplify polyline in X,Y 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: TChartValue;    // Maximum value squared
  CU, CW, B: TChartValue;
  DV2: TChartValue;
  P0, P1, PB, U, W: TPointFloat2D;
begin
  // Is there anything to simplify?
  if k <= j + 1 then exit;

  P0.X := X[j];
  P0.Y := Y[j];

  P1.X := X[k];
  P1.Y := Y[k];
  U  := VecMin2D(P1.X,p1.Y, P0.X, P0.Y); // Segment vector
  CU := DotProd2D(U.X, U.Y, U.X, U.Y); // 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  := VecMin2D(X[i],Y[i], P0.X, P0.Y);
    CW := DotProd2D(W.X,W.Y, U.X, U.Y);

    // Distance of point Orig[i] from segment
    if CW <= 0 then
    begin
      // Before segment
      DV2 := DistSquared2D(X[i], Y[i], P0.X, P0.X)
    end
    else
    begin
      if CW > CU then
      begin
        // Past segment
        DV2 := DistSquared2D(X[i], Y[i], P1.X, P1.Y);
      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 := DistSquared2D(X[i], Y[i], PB.X, PB.Y);
      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]
    Simplify2D(Tol2, X,Y, Marker, j, MaxI); // polyline Orig[j] to Orig[maxi]
    Simplify2D(Tol2, X,Y, Marker, MaxI, k); // polyline Orig[maxi] to Orig[k]
  end;
end;

{ PolySimplify2D:
  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
          X,Y   = polyline array of vertex points
  Output: rX,rY = simplified polyline vertices. This array must initially
                     have the same length as X,Y
  Return: the number of points in rX or rY
}
function PolySimplify2D(const Tol: TChartValue; const X,Y: TChartValues; N: Integer;
                        out RX,RY: Array of TChartValue): Integer;
var Marker : Array of Boolean;
    i : Integer;
begin
  result := 0;

  if Length(X)>=2 then
  begin
    // Create a marker array initialized with False values
    SetLength(Marker, N);

    // Include first and last point
    Marker[0]     := True;
    Marker[N - 1] := True;

    // Simplify
    Simplify2D(Sqr(Tol), X,Y, Marker, 0, N - 1);

    // Copy to resulting list
    for i := 0 to N - 1 do
    if Marker[i] then
    begin
      RX[result] := X[i];
      RY[result] := Y[i];
      Inc(result);
    end;
  end;
end;

{ TDownSamplingFunction }

Constructor TDownSamplingFunction.Create(AOwner: TComponent);
begin
  inherited;

  CanUsePeriod:=False;
  SingleSource:=True;
  FTolerance := 1.0;
  InternalSetPeriod(1);
end;

class function TDownSamplingFunction.GetEditorClass: String;
begin
  result:='TDownSamplingEditor'; // Do not localize
end;

Procedure TDownSamplingFunction.PrepareForGallery(const Chart:TCustomAxisPanel);
begin
  inherited;

  ParentSeries.VertAxis:=aRightAxis;
  ParentSeries.GetVertAxis.StartPosition:=50;
  Chart[0].GetVertAxis.EndPosition:=50;

  Chart[0].FillSampleValues(2000);

  Chart[0].GetVertAxis.Title.Caption:=IntToStr(Chart[0].Count)+' '+TeeMsg_Points;
  ParentSeries.GetVertAxis.Title.Caption:=IntToStr(ParentSeries.Count)+' '+TeeMsg_Points;
end;

function TDownSamplingFunction.Reduce(Method: TDownSampleMethod; const Tol: TChartValue;
                const lbound, ubound: Integer;
                const s: TChartSeries;
                var RX, RY: TChartValues;
                var colors: array of TColor;
                const color: TColor): Integer;

  procedure InitColors;
  var k : Integer;
  begin
    for k:=0 to High(colors)-1 do
        colors[k]:=color;
  end;

var
  rSize : Integer;

  procedure AddPoint(const AX,AY:TChartValue);
  begin
    rX[rSize]:=AX;
    rY[rSize]:=AY;

    if AY=s.DefaultNullValue then
       colors[rSize]:=clNone;

    inc(rSize);
  end;

var
  tmpMax,
  tmpMin,
  tmpFirst,
  tmpLast : TChartValue;

  tmpXfirst,
  tmpXmax,
  tmpXmin,
  tmpXlast : TChartValue;

  procedure AddXRangePoints;
  begin
    if tmpXMax<=tmpXMin then
    begin
      rX[rSize+1]:=tmpXMax;
      rX[rSize+2]:=tmpXMin;
      rY[rSize+1]:=tmpMax;
      rY[rSize+2]:=tmpMin;
    end
    else
    begin
      rX[rSize+1]:=tmpXmin;
      rX[rSize+2]:=tmpXmax;
      rY[rSize+1]:=tmpMin;
      rY[rSize+2]:=tmpMax;
    end;

    rX[rSize+3]:=tmpXLast;
    rY[rSize+3]:=tmpLast;

    Inc(rSize,4);
  end;

  procedure AddDifferentXPoints;
  begin
    if tmpXMax<tmpXMin then
    begin
      if tmpXFirst<>tmpXMax then
         AddPoint(tmpXMax,tmpMax);

      if (tmpXFirst<>tmpXMin) and (tmpXMin<tmpXMax) then
         AddPoint(tmpXMin,tmpMin);
    end
    else
    begin
      if tmpXFirst<>tmpXMin then
         AddPoint(tmpXMin,tmpMin);

      if (tmpXFirst<>tmpXMax) and (tmpXMin<>tmpXMax) then
         AddPoint(tmpXMax,tmpMax);
    end;

    if (tmpXLast<>tmpXFirst) and (tmpXLast<>tmpXMin) and (tmpXLast<>tmpXMax) then
        AddPoint(tmpXLast,tmpLast);
  end;

  procedure AddThreePoints;
  var one, two, three,
      oneX, twoX, threeX: Double;
  begin
    if (tmpFirst=tmpMin) or (tmpLast=tmpMin) then
    begin
      one:=tmpFirst;
      oneX:=tmpXFirst;
      two:=tmpMax;
      twoX:=tmpXMax;
      three:=tmpLast;
      threeX:=tmpXLast;
    end
    else
    if (tmpFirst=tmpMax) or (tmpLast=tmpMax) then
    begin
      one:=tmpFirst;
      oneX:=tmpXFirst;
      two:=tmpMin;
      twoX:=tmpXMin;
      three:=tmpLast;
      threeX:=tmpXLast;
    end
    else
    begin
      one:=0;
      two:=0;
      three:=0;
      oneX:=0;
      twoX:=0;
      threeX:=0;
    end;

    rX[rSize]:=oneX;
    rX[rSize+1]:=twoX;
    rX[rSize+2]:=threeX;

    rY[rSize]:=one;
    rY[rSize+1]:=two;
    rY[rSize+2]:=three;

    Inc(rSize,3);
  end;

  procedure CheckMinMax(const X,Y:TChartValue);
  begin
    if Y>tmpMax then
    begin
      tmpMax:=Y;
      tmpXMax:=X;
    end;

    if Y<tmpMin then
    begin
      tmpMin:=Y;
      tmpXMin:=X;
    end;

    tmpLast:=Y;
    tmpXLast:=X;
  end;

var tmpSum : TChartValue;

    i,j,
    nullCount,iTol,count: Integer;

    X,Y: TChartValues;
    Xj,Yj : TChartValue;

    IsMinMaxAverage : Boolean;
begin
  X:=s.XValues.Value;
  Y:=s.YValues.Value;

  InitColors;

  if Method=dsDouglas then
     result:=PolySimplify2D(Tol,X,Y,ubound-lbound,rX,rY)
  else
  begin
    IsMinMaxAverage:=(Method=dsMin) or (Method=dsMax) or (Method=dsAverage);

    rSize:=0;
    i:=lbound;
    nullCount:=0;
    iTol:=Round(Tol);

    while i<ubound do
    begin
      j:=i;

      if Method=dsMinMaxFirstLastNull then
         if s.IsNull(i) then
            nullCount:=1
         else
            nullCount:=0;

      tmpMax:=Y[i];
      tmpMin:=tmpMax;
      tmpFirst:=tmpMax;
      tmpLast:=tmpMax;
      tmpSum:=tmpMax;

      tmpXMax:=X[i];
      tmpXFirst:=tmpXMax;
      tmpXMin:=tmpXMax;
      tmpXLast:=tmpXMax;

      // Calculate ranges

      if DisplayedPointCount>0 then
      begin
        count:=1;

        while ((j+1)<ubound) and (count<iTol) do
        begin
          inc(j);

          Yj:=Y[j];
          Xj:=X[j];

          if Method=dsMinMaxFirstLastNull then
          begin
            if s.IsNull(j) then inc(nullCount);

            if nullCount<=1 then
            begin
              if not s.IsNull(j) then
                 if tmpSum <> s.DefaultNullValue then
                    tmpSum:=tmpSum + Yj
                 else
                    tmpSum:=Yj;

              CheckMinMax(Xj,Yj);
            end;
          end
          else
          begin
            tmpSum:=tmpSum+Yj;

            CheckMinMax(Xj,Yj);
          end;

          inc(count);
        end;
      end
      else
        while ((j+1)<ubound) and (Abs(X[j+1]-X[i])<Tol) do
        begin
          Inc(j);

          Yj:=Y[j];

          tmpSum:=tmpSum+Yj;

          CheckMinMax(X[j],Yj);
        end;

      // Add points

      if IsMinMaxAverage then
      begin
        rX[rSize]:=(X[j] + X[i])*0.5; // x is average of group

        if Method=dsAverage then
           rY[rSize]:=tmpSum/(j-i+1)
        else
        if Method=dsMax then
           rY[rSize]:=tmpMax
        else
        if Method=dsMin then
           rY[rSize]:=tmpMin;

        inc(rsize);
      end
      else
      if rSize<=High(rX) then // safeguard in case somebody tries to do this on very short array
      begin
        if (j-i)=0 then  //CDI only drawing half the plot when tolerance = 0
        begin
          rX[rSize]:=X[i];
          rY[rSize]:=tmpMin;

          if Method=dsMinMaxFirstLastNull then
             if tmpMin=s.DefaultNullValue then
                colors[rSize]:=clNone;

          inc(rSize);
        end
        else
          if Method=dsMinMax then
          begin
            rX[rSize]:= X[i];
            rX[rSize+1]:= X[j];
            rY[rSize]:= tmpMin;
            rY[rSize+1]:= tmpMax;

            Inc(rSize,2);
          end
          else
            if (j-1)>2 then
            begin
              rX[rSize]:=tmpXFirst;
              rY[rSize]:=tmpFirst;

              if Method=dsMinMaxFirstLastNull then
              begin
                if tmpFirst=s.DefaultNullValue then
                   colors[rSize]:=clNone;

                inc(rSize);

                AddDifferentXPoints;
              end
              else
                AddXRangePoints;
            end
            else
            if (j-1)<2 then
            begin
              rX[rSize]:=X[i];
              rX[rSize+1]:=X[j];
              rX[rSize]:=Y[i];
              rX[rSize+1]:=Y[i];

              Inc(rSize,2);
            end
            else
              AddThreePoints;
          end;

        i:=Succ(j);
      end;

    result:=rSize;
  end;
end;

procedure TDownSamplingFunction.AddPoints(Source: TChartSeries);

  function FixedEmptyChartRect(const AChart:TCustomAxisPanel;
                               out N,lBound:Integer):Boolean;
  var tmpBitmap : TBitmap;
  begin
    result:=Assigned(AChart) and IsRectEmpty(AChart.ChartRect)
             // Avoid Parent Window exception when drawing onto bitmap
             and Assigned(AChart.Parent);

    if result then
    begin
      tmpBitmap:=AChart.TeeCreateBitmap;
      try
        if Source.ParentChart<>nil then
           Source.CalcFirstLastVisibleIndex;

        N:=Source.LastValueIndex - Source.FirstValueIndex;
        lbound:=Source.FirstValueIndex;
      finally
        tmpBitmap.Free;
      end;
    end;
  end;

  procedure ChangeOrders(const IsYMandatory:Boolean);
  begin
    if IsYMandatory then
    begin
      ParentSeries.NotMandatoryValueList.Order:=loAscending;
      ParentSeries.MandatoryValueList.Order:=loNone;
    end
    else
    begin
      ParentSeries.NotMandatoryValueList.Order:=loNone;
      ParentSeries.MandatoryValueList.Order:=loAscending;
    end;
  end;

  function CalcTolerance(const N:Integer):Double;
  begin
    result:=FTolerance;

    if DisplayedPointCount>0 then
       result:=4*N/DisplayedPointCount;
  end;

var
  rX,rY : TChartValues;
  colors: Array of TColor;

  procedure DoAddPoints;
  var t : Integer;
  begin
    if Source.YMandatory then
       for t:=0 to ReducedSize-1 do
           ParentSeries.AddXY(rX[t],rY[t],'',colors[t])
    else
       for t:=0 to ReducedSize-1 do
           ParentSeries.AddXY(rY[t],rX[t],'',colors[t]);
  end;

var N     : Integer;
    lbound: Integer;
begin
  if (not IUpdating) and Assigned(Source) and (Source.Count>0) then
  begin
    ParentSeries.Clear;

    if not FixedEmptyChartRect(ParentSeries.ParentChart,N,lbound) then
    begin
      N:=Source.Count;
      lbound:=0;
    end;

    if N>0 then
    begin
      ChangeOrders(Source.YMandatory=ParentSeries.YMandatory);

      // step 1 : setup result arrays
      SetLength(rX,N);
      SetLength(rY,N);
      SetLength(colors,N);

      // step 2 : do downsample
      ReducedSize := Reduce(FDownSampleMethod,
                            CalcTolerance(N),
                            lbound,
                            lbound+N,
                            Source,
                            rX,rY,
                            colors,
                            ParentSeries.Color);

      // step 3 : populate target series with results
      DoAddPoints;
    end;
  end;
end;

procedure TDownSamplingFunction.SetDownSampleMethod(const Value: TDownSampleMethod);
begin
  if FDownSampleMethod<>Value then
  begin
    FDownSampleMethod := Value;
    Recalculate;
  end;
end;

procedure TDownSamplingFunction.SetTolerance(const Value:TChartValue);
begin
  if FTolerance<>Value then
  begin
    if Value<0 then
       FTolerance := 0
    else
       FTolerance := Value;

    Recalculate;
  end;
end;

initialization
  RegisterTeeFunction( TDownSamplingFunction, @TeeMsg_FunctionDownSampling,
                                              @TeeMsg_GalleryExtended );
finalization
  UnRegisterTeeFunctions([TDownSamplingFunction]);
end.
