{*********************************************}
{  TSmoothingFunction                         }
{  Copyright (c) 2002-2020 by Steema Software }
{  With permission from M. v. Engeland        }
{*********************************************}
unit TeeSpline;
{$I TeeDefs.inc}

{------------------------------------------------------------------------------}
{                                                                              }
{ This code was written by : M. v. Engeland                                    }
{                                                                              }
{ This code is copyright 2000 by                                               }
{       M. v. Engeland                                                         }
{                                                                              }
{                                                                              }
{   This is the followup of the B-splines component posted in 1997.            }
{   Some bugs were fixed, and some changes were made. The splines now use      }
{   dynamic allocation of memory when needed when adding points to the spline. }
{   If the number of points to be added is known in advance, the the precise   }
{   amount of memory can be allocated setting the capacity property.           }
{                                                                              }
{   The conditions for using this component however are not changed:           }
{   You may use this component in any way you want to, whether it is for       }
{   commercial purposes or not, as long as you don't hold me responsible for   }
{   any disasters it may cause. But I would appreciate it if you would send me }
{   an e-mail (martijn@dutw38.wbmt.tudelft.nl) to tell me what you use it for  }
{   because over the last three years I've learned that it has been used for   }
{   a wide variaty of interesting applications which it was initially never    }
{   intended for. Besides, it's nice to know how your offspring's doing.       }
{   Also all comments and/or remarks are welcome.                              }
{                                                                              }
{   See the demo program for how to use the component.                         }
{   Special attention is payed on the possibility of interpolating the         }
{   vertices. As you may or may not know, B-Splines normally do not            }
{   interpolate the used controlpoints. Thanks to simple matrix calculation    }
{   however it is possible to interpolate the controlpoints by calculating     }
{   new vertices in such a way that the spline interpolates the original ones. }
{------------------------------------------------------------------------------}

interface

uses {$IFNDEF LINUX}
     {$IFNDEF FMX}
     Windows,
     {$ENDIF}
     {$ENDIF}
     Classes,

     {$IFDEF D18}
     System.Generics.Collections,
     {$ENDIF}
     
     {$IFDEF FMX}
     FMXTee.Canvas, FMXTee.Engine
     {$ELSE}
     TeCanvas, TeEngine
     {$ENDIF}
     ;

const
  MaxFragments            = 600; // The maximum of straight line segments allowed for drawing the spline
  MaxResults              = MaxFragments+10; // Max. number of calculated intersections

  // Removed max limit. No longer necessary.
  // Warning: More than 500 points takes huge CPU time.
  // MaxInterpolatedVertices :Integer= 250; // The maximum number of vertices that can be interpolated, up to 16000 allowed

  MaxCalcSteps            = 150; // Number of steps for numerical intersection calculating
  MaxError                = 1e-5;// Max error for intersection calculating
  MaxIterations           = 80;
  VerticesIncrement       = 5;  // Number of vertices to allocate memory for when the count property exceeds the current capacity

// BCC cannot handle constant floating-point numbers in headers that are precompiled.
// So 'MaxError' above makes TeeSpline.hpp (and any headers that use it) fail precompiled headers.
// To remedy we will avoid generating as a constant - and opt for a macro declaration instead
{$EXTERNALSYM MaxError}
{$HPPEMIT '#define TeeSpline_MaxError  ((System::Extended)1.000000E-05)'}

type
  TDataType=Double;

  TVertex=record
    X,Y : TDataType;
  end;

  // The following dynamic array is used to store the desired user-specified controlpoints
  T2DPointList = array of TVertex;

  // The vertexlist is used internally to make the spline interpolate the controlpoints
  TVertexList  = array of TVertex;

  // The knuckle list stores a flag to see whether a point is a knuckle or not
  TKnuckleList = array of Boolean;

  // The following tpes are used for the interpolation routines
  TMatrixRow = array of TDataType;

   // 2D B-spline class
  TBSpline=class
  private
    Matrix         : Array of TMatrixRow;
    FNoPoints      : Integer;
    FCapacity      : Integer;
    FPointList     : T2DPointList;
    FVertexList    : TVertexList;
    FKnuckleList   : TKnuckleList;
    FBuild         : Boolean;
    FNoVertices    : Integer;
    FInterpolated  : Boolean;
    FFragments     : Integer;

    procedure ClearVertexList;
    function CreateNewVertexList:TVertexList;
    function GetKnuckle(Index:Integer):Boolean;
    function GetNumberOfVertices:Integer;
    function GetPoint(Index:Integer):TVertex;
    procedure Interpolate;
    procedure PhantomPoints;
    procedure SetBuild(Val:Boolean);
    procedure SetCapacity(NewCapacity:Integer);
    procedure SetFragments(Const Value:Integer);
    procedure SetInterpolated(const Value:Boolean);
    procedure SetKnuckle(Index:Integer; Value:Boolean);
    procedure SetPoint(Index:Integer; const Value:TVertex);
   public
    Constructor Create; overload;
    Constructor Create(const P:Array of TPointFloat); overload;

    Destructor Destroy; override;

    procedure AddPoint(const X,Y:TDataType);
    procedure Clear;
    procedure Rebuild;
    function Value(const Parameter:TDataType):TVertex;

    property Build: Boolean read FBuild write SetBuild;
    property Count:Integer read FNoPoints;
    property Fragments: Integer read FFragments write SetFragments;
    property Knuckle[Index:integer]: Boolean read GetKnuckle write SetKnuckle;
    property Interpolated: Boolean read FInterpolated write SetInterpolated;
    property NumberOfVertices: Integer read GetNumberOfVertices;
    property Point[Index:Integer]: TVertex read GetPoint write SetPoint;
  end;

  TSmoothingFunction=class(TTeeFunction)
  private
    FInterpolate: Boolean;
    FFactor: Integer;

    procedure SetFactor(const Value: Integer);
    procedure SetInterpolate(const Value: Boolean);
  protected
    class function GetEditorClass: String; override;
    procedure PrepareForGallery(const Chart:TCustomAxisPanel); override;
  public
    Constructor Create(AOwner: TComponent); override;
    procedure AddPoints(Source:TChartSeries); override;
  published
    property Interpolate: Boolean read FInterpolate write SetInterpolate default True;
    property Factor:Integer read FFactor write SetFactor default 4;
  end;

implementation

uses
  Math, SysUtils,

  {$IFDEF FMX}
  FMXTee.Chart, FMXTee.Constants
  {$ELSE}
  Chart, TeeConst
  {$ENDIF}
  ;

{ TBSpline }

Constructor TBSpline.Create;
begin
  inherited;
  Clear;
end;

Constructor TBSpline.Create(const P: array of TPointFloat);
var t : Integer;
begin
  Create;

  for t:=Low(P) to High(P) do
  begin
    AddPoint(P[t].X,P[t].Y);
    Knuckle[t]:=False;
  end;
end;

Destructor TBSpline.Destroy;
begin
  Clear;
  inherited;
end;

procedure TBSpline.SetBuild(Val:Boolean);
begin
  if not Val then
  begin
    // Release allocated memory for vertices
    if FBuild then
       ClearVertexList;

    FNoVertices:=0;
  end;

  FBuild:=Val;
end;

procedure TBSpline.SetCapacity(NewCapacity:Integer);
var t : Integer;
begin
  if NewCapacity<>FCapacity then
  begin
    if NewCapacity>0 then
    begin
      SetLength(FPointList,NewCapacity+1 {*SizeOf(TVertex)} );
      SetLength(FKnuckleList,NewCapacity);

      for t:=0 to NewCapacity-1 do FKnuckleList[t]:=False;
    end;

    FCapacity:=NewCapacity;
  end;
end;

procedure TBSpline.SetFragments(Const Value:integer);
begin
  if FFragments<>value then
  begin
    FFragments:=Value;

    if FFragments>MaxFragments then
       FFragments:=MaxFragments;
 end;
end;

procedure TBSpline.SetInterpolated(const Value:boolean);
begin
  if Value<>FInterpolated then
  begin
    FInterpolated:=Value;
    Build:=False;
  end;
end;

function TBSpline.GetPoint(Index:Integer):TVertex;
begin
  Result:=FPointList[Index];
end;

procedure TBSpline.SetPoint(Index:Integer; const Value:TVertex);
begin
  FPointList[Index]:=Value;
  Build:=False;
end;

function TBSpline.GetKnuckle(Index:integer):Boolean;
begin
  if (Index=1) or (Index=FNoPoints) then
     result:=False
  else
     result:=FKnuckleList[Index]
end;

procedure TBSpline.SetKnuckle(Index:Integer; Value:Boolean);
begin
  FKnuckleList[Index]:=Value;
  Build:=False;
end;

function TBSpline.GetNumberOfVertices:integer;
begin
  if not FBuild then
     Rebuild;

  result:=FNoVertices;
end;

procedure TBSpline.Rebuild;

  procedure FillMatrix;
  Const MinLimit=1e-5;
  var I,J : Integer;
  begin
    if FNoVertices>2 then
    begin
      for i:=2 to FNoVertices-1 do
      begin
        Matrix[I][I-1]:=1/6;
        Matrix[I][I]:=2/3;
        Matrix[I][I+1]:=1/6;
      end;

      Matrix[1][1]:=1;
      Matrix[FNoVertices][FNoVertices]:=1;

      I:=3;
      while I<FNoVertices-1 do
      begin
        if (Abs(FVertexList[I].X-FVertexList[I-1].X)<MinLimit) and
           (Abs(FVertexList[I+1].X-FVertexList[I].X)<MinLimit) and
           (Abs(FVertexList[I].Y-FVertexList[I-1].Y)<MinLimit) and
           (Abs(FVertexList[I+1].Y-FVertexList[I].Y)<MinLimit) then
        begin
          for J:=I-1 to I+1 do
          begin
            Matrix[J][J-1]:=0;
            Matrix[J][J]:=1;
            Matrix[J][J+1]:=0;
          end;

          Inc(I,2);
        end
        else Inc(I);
      end;
    end;
  end;

var I,J,t    : integer;
    Vertex2D : TVertex;
begin
  if FNoPoints>1 then
  begin
    ClearVertexList;
    FNoVertices:=0;

    for i:=1 to FNoPoints do
        if Knuckle[I] then Inc(FNoVertices,3)
                      else Inc(FNoVertices,1);

    SetLength(FVertexList,FNoVertices+2);

    J:=0;
    for i:=1 to FNoPoints do
    begin
      Vertex2D:=FPointList[I];

      if Knuckle[I] then
      begin
        FVertexList[J+1]:=Vertex2D;
        FVertexList[J+2]:=Vertex2D;
        Inc(J,2);
      end;

      FVertexList[J+1]:=FPointList[I];
      Inc(J);
    end;

    if Interpolated then
    begin
      // Init Matrix
      SetLength(Matrix,FNoVertices+1);

      for i:=1 to FNoVertices do
      begin
        SetLength(Matrix[I],FNoVertices+1);

        for t:=0 to FNoVertices-1 do
            Matrix[i][t]:=0;
      end;

      FillMatrix;
      Interpolate;

      // Release memory
      for i:=1 to FNoVertices do Matrix[I]:=nil;

      Matrix:=nil;
    end;
  end;

  FBuild:=true;
  PhantomPoints;
end;

function TBSpline.CreateNewVertexList:TVertexList;
var tmpI,
    tmpINext   : TDataType;
    tmpMatrixI : TMatrixRow;

    I : Integer;
begin
  SetLength(result,FNoVertices+2);

  tmpI:=1/Matrix[FNoVertices][FNoVertices];

  result[FNoVertices].x:=FVertexList[FNoVertices].x*tmpI;
  result[FNoVertices].y:=FVertexList[FNoVertices].y*tmpI;

  for I:=FNoVertices-1 downto 1 do
  begin
    tmpMatrixI:=Matrix[I];

    tmpI:=1/tmpMatrixI[I];

    tmpINext:=tmpMatrixI[I+1];

    result[I].x:=tmpI*(FVertexList[I].x-tmpINext*result[I+1].x);
    result[I].y:=tmpI*(FVertexList[I].y-tmpINext*result[I+1].y);
  end;
end;

procedure TBSpline.Interpolate;
var I,J,K  : Integer;
    tmpI : TDataType;
    tmpMatrixI,
    tmpJ   : TMatrixRow;
    Factor : TDataType;
begin
  if FNoVertices>2 then
  begin
    for i:=1 to FNoVertices do
    begin
      tmpMatrixI:=Matrix[I];

      tmpI:=1/tmpMatrixI[I];

      for J:=I+1 to FNoVertices do
      begin
        tmpJ:=Matrix[J];

        factor:=tmpJ[I]*tmpI;

        for K:=1 to FNoVertices do
            tmpJ[K]:=tmpJ[K]-(factor*tmpMatrixI[K]);

        FVertexList[J].x:=FVertexList[J].x-factor*FVertexList[J-1].x;
        FVertexList[J].y:=FVertexList[J].y-factor*FVertexList[J-1].y;
      end;
    end;

    FVertexList:=CreateNewVertexList;
  end;
end;

procedure TBSpline.AddPoint(const X,Y:TDataType);
begin
  if FNoPoints=FCapacity then
     SetCapacity(FCapacity+VerticesIncrement);

  Inc(FNoPoints);
  FPointList[FNoPoints].X:=X;
  FPointList[FNoPoints].Y:=Y;
  Build:=False;
end;

procedure TBSpline.ClearVertexList;
begin
  FVertexList:=nil;
end;

procedure TBSpline.Clear;
begin
  if NumberOfVertices>0 then ClearVertexList;
  FNoPoints:=0;
  FNoVertices:=0;
  Build:=False;
  SetCapacity(0);
  FInterpolated:=False;
  FFragments:=100;
end;

procedure TBSpline.PhantomPoints;
var I : integer;
begin
  if NumberOfVertices>1 then
  begin
    I:=0;
    FVertexList[I].X:=2*FVertexList[I+1].X-FVertexList[I+2].X;
    FVertexList[I].Y:=2*FVertexList[I+1].Y-FVertexList[I+2].Y;

    FVertexList[NumberOfVertices+1].X:=2*FVertexList[NumberOfVertices].X-FVertexList[NumberOfVertices-1].X;
    FVertexList[NumberOfVertices+1].Y:=2*FVertexList[NumberOfVertices].Y-FVertexList[NumberOfVertices-1].Y;
  end;
end;

function TBSpline.Value(const Parameter:TDataType):TVertex;
const FourSixths=4/6;
      OneSixth=1/6;
var c,S,E : Integer;
    Dist  : TDataType;
    Mix   : TDataType;
    Mid   : TDataType;
    tmp   : TDataType;
begin
  result.X:=0;
  result.Y:=0;

  if FNoPoints<2 then Exit;

  if not FBuild then Rebuild;

  Mid:=(FNoVertices-1)*Parameter+1;
  S:=Trunc(Mid-1);
  if S<0 then S:=0;
  E:=S+3;
  if S>FNoVertices+1 then S:=FNoVertices+1;

  for c:=S to E do
  begin
    dist:=Abs(c-Mid);

    if dist<2 then
    begin
      if dist<1 then
         mix:=FourSixths-dist*dist+0.5*dist*dist*dist
      else
      begin
        tmp:=2-dist;
        mix:=tmp*tmp*tmp*OneSixth;
      end;

      result.x:=result.x+FVertexList[c].x*mix;
      result.y:=result.y+FVertexList[c].y*mix;
    end;
  end;
end;

{ TSmoothingFunction }
Constructor TSmoothingFunction.Create(AOwner: TComponent);
begin
  inherited;
  CanUsePeriod:=False;
  SingleSource:=True;
  FInterpolate:=True;
  FFactor:=4;
  InternalSetPeriod(1);
end;

procedure TSmoothingFunction.AddPoints(Source: TChartSeries);
var BSpline : TBSpline;
    t       : Integer;
    tmpList : TChartValueList;
    tmp     : TDataType;
begin
  with ParentSeries do
  begin
    Clear;

    // Improvement:  "if any of the ValuesLists has Order"...
    if not HasZValues then
    begin
      if YMandatory=Source.YMandatory then
      begin
        NotMandatoryValueList.Order:=loAscending;
        MandatoryValueList.Order:=loNone;
        CalcVisiblePoints:=True;
      end
      else
      begin
        NotMandatoryValueList.Order:=loNone;
        MandatoryValueList.Order:=loAscending;
        CalcVisiblePoints:=False;
      end;
    end;
  end;

  if Source.Count>0 then
  begin
    tmpList:=ValueList(Source);

    BSpline:=TBSpline.Create;
    try
      With Source do
      for t:=0 to Count-1 do
      begin
        if ParentSeries.YMandatory then
           BSpline.AddPoint(XValues.Value[t],tmpList.Value[t])
        else
           BSpline.AddPoint(YValues.Value[t],tmpList.Value[t]);

        BSpline.Knuckle[t]:=False;
      end;

      BSpline.Interpolated:=Interpolate;
      BSpline.Fragments:=Source.Count*Factor;

      with BSpline do
      begin
        tmp:=1/Fragments;

        for t:=0 to Fragments do
        with Value(t*tmp) do
             if ParentSeries.YMandatory then ParentSeries.AddXY(X,Y)
                                        else ParentSeries.AddXY(Y,X);
      end;
    finally
      BSpline.Free;
    end;
  end;
end;

procedure TSmoothingFunction.SetFactor(const Value: Integer);
begin
  if FFactor<>Value then
  begin
    FFactor:=Math.Max(1,Value);
    Recalculate;
  end;
end;

procedure TSmoothingFunction.SetInterpolate(const Value: Boolean);
begin
  if FInterpolate<>Value then
  begin
    FInterpolate:=Value;
    Recalculate;
  end;
end;

class function TSmoothingFunction.GetEditorClass: String;
begin
  result:='TSmoothFuncEditor';
end;

procedure TSmoothingFunction.PrepareForGallery(const Chart:TCustomAxisPanel);
begin
  inherited;

  Chart[0].FillSampleValues(10);
  ParentSeries.Transparency:=70;
  ParentSeries.Pen.Width:=10;
  Factor:=8;
end;

initialization
  RegisterTeeFunction( TSmoothingFunction, @TeeMsg_FunctionSmooth,
                                           @
     {$IFDEF TEELITE}TeeMsg_GalleryStandard{$ELSE}TeeMsg_GalleryExtended{$ENDIF});
finalization
  UnRegisterTeeFunctions([TSmoothingFunction]);
end.
