castle-game-engine-src 5.2.0-2 / usr / src / castle-game-engine-5.2.0 / opengl / castlecurves.pas

{
  Copyright 2004-2014 Michalis Kamburelis.

  This file is part of "Castle Game Engine".

  "Castle Game Engine" is free software; see the file COPYING.txt,
  included in this distribution, for details about the copyright.

  "Castle Game Engine" is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

  ----------------------------------------------------------------------------
}

{ 3D curves (TCurve and basic descendants). }
unit CastleCurves;

{$I castleconf.inc}

interface

uses Classes, FGL, CastleVectors, CastleBoxes, CastleUtils, CastleScript,
  CastleClassUtils, Castle3D, CastleFrustum, CastleColors;

type
  { 3D curve, a set of points defined by a continous function @link(Point)
    for arguments within [TBegin, TEnd].

    Note that some descendants return only an approximate BoundingBox result,
    it may be too small or too large sometimes.
    (Maybe at some time I'll make this more rigorous, as some code may require
    that it's a proper bounding box, maybe too large but never too small.) }
  TCurve = class(T3D)
  private
    FColor: TCastleColor;
    FLineWidth: Single;
    FTBegin, FTEnd: Single;
    FDefaultSegments: Cardinal;
  public
    { The valid range of curve function argument. Must be TBegin <= TEnd.
      @groupBegin }
    property TBegin: Single read FTBegin write FTBegin default 0;
    property TEnd: Single read FTEnd write FTEnd default 1;
    { @groupEnd }

    { Curve function, for each parameter value determine the 3D point.
      This determines the actual shape of the curve. }
    function Point(const t: Float): TVector3Single; virtual; abstract;

    { Curve function to work with rendered line segments begin/end points.
      This is simply a more specialized version of @link(Point),
      it scales the argument such that you get Point(TBegin) for I = 0
      and you get Point(TEnd) for I = Segments. }
    function PointOfSegment(i, Segments: Cardinal): TVector3Single;

    { Render curve by dividing it into a given number of line segments.
      So actually @italic(every) curve will be rendered as a set of straight lines.
      You should just give some large number for Segments to have something
      that will be really smooth.

      This does direct OpenGL rendering right now, setting GL color
      and then rendering a line strip. }
    procedure Render(Segments: Cardinal); deprecated 'Do not render curve directly by this method, instead add the curve to SceneManager.Items to have it rendered automatically.';

    { Curve rendering color. White by default. }
    property Color: TCastleColor read FColor write FColor;

    property LineWidth: Single read FLineWidth write FLineWidth default 1;

    { Default number of segments, used when rendering by T3D interface
      (that is, @code(Render(Frustum, TransparentGroup...)) method.) }
    property DefaultSegments: Cardinal
      read FDefaultSegments write FDefaultSegments default 10;

    procedure Render(const Frustum: TFrustum;
      const Params: TRenderParams); override;

    constructor Create(AOwner: TComponent); override;
  end;

  TCurveList = specialize TFPGObjectList<TCurve>;

  { Curve defined by explicitly giving functions for
    Point(t) = x(t), y(t), z(t) as CastleScript expressions. }
  TCasScriptCurve = class(TCurve)
  private
    FSegmentsForBoundingBox: Cardinal;
    procedure SetSegmentsForBoundingBox(AValue: Cardinal);
    procedure SetTVariable(AValue: TCasScriptFloat);
  protected
    FTVariable: TCasScriptFloat;
    FFunction: array [0..2] of TCasScriptExpression;
    FBoundingBox: TBox3D;
    function GetFunction(const Index: Integer): TCasScriptExpression;
    procedure SetFunction(const Index: Integer; const Value: TCasScriptExpression);
    procedure UpdateBoundingBox;
  public
    function Point(const t: Float): TVector3Single; override;

    { XFunction, YFunction, ZFunction are functions based on variable 't'.
      Once set, these instances become owned by this class, do not free
      them yourself!
      @groupBegin }
    property XFunction: TCasScriptExpression index 0 read GetFunction write SetFunction;
    property YFunction: TCasScriptExpression index 1 read GetFunction write SetFunction;
    property ZFunction: TCasScriptExpression index 2 read GetFunction write SetFunction;
    { @groupEnd }

    { This is the variable controlling 't' value, embedded also in
      XFunction, YFunction, ZFunction. This is NOT owned by this class,
      make sure to free it yourself! }
    property TVariable: TCasScriptFloat read FTVariable write SetTVariable;

    property SegmentsForBoundingBox: Cardinal
      read FSegmentsForBoundingBox write SetSegmentsForBoundingBox default 100;

    { Simple bounding box. It is simply
      a BoundingBox of Point(i, SegmentsForBoundingBox)
      for i in [0 .. SegmentsForBoundingBox].
      Subclasses may override this to calculate something more accurate. }
    function BoundingBox: TBox3D; override;

    constructor Create(AOwner: TComponent); override;

    destructor Destroy; override;
  end;

  { A basic abstract class for curves determined my some set of ControlPoints.
    Note: it is @italic(not) defined in this class any correspondence between
    values of T (argument for Point function) and ControlPoints. }
  TControlPointsCurve = class(TCurve)
  private
    FBoundingBox: TBox3D;
    FControlPointsColor: TCastleColor;
    FConvexHullColor: TCastleColor;
  protected
    { Using these function you can control how Convex Hull (for RenderConvexHull)
      is calculated: CreateConvexHullPoints should return points that must be
      in convex hull (we will run ConvexHull function on those points),
      DestroyConvexHullPoints should finalize them.

      This way you can create new object in CreateConvexHullPoints and free it in
      DestroyConvexHullPoints, but you can also give in CreateConvexHullPoints
      reference to some already existing object and do nothing in
      DestroyConvexHullPoints. (we will not modify object given as
      CreateConvexHullPoints in any way)

      Default implementation in this class returns ControlPoints as
      CreateConvexHullPoints. (and does nothing in DestroyConvexHullPoints) }
    function CreateConvexHullPoints: TVector3SingleList; virtual;
    procedure DestroyConvexHullPoints(Points: TVector3SingleList); virtual;
  public
    ControlPoints: TVector3SingleList;

    property ControlPointsColor: TCastleColor read FControlPointsColor write FControlPointsColor;

    { Render control points, using ControlPointsColor. }
    procedure RenderControlPoints;

    { This class provides implementation for BoundingBox: it is simply
      a BoundingBox of ControlPoints. Subclasses may (but don't have to)
      override this to calculate better (more accurate) BoundingBox. }
    function BoundingBox: TBox3D; override;

    { Always after changing ControlPoints or TBegin or TEnd and before calling Point,
      BoundingBox (and anything that calls them, e.g. Render calls Point)
      call this method. It recalculates necessary things.
      ControlPoints.Count must be >= 2.

      When overriding: always call inherited first. }
    procedure UpdateControlPoints; virtual;

    { Nice class name, with spaces, starts with a capital letter.
      Much better than ClassName. Especially under FPC 1.0.x where
      ClassName is always uppercased. }
    class function NiceClassName: string; virtual; abstract;

    property ConvexHullColor: TCastleColor read FConvexHullColor write FConvexHullColor;

    { Render convex hull polygon, using ConvexHullColor.
      Ignores Z-coord of ControlPoints. }
    procedure RenderConvexHull;

    { Constructor.
      It has to be virtual because it's called by CreateDivideCasScriptCurve. }
    constructor Create(AOwner: TComponent); override;

    { Calculates ControlPoints taking Point(i, ControlPointsCount-1)
      for i in [0 .. ControlPointsCount-1] from CasScriptCurve.
      TBegin and TEnd are copied from CasScriptCurve. }
    constructor CreateDivideCasScriptCurve(CasScriptCurve: TCasScriptCurve;
      ControlPointsCount: Cardinal);

    destructor Destroy; override;
  end;

  TControlPointsCurveClass = class of TControlPointsCurve;

  TControlPointsCurveList = specialize TFPGObjectList<TControlPointsCurve>;

  { Curve that passes exactly through it's ControlPoints.x
    I.e. for each ControlPoint[i] there exists some value Ti
    that Point(Ti) = ControlPoint[i] and
    TBegin = T0 <= .. Ti-1 <= Ti <= Ti+1 ... <= Tn = TEnd
    (i.e. Point(TBegin) = ControlPoints[0],
          Point(TEnd) = ControlPoints[n]
     and all Ti are ordered,
     n = ControlPoints.High) }
  TInterpolatedCurve = class(TControlPointsCurve)
    { This can be overriden in subclasses.
      In this class this provides the most common implementation:
      equally (uniformly) spaced Ti values. }
    function ControlPointT(i: Integer): Float; virtual;
  end;

  { Curve defined as [Lx(t), Ly(t), Lz(t)] where
    L?(t) are Lagrange's interpolation polynomials.
    Lx(t) crosses points (ti, xi) (i = 0..ControlPoints.Count-1)
    where ti = TBegin + i/(ControlPoints.Count-1) * (TEnd-TBegin)
    and xi = ControlPoints[i, 0].
    Similarly for Ly and Lz.

    Later note: in fact, you can override ControlPointT to define
    function "ti" as you like. }
  TLagrangeInterpolatedCurve = class(TInterpolatedCurve)
  private
    { Values for Newton's polynomial form of Lx, Ly and Lz.
      Will be calculated in UpdateControlPoints. }
    Newton: array[0..2]of TFloatList;
  public
    procedure UpdateControlPoints; override;
    function Point(const t: Float): TVector3Single; override;

    class function NiceClassName: string; override;

    constructor Create(AOwner: TComponent); override;
    destructor Destroy; override;
  end deprecated 'Rendering of this is not portable to OpenGLES, and this is not really a useful curve for most practical game uses. For portable and fast curves consider using X3D NURBS nodes (wrapped in a TCastleScene) instead.';

  { Natural cubic spline (1D).
    May be periodic or not. }
  TNaturalCubicSpline = class
  private
    FMinX, FMaxX: Float;
    FOwnsX, FOwnsY: boolean;
    FPeriodic: boolean;
    FX, FY: TFloatList;
    M: TFloatList;
  public
    property MinX: Float read FMinX;
    property MaxX: Float read FMaxX;
    property Periodic: boolean read FPeriodic;

    { Constructs natural cubic spline such that for every i in [0; X.Count-1]
      s(X[i]) = Y[i]. Must be X.Count = Y.Count.
      X must be already sorted.
      MinX = X[0], MaxX = X[X.Count-1].

      Warning: we will copy references to X and Y ! So make sure that these
      objects are available for the life of this object.
      We will free in destructor X if OwnsX and free Y if OwnsY. }
    constructor Create(X, Y: TFloatList; AOwnsX, AOwnsY, APeriodic: boolean);
    destructor Destroy; override;

    { Evaluate value of natural cubic spline at x.
      Must be MinX <= x <= MaxX. }
    function Evaluate(x: Float): Float;
  end;

  { 3D curve defined by three 1D natural cubic splines.
    Works just like TLagrangeInterpolatedCurve, only the interpolation
    is different now. }
  TNaturalCubicSplineCurve_Abstract = class(TInterpolatedCurve)
  protected
    { Is the curve closed. May depend on ControlPoints,
      it will be recalculated in UpdateControlPoints. }
    function Closed: boolean; virtual; abstract;
  private
    { Created/Freed in UpdateControlPoints, Freed in Destroy }
    Spline: array[0..2]of TNaturalCubicSpline;
  public
    procedure UpdateControlPoints; override;
    function Point(const t: Float): TVector3Single; override;

    constructor Create(AOwner: TComponent); override;
    destructor Destroy; override;
  end;

  { 3D curve defined by three 1D natural cubic splines, automatically
    closed if first and last points match. This is the most often suitable
    non-abstract implementation of TNaturalCubicSplineCurve_Abstract. }
  TNaturalCubicSplineCurve = class(TNaturalCubicSplineCurve_Abstract)
  protected
    function Closed: boolean; override;
  public
    class function NiceClassName: string; override;
  end deprecated 'Rendering of this is not portable to OpenGLES, and this is not really a useful curve for most practical game uses. For portable and fast curves consider using X3D NURBS nodes (wrapped in a TCastleScene) instead.';

  { 3D curve defined by three 1D natural cubic splines, always treated as closed. }
  TNaturalCubicSplineCurveAlwaysClosed = class(TNaturalCubicSplineCurve_Abstract)
  protected
    function Closed: boolean; override;
  public
    class function NiceClassName: string; override;
  end deprecated 'Rendering of this is not portable to OpenGLES, and this is not really a useful curve for most practical game uses. For portable and fast curves consider using X3D NURBS nodes (wrapped in a TCastleScene) instead.';

  { 3D curve defined by three 1D natural cubic splines, never treated as closed. }
  TNaturalCubicSplineCurveNeverClosed = class(TNaturalCubicSplineCurve_Abstract)
  protected
    function Closed: boolean; override;
  public
    class function NiceClassName: string; override;
  end deprecated 'Rendering of this is not portable to OpenGLES, and this is not really a useful curve for most practical game uses. For portable and fast curves consider using X3D NURBS nodes (wrapped in a TCastleScene) instead.';

  { Rational Bezier curve (Bezier curve with weights).
    Note: for Bezier Curve ControlPoints.Count MAY be 1.
    (For TControlPointsCurve it must be >= 2) }
  TRationalBezierCurve = class(TControlPointsCurve)
  public
    { Splits this curve using Casteljau algorithm.

      Under B1 and B2 returns two new, freshly created, bezier curves,
      such that if you concatenate them - they will create this curve.
      Proportion is something from (0; 1).
      B1 will be equal to Self for T in TBegin .. TMiddle,
      B2 will be equal to Self for T in TMiddle .. TEnd,
      where TMiddle = TBegin + Proportion * (TEnd - TBegin).

      B1.ControlPoints.Count = B2.ControlPoints.Count =
        Self.ControlPoints.Count. }
    procedure Split(const Proportion: Float; var B1, B2: TRationalBezierCurve);

    function Point(const t: Float): TVector3Single; override;
    class function NiceClassName: string; override;
  public
    { Curve weights.
      Must always be Weights.Count = ControlPoints.Count.
      After changing Weights you also have to call UpdateControlPoints.}
    Weights: TFloatList;

    procedure UpdateControlPoints; override;

    constructor Create(AOwner: TComponent); override;
    destructor Destroy; override;
  end deprecated 'Rendering of TRationalBezierCurve is not portable to OpenGLES (that is: Android and iOS) and not very efficient. For portable and fast curves consider using X3D NURBS nodes (wrapped in a TCastleScene) instead.';

  {$warnings off} { Consciously using deprecated stuff. }
  TRationalBezierCurveList = specialize TFPGObjectList<TRationalBezierCurve>;
  {$warnings on}

  { Smooth interpolated curve, each segment (ControlPoints[i]..ControlPoints[i+1])
    is converted to a rational Bezier curve (with 4 control points)
    when rendering.

    You can also explicitly convert it to a list of bezier curves using
    ToRationalBezierCurves.

    Here too ControlPoints.Count MAY be 1.
    (For TControlPointsCurve it must be >= 2). }
  TSmoothInterpolatedCurve = class(TInterpolatedCurve)
  private
    BezierCurves: TRationalBezierCurveList;
    ConvexHullPoints: TVector3SingleList;
  protected
    function CreateConvexHullPoints: TVector3SingleList; override;
    procedure DestroyConvexHullPoints(Points: TVector3SingleList); override;
  public
    function Point(const t: Float): TVector3Single; override;

    { convert this to a list of TRationalBezierCurve.

      From each line segment ControlPoint[i] ... ControlPoint[i+1]
      you get one TRationalBezierCurve with 4 control points,
      where ControlPoint[0] and ControlPoint[3] are taken from
      ours ControlPoint[i] ... ControlPoint[i+1] and the middle
      ControlPoint[1], ControlPoint[2] are calculated so that all those
      bezier curves join smoothly.

      All Weights are set to 1.0 (so actually these are all normal
      Bezier curves; but I'm treating normal Bezier curves as Rational
      Bezier curves everywhere here) }
    function ToRationalBezierCurves(ResultOwnsCurves: boolean): TRationalBezierCurveList;

    procedure UpdateControlPoints; override;

    class function NiceClassName: string; override;

    constructor Create(AOwner: TComponent); override;
    destructor Destroy; override;
  end deprecated 'Rendering of TSmoothInterpolatedCurve is not portable to OpenGLES (that is: Android and iOS) and not very efficient. For portable and fast curves consider using X3D NURBS nodes (wrapped in a TCastleScene) instead.';

implementation

uses SysUtils, CastleGL, CastleConvexHull, CastleGLUtils;

{ TCurve ------------------------------------------------------------ }

function TCurve.PointOfSegment(i, Segments: Cardinal): TVector3Single;
begin
  Result := Point(TBegin + (i/Segments) * (TEnd-TBegin));
end;

procedure TCurve.Render(Segments: Cardinal);
var i: Integer;
begin
  {$ifndef OpenGLES} //TODO-es
  glColorv(Color);
  glLineWidth(LineWidth);
  glBegin(GL_LINE_STRIP);
  for i := 0 to Segments do glVertexv(PointOfSegment(i, Segments));
  glEnd;
  {$endif}
end;

procedure TCurve.Render(const Frustum: TFrustum;
  const Params: TRenderParams);
begin
  if GetExists and (not Params.Transparent) and Params.ShadowVolumesReceivers then
  begin
    {$ifndef OpenGLES} //TODO-es
    if not Params.RenderTransformIdentity then
    begin
      glPushMatrix;
      glMultMatrix(Params.RenderTransform);
    end;
    {$endif}

    {$warnings off}
    Render(DefaultSegments);
    {$warnings on}

    {$ifndef OpenGLES}
    if not Params.RenderTransformIdentity then
      glPopMatrix;
    {$endif}
  end;
end;

constructor TCurve.Create(AOwner: TComponent);
begin
  inherited;
  FTBegin := 0;
  FTEnd := 1;
  FDefaultSegments := 10;
  FLineWidth := 1;
  FColor := White;
end;

{ TCasScriptCurve ------------------------------------------------------------ }

procedure TCasScriptCurve.SetTVariable(AValue: TCasScriptFloat);
begin
  if FTVariable = AValue then Exit;
  FTVariable := AValue;
  UpdateBoundingBox;
end;

procedure TCasScriptCurve.SetSegmentsForBoundingBox(AValue: Cardinal);
begin
  if FSegmentsForBoundingBox = AValue then Exit;
  FSegmentsForBoundingBox := AValue;
  UpdateBoundingBox;
end;

function TCasScriptCurve.GetFunction(const Index: Integer): TCasScriptExpression;
begin
  Result := FFunction[Index];
end;

procedure TCasScriptCurve.SetFunction(const Index: Integer;
  const Value: TCasScriptExpression);
begin
  if FFunction[Index] = Value then Exit;

  if FFunction[Index] <> nil then
    FFunction[Index].FreeByParentExpression;

  FFunction[Index] := Value;
  UpdateBoundingBox;
end;

procedure TCasScriptCurve.UpdateBoundingBox;
var
  i, k: Integer;
  P: TVector3Single;
begin
  if (XFunction = nil) or
     (YFunction = nil) or
     (ZFunction = nil) or
     (TVariable = nil) then
    FBoundingBox := EmptyBox3D else
  begin
    { calculate FBoundingBox }
    P := PointOfSegment(0, SegmentsForBoundingBox); { = Point(TBegin) }
    FBoundingBox.Data[0] := P;
    FBoundingBox.Data[1] := P;
    for i := 1 to SegmentsForBoundingBox do
    begin
      P := PointOfSegment(i, SegmentsForBoundingBox);
      for k := 0 to 2 do
      begin
        FBoundingBox.Data[0, k] := Min(FBoundingBox.Data[0, k], P[k]);
        FBoundingBox.Data[1, k] := Max(FBoundingBox.Data[1, k], P[k]);
      end;
    end;
  end;
end;

function TCasScriptCurve.Point(const t: Float): TVector3Single;
var
  I: Integer;
begin
  TVariable.Value := T;
  for I := 0 to 2 do
    Result[I] := (FFunction[I].Execute as TCasScriptFloat).Value;

  {test: Writeln('Point at t = ',FloatToNiceStr(Single(t)), ' is (',
    VectorToNiceStr(Result), ')');}
end;

function TCasScriptCurve.BoundingBox: TBox3D;
begin
  Result := FBoundingBox;
end;

constructor TCasScriptCurve.Create(AOwner: TComponent);
begin
  inherited;
  FSegmentsForBoundingBox := 100;
end;

destructor TCasScriptCurve.Destroy;
var
  I: Integer;
begin
  for I := 0 to 2 do
    if FFunction[I] <> nil then
    begin
      FFunction[I].FreeByParentExpression;
      FFunction[I] := nil;
    end;
  inherited;
end;

{ TControlPointsCurve ------------------------------------------------ }

procedure TControlPointsCurve.RenderControlPoints;
var
  i: Integer;
begin
  {$ifndef OpenGLES} //TODO-es
  glColorv(ControlPointsColor);
  glBegin(GL_POINTS);
  for i := 0 to ControlPoints.Count-1 do glVertexv(ControlPoints.L[i]);
  glEnd;
  {$endif}
end;

function TControlPointsCurve.BoundingBox: TBox3D;
begin
  Result := FBoundingBox;
end;

procedure TControlPointsCurve.UpdateControlPoints;
begin
  FBoundingBox := CalculateBoundingBox(PVector3Single(ControlPoints.List),
    ControlPoints.Count, 0);
end;

function TControlPointsCurve.CreateConvexHullPoints: TVector3SingleList;
begin
  Result := ControlPoints;
end;

procedure TControlPointsCurve.DestroyConvexHullPoints(Points: TVector3SingleList);
begin
end;

procedure TControlPointsCurve.RenderConvexHull;
var
  CHPoints: TVector3SingleList;
  CH: TIntegerList;
  i: Integer;
begin
  CHPoints := CreateConvexHullPoints;
  try
    CH := ConvexHull(CHPoints);
    try
      {$ifndef OpenGLES} //TODO-es
      glColorv(ConvexHullColor);
      glBegin(GL_POLYGON);
      try
        for i := 0 to CH.Count-1 do
          glVertexv(CHPoints.L[CH[i]]);
      finally glEnd end;
      {$endif}
    finally CH.Free end;
  finally DestroyConvexHullPoints(CHPoints) end;
end;

constructor TControlPointsCurve.Create(AOwner: TComponent);
begin
  inherited;
  ControlPoints := TVector3SingleList.Create;
  { DON'T call UpdateControlPoints from here - UpdateControlPoints is virtual !
    So we set FBoundingBox by hand. }
  FBoundingBox := EmptyBox3D;
  FControlPointsColor := White;
  FConvexHullColor := White;
end;

constructor TControlPointsCurve.CreateDivideCasScriptCurve(
  CasScriptCurve: TCasScriptCurve; ControlPointsCount: Cardinal);
var
  i: Integer;
begin
  Create(nil);
  TBegin := CasScriptCurve.TBegin;
  TEnd := CasScriptCurve.TEnd;
  ControlPoints.Count := ControlPointsCount;
  for i := 0 to ControlPointsCount-1 do
    ControlPoints.L[i] := CasScriptCurve.PointOfSegment(i, ControlPointsCount-1);
  UpdateControlPoints;
end;

destructor TControlPointsCurve.Destroy;
begin
  FreeAndNil(ControlPoints);
  inherited;
end;

{ TInterpolatedCurve ----------------------------------------------- }

function TInterpolatedCurve.ControlPointT(i: Integer): Float;
begin
  Result := TBegin + (i/(ControlPoints.Count-1)) * (TEnd-TBegin);
end;

{ TLagrangeInterpolatedCurve ----------------------------------------------- }

procedure TLagrangeInterpolatedCurve.UpdateControlPoints;
var
  i, j, k, l: Integer;
begin
  inherited;

  for i := 0 to 2 do
  begin
    Newton[i].Count := ControlPoints.Count;
    for j := 0 to ControlPoints.Count-1 do
      Newton[i].L[j] := ControlPoints.L[j, i];

    { licz kolumny tablicy ilorazow roznicowych in place, overriding Newton[i] }
    for k := 1 to ControlPoints.Count-1 do
      { licz k-ta kolumne }
      for l := ControlPoints.Count-1 downto k do
        { licz l-ty iloraz roznicowy w k-tej kolumnie }
        Newton[i].L[l]:=
          (Newton[i].L[l] - Newton[i].L[l-1]) /
          (ControlPointT(l) - ControlPointT(l-k));
  end;
end;

function TLagrangeInterpolatedCurve.Point(const t: Float): TVector3Single;
var
  i, k: Integer;
  f: Float;
begin
  for i := 0 to 2 do
  begin
    { Oblicz F przy pomocy uogolnionego schematu Hornera z Li(t).
      Wspolczynniki b_k sa w tablicy Newton[i].L[k],
      wartosci t_k sa w ControlPointT(k). }
    F := Newton[i].L[ControlPoints.Count-1];
    for k := ControlPoints.Count-2 downto 0 do
      F := F*(t-ControlPointT(k)) + Newton[i].L[k];
    { Dopiero teraz przepisz F do Result[i]. Dzieki temu obliczenia wykonujemy
      na Floatach. Tak, to naprawde pomaga -- widac ze kiedy uzywamy tego to
      musimy miec wiecej ControlPoints zeby dostac Floating Point Overflow. }
    Result[i] := F;
  end;
end;

class function TLagrangeInterpolatedCurve.NiceClassName: string;
begin
  Result := 'Lagrange interpolated curve';
end;

constructor TLagrangeInterpolatedCurve.Create(AOwner: TComponent);
var
  i: Integer;
begin
  inherited;
  for i := 0 to 2 do Newton[i] := TFloatList.Create;
end;

destructor TLagrangeInterpolatedCurve.Destroy;
var
  i: Integer;
begin
  for i := 0 to 2 do FreeAndNil(Newton[i]);
  inherited;
end;

{ TNaturalCubicSpline -------------------------------------------------------- }

constructor TNaturalCubicSpline.Create(X, Y: TFloatList;
  AOwnsX, AOwnsY, APeriodic: boolean);

{ Based on SLE (== Stanislaw Lewanowicz) notes on ii.uni.wroc.pl lecture. }

var
  { n = X.High. Integer, not Cardinal, to avoid some overflows in n-2. }
  n: Integer;

  { [Not]PeriodicDK licza wspolczynnik d_k. k jest na pewno w [1; n-1] }
  function PeriodicDK(k: Integer): Float;
  var
    h_k: Float;
    h_k1: Float;
  begin
    h_k  := X[k] - X[k-1];
    h_k1 := X[k+1] - X[k];
    Result := ( 6 / (h_k + h_k1) ) *
              ( (Y[k+1] - Y[k]) / h_k1 -
                (Y[k] - Y[k-1]) / h_k
              );
  end;

  { special version, works like PeriodicDK(n) should work
    (but does not, PeriodicDK works only for k < n.) }
  function PeriodicDN: Float;
  var
    h_n: Float;
    h_n1: Float;
  begin
    h_n  := X[n] - X[n-1];
    h_n1 := X[1] - X[0];
    Result := ( 6 / (h_n + h_n1) ) *
              ( (Y[1] - Y[n]) / h_n1 -
                (Y[n] - Y[n-1]) / h_n
              );
  end;

  function IlorazRoznicowy(const Start, Koniec: Cardinal): Float;
  { liczenie pojedynczego ilorazu roznicowego wzorem rekurencyjnym.
    Poniewaz do algorytmu bedziemy potrzebowali tylko ilorazow stopnia 3
    (lub mniej) i to tylko na chwile - wiec taka implementacja
    (zamiast zabawa w tablice) bedzie prostsza i wystarczajaca. }
  begin
    if Start = Koniec then
      Result := Y[Start] else
      Result := (IlorazRoznicowy(Start + 1, Koniec) -
                 IlorazRoznicowy(Start, Koniec - 1))
                 / (X[Koniec] - X[Start]);
  end;

  function NotPeriodicDK(k: Integer): Float;
  begin
    Result := 6 * IlorazRoznicowy(k-1, k+1);
  end;

var
  u, q, s, t, v: TFloatList;
  hk, dk, pk, delta_k, delta_n, h_n, h_n1: Float;
  k: Integer;
begin
  inherited Create;
  Assert(X.Count = Y.Count);
  FMinX := X.First;
  FMaxX := X.Last;
  FOwnsX := AOwnsX;
  FOwnsY := AOwnsY;
  FX := X;
  FY := Y;
  FPeriodic := APeriodic;

  { prepare to calculate M }
  n := X.Count - 1;
  M := TFloatList.Create;
  M.Count := n+1;

  { Algorytm obliczania wartosci M[0..n] z notatek SLE, te same oznaczenia.
    Sa tutaj polaczone algorytmy na Periodic i not Perdiodic, zeby mozliwie
    nie duplikowac kodu (i uniknac pomylek z copy&paste).
    Tracimy przez to troche czasu (wielokrotne testy "if Periodic ..."),
    ale kod jest prostszy i bardziej elegancki.

    Notka: W notatkach SLE dla Periodic = true w jednym miejscu uzyto
    M[n+1] ale tak naprawde nie musimy go liczyc ani uzywac. }

  u := nil;
  q := nil;
  s := nil;
  try
    u := TFloatList.Create; U.Count := N;
    q := TFloatList.Create; Q.Count := N;
    if Periodic then begin s := TFloatList.Create; S.Count := N; end;

    { calculate u[0], q[0], s[0] }
    u[0] := 0;
    q[0] := 0;
    if Periodic then s[0] := 1;

    for k := 1 to n - 1 do
    begin
      { calculate u[k], q[k], s[k] }

      hk := X[k] - X[k-1];
      { delta[k] = h[k] / (h[k] + h[k+1])
                 = h[k] / (X[k] - X[k-1] + X[k+1] - X[k])
                 = h[k] / (X[k+1] - X[k-1])
      }
      delta_k := hk / (X[k+1] - X[k-1]);
      pk := delta_k * q[k-1] + 2;
      q[k]:=(delta_k - 1) / pk;
      if Periodic then s[k] := - delta_k * s[k-1] / pk;
      if Periodic then
        dk := PeriodicDK(k) else
        dk := NotPeriodicDK(k);
      u[k]:=(dk - delta_k * u[k-1]) / pk;
    end;

    { teraz wyliczamy wartosci M[0..n] }
    if Periodic then
    begin
      t := nil;
      v := nil;
      try
        t := TFloatList.Create; T.Count := N + 1;
        v := TFloatList.Create; V.Count := N + 1;

        t[n] := 1;
        v[n] := 0;

        { z notatek SLE wynika ze t[0], v[0] nie sa potrzebne (i k moze robic
          "downto 1" zamiast "downto 0") ale t[0], v[0] MOGA byc potrzebne:
          przy obliczaniu M[n] dla n = 1. }
        for k := n-1 downto 0 do
        begin
          t[k] := q[k] * t[k+1] + s[k];
          v[k] := q[k] * v[k+1] + u[k];
        end;

        h_n  := X[n] - X[n-1];
        h_n1 := X[1] - X[0];
        delta_n := h_n / (h_n + h_n1);
        M[n] := (PeriodicDN - (1-delta_n) * v[1] - delta_n * v[n-1]) /
                (2          + (1-delta_n) * t[1] + delta_n * t[n-1]);
        M[0] := M[n];
        for k := n-1 downto 1 do  M[k] := v[k] + t[k] * M[n];
      finally
        t.Free;
        v.Free;
      end;
    end else
    begin
      { zawsze M[0] = M[n] = 0, zeby latwo bylo zapisac obliczenia w Evaluate }
      M[0] := 0;
      M[n] := 0;
      M[n-1] := u[n-1];
      for k := n-2 downto 1 do M[k] := u[k] + q[k] * M[k+1];
    end;
  finally
    u.Free;
    q.Free;
    s.Free;
  end;
end;

destructor TNaturalCubicSpline.Destroy;
begin
  if FOwnsX then FX.Free;
  if FOwnsY then FY.Free;
  M.Free;
  inherited;
end;

function TNaturalCubicSpline.Evaluate(x: Float): Float;

  function Power3rd(const a: Float): Float;
  begin
    Result := a * a * a;
  end;

var
  k, KMin, KMax, KMiddle: Cardinal;
  hk: Float;
begin
  Clamp(x, MinX, MaxX);

  { calculate k: W ktorym przedziale x[k-1]..x[k] jest argument ?
    TODO: nalezoloby pomyslec o wykorzystaniu faktu
    ze czesto wiadomo iz wezly x[i] sa rownoodlegle. }
  KMin := 1;
  KMax := FX.Count - 1;
  repeat
    KMiddle:=(KMin + KMax) div 2;
    { jak jest ulozony x w stosunku do przedzialu FX[KMiddle-1]..FX[KMiddle] ? }
    if x < FX[KMiddle-1] then KMax := KMiddle-1 else
      if x > FX[KMiddle] then KMin := KMiddle+1 else
        begin
          KMin := KMiddle; { set only KMin, KMax is meaningless from now }
          break;
        end;
  until KMin = KMax;

  k := KMin;

  Assert(Between(x, FX[k-1], FX[k]));

  { obliczenia uzywaja tych samych symboli co w notatkach SLE }
  { teraz obliczam wartosc s(x) gdzie s to postac funkcji sklejanej
    na przedziale FX[k-1]..FX[k] }
  hk := FX[k] - FX[k-1];
  Result := ( ( M[k-1] * Power3rd(FX[k] - x) + M[k] * Power3rd(x - FX[k-1]) )/6 +
              ( FY[k-1] - M[k-1]*Sqr(hk)/6 )*(FX[k] - x)                        +
              ( FY[k]   - M[k  ]*Sqr(hk)/6 )*(x - FX[k-1])
            ) / hk;
end;

{ TNaturalCubicSplineCurve_Abstract ------------------------------------------- }

procedure TNaturalCubicSplineCurve_Abstract.UpdateControlPoints;
var
  i, j: Integer;
  SplineX, SplineY: TFloatList;
begin
  inherited;

  { calculate SplineX.
    Spline[0] and Spline[1] and Spline[2] will share the same reference to X.
    Only Spline[2] will own SplineX. (Spline[2] will be always Freed as the
    last one, so it's safest to set OwnsX in Spline[2]) }
  SplineX := TFloatList.Create;
  SplineX.Count := ControlPoints.Count;
  for i := 0 to ControlPoints.Count-1 do SplineX[i] := ControlPointT(i);

  for i := 0 to 2 do
  begin
    FreeAndNil(Spline[i]);

    { calculate SplineY }
    SplineY := TFloatList.Create;
    SplineY.Count := ControlPoints.Count;
    for j := 0 to ControlPoints.Count-1 do SplineY[j] := ControlPoints.L[j, i];

    Spline[i] := TNaturalCubicSpline.Create(SplineX, SplineY, i = 2, true, Closed);
  end;
end;

function TNaturalCubicSplineCurve_Abstract.Point(const t: Float): TVector3Single;
var
  i: Integer;
begin
  for i := 0 to 2 do Result[i] := Spline[i].Evaluate(t);
end;

constructor TNaturalCubicSplineCurve_Abstract.Create(AOwner: TComponent);
begin
  inherited;
end;

destructor TNaturalCubicSplineCurve_Abstract.Destroy;
var
  i: Integer;
begin
  for i := 0 to 2 do FreeAndNil(Spline[i]);
  inherited;
end;

{ TNaturalCubicSplineCurve -------------------------------------------------- }

class function TNaturalCubicSplineCurve.NiceClassName: string;
begin
  Result := 'Natural cubic spline curve';
end;

function TNaturalCubicSplineCurve.Closed: boolean;
begin
  Result := VectorsEqual(ControlPoints.First,
                         ControlPoints.Last);
end;

{ TNaturalCubicSplineCurveAlwaysClosed -------------------------------------- }

class function TNaturalCubicSplineCurveAlwaysClosed.NiceClassName: string;
begin
  Result := 'Natural cubic spline curve (closed)';
end;

function TNaturalCubicSplineCurveAlwaysClosed.Closed: boolean;
begin
  Result := true;
end;

{ TNaturalCubicSplineCurveNeverClosed ---------------------------------------- }

class function TNaturalCubicSplineCurveNeverClosed.NiceClassName: string;
begin
  Result := 'Natural cubic spline curve (not closed)';
end;

function TNaturalCubicSplineCurveNeverClosed.Closed: boolean;
begin
  Result := false;
end;

{ TRationalBezierCurve ----------------------------------------------- }

{$define DE_CASTELJAU_DECLARE:=
var
  W: TVector3SingleList;
  Wgh: TFloatList;
  i, k, n, j: Integer;}

{ This initializes W and Wgh (0-th step of de Casteljau algorithm).
  It uses ControlPoints, Weights. }
{$define DE_CASTELJAU_BEGIN:=
  n := ControlPoints.Count - 1;

  W := nil;
  Wgh := nil;
  try
    // using nice FPC memory manager should make this memory allocating
    // (in each call to Point) painless. So I don't care about optimizing
    // this by moving W to private class-scope.
    W := TVector3SingleList.Create;
    W.Assign(ControlPoints);
    Wgh := TFloatList.Create;
    Wgh.Assign(Weights);
}

{ This caculates in W and Wgh k-th step of de Casteljau algorithm.
  This assumes that W and Wgh already contain (k-1)-th step.
  Uses u as the target point position (in [0; 1]) }
{$define DE_CASTELJAU_STEP:=
begin
  for i := 0 to n - k do
  begin
    for j := 0 to 2 do
      W.L[i][j]:=(1-u) * Wgh[i  ] * W.L[i  ][j] +
                         u * Wgh[i+1] * W.L[i+1][j];
    Wgh.L[i]:=(1-u) * Wgh[i] + u * Wgh[i+1];
    for j := 0 to 2 do
      W.L[i][j] /= Wgh.L[i];
  end;
end;}

{ This frees W and Wgh. }
{$define DE_CASTELJAU_END:=
  finally
    Wgh.Free;
    W.Free;
  end;}

procedure TRationalBezierCurve.Split(const Proportion: Float; var B1, B2: TRationalBezierCurve);
var TMiddle, u: Float;
DE_CASTELJAU_DECLARE
begin
  TMiddle := TBegin + Proportion * (TEnd - TBegin);
  {$warnings off} { Consciously using deprecated stuff. }
  B1 := TRationalBezierCurve.Create(nil);
  B2 := TRationalBezierCurve.Create(nil);
  {$warnings on}
  B1.TBegin := TBegin;
  B1.TEnd := TMiddle;
  B2.TBegin := TMiddle;
  B2.TEnd := TEnd;
  B1.ControlPoints.Count := ControlPoints.Count;
  B2.ControlPoints.Count := ControlPoints.Count;
  B1.Weights.Count := Weights.Count;
  B2.Weights.Count := Weights.Count;

  { now we do Casteljau algorithm, similiar to what we do in Point.
    But this time our purpose is to update B1.ControlPoints/Weights and
    B2.ControlPoints/Weights. }

  u := Proportion;

  DE_CASTELJAU_BEGIN
    B1.ControlPoints.L[0] := ControlPoints.L[0];
    B1.Weights      .L[0] := Wgh          .L[0];
    B2.ControlPoints.L[n] := ControlPoints.L[n];
    B2.Weights      .L[n] := Wgh          .L[n];

    for k := 1 to n do
    begin
      DE_CASTELJAU_STEP

      B1.ControlPoints.L[k]   := W  .L[0];
      B1.Weights      .L[k]   := Wgh.L[0];
      B2.ControlPoints.L[n-k] := W  .L[n-k];
      B2.Weights      .L[n-k] := Wgh.L[n-k];
    end;
  DE_CASTELJAU_END
end;

function TRationalBezierCurve.Point(const t: Float): TVector3Single;
var
  u: Float;
DE_CASTELJAU_DECLARE
begin
  { u := t normalized to [0; 1] }
  u := (t - TBegin) / (TEnd - TBegin);

  DE_CASTELJAU_BEGIN
    for k := 1 to n do DE_CASTELJAU_STEP
    Result := W.L[0];
  DE_CASTELJAU_END
end;

class function TRationalBezierCurve.NiceClassName: string;
begin
  Result := 'Rational Bezier curve';
end;

procedure TRationalBezierCurve.UpdateControlPoints;
begin
  inherited;
  Assert(Weights.Count = ControlPoints.Count);
end;

constructor TRationalBezierCurve.Create(AOwner: TComponent);
begin
  inherited;
  Weights := TFloatList.Create;
  Weights.Count := ControlPoints.Count;
end;

destructor TRationalBezierCurve.Destroy;
begin
  Weights.Free;
  inherited;
end;

{ TSmoothInterpolatedCurve ------------------------------------------------------------ }

function TSmoothInterpolatedCurve.CreateConvexHullPoints: TVector3SingleList;
begin
  Result := ConvexHullPoints;
end;

procedure TSmoothInterpolatedCurve.DestroyConvexHullPoints(Points: TVector3SingleList);
begin
end;

function TSmoothInterpolatedCurve.Point(const t: Float): TVector3Single;
var
  i: Integer;
begin
  if ControlPoints.Count = 1 then
    Exit(ControlPoints.L[0]);

  for i := 0 to BezierCurves.Count-1 do
    if t <= BezierCurves[i].TEnd then Break;

  Result := BezierCurves[i].Point(t);
end;

function TSmoothInterpolatedCurve.ToRationalBezierCurves(ResultOwnsCurves: boolean): TRationalBezierCurveList;
var
  S: TVector3SingleList;

  function MiddlePoint(i, Sign: Integer): TVector3Single;
  begin
    Result := ControlPoints.L[i];
    VectorAddTo1st(Result,
      VectorScale(S.L[i], Sign * (ControlPointT(i) - ControlPointT(i-1)) / 3));
  end;

var
  C: TVector3SingleList;
  i: Integer;
  {$warnings off} { Consciously using deprecated stuff. }
  NewCurve: TRationalBezierCurve;
  {$warnings on}
begin
  Result := TRationalBezierCurveList.Create(ResultOwnsCurves);
  try
    if ControlPoints.Count <= 1 then Exit;

    if ControlPoints.Count = 2 then
    begin
      { Normal calcualtions (based on SLE mmgk notes) cannot be done when
        ControlPoints.Count = 2:
        C.Count would be 1, S.Count would be 2,
        S[0] would be calculated based on C[0] and S[1],
        S[1] would be calculated based on C[0] and S[0].
        So I can't calculate S[0] and S[1] using given equations when
        ControlPoints.Count = 2. So I must implement a special case for
        ControlPoints.Count = 2. }
      {$warnings off} { Consciously using deprecated stuff. }
      NewCurve := TRationalBezierCurve.Create(nil);
      {$warnings on}
      NewCurve.TBegin := ControlPointT(0);
      NewCurve.TEnd := ControlPointT(1);
      NewCurve.ControlPoints.Add(ControlPoints.L[0]);
      NewCurve.ControlPoints.Add(Lerp(1/3, ControlPoints.L[0], ControlPoints.L[1]));
      NewCurve.ControlPoints.Add(Lerp(2/3, ControlPoints.L[0], ControlPoints.L[1]));
      NewCurve.ControlPoints.Add(ControlPoints.L[1]);
      NewCurve.Weights.AddArray([1.0, 1.0, 1.0, 1.0]);
      NewCurve.UpdateControlPoints;
      Result.Add(NewCurve);

      Exit;
    end;

    { based on SLE mmgk notes, "Krzywe Beziera" page 4 }
    C := nil;
    S := nil;
    try
      C := TVector3SingleList.Create;
      C.Count := ControlPoints.Count-1;
      { calculate C values }
      for i := 0 to C.Count-1 do
      begin
        C.L[i] := VectorSubtract(ControlPoints.L[i+1], ControlPoints.L[i]);
        VectorScaleTo1st(C.L[i],
          1/(ControlPointT(i+1) - ControlPointT(i)));
      end;

      S := TVector3SingleList.Create;
      S.Count := ControlPoints.Count;
      { calculate S values }
      for i := 1 to S.Count-2 do
        S.L[i] := Lerp( (ControlPointT(i+1) - ControlPointT(i))/
                            (ControlPointT(i+1) - ControlPointT(i-1)),
                            C.L[i-1], C.L[i]);
      S.L[0        ] := VectorSubtract(VectorScale(C.L[0        ], 2), S.L[1        ]);
      S.L[S.Count-1] := VectorSubtract(VectorScale(C.L[S.Count-2], 2), S.L[S.Count-2]);

      for i := 1 to ControlPoints.Count-1 do
      begin
        {$warnings off} { Consciously using deprecated stuff. }
        NewCurve := TRationalBezierCurve.Create(nil);
        {$warnings on}
        NewCurve.TBegin := ControlPointT(i-1);
        NewCurve.TEnd := ControlPointT(i);
        NewCurve.ControlPoints.Add(ControlPoints.L[i-1]);
        NewCurve.ControlPoints.Add(MiddlePoint(i-1, +1));
        NewCurve.ControlPoints.Add(MiddlePoint(i  , -1));
        NewCurve.ControlPoints.Add(ControlPoints.L[i]);
        NewCurve.Weights.AddArray([1.0, 1.0, 1.0, 1.0]);
        NewCurve.UpdateControlPoints;
        Result.Add(NewCurve);
      end;
    finally
      C.Free;
      S.Free;
    end;
  except Result.Free; raise end;
end;

class function TSmoothInterpolatedCurve.NiceClassName: string;
begin
  Result := 'Smooth Interpolated curve';
end;

procedure TSmoothInterpolatedCurve.UpdateControlPoints;
var
  i: Integer;
begin
  inherited;
  FreeAndNil(BezierCurves);

  BezierCurves := ToRationalBezierCurves(true);

  ConvexHullPoints.Clear;
  ConvexHullPoints.AddList(ControlPoints);
  for i := 0 to BezierCurves.Count-1 do
  begin
    ConvexHullPoints.Add(BezierCurves[i].ControlPoints.L[1]);
    ConvexHullPoints.Add(BezierCurves[i].ControlPoints.L[2]);
  end;
end;

constructor TSmoothInterpolatedCurve.Create(AOwner: TComponent);
begin
  inherited;
  ConvexHullPoints := TVector3SingleList.Create;
end;

destructor TSmoothInterpolatedCurve.Destroy;
begin
  FreeAndNil(BezierCurves);
  FreeAndNil(ConvexHullPoints);
  inherited;
end;

end.