GLScene/Externals/Graphics32/Source/GR32_Transforms.pas

unit GR32_Transforms;

(* ***** BEGIN LICENSE BLOCK *****
 * Version: MPL 1.1 or LGPL 2.1 with linking exception
 *
 * 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.
 *
 * Alternatively, the contents of this file may be used under the terms of the
 * Free Pascal modified version of the GNU Lesser General Public License
 * Version 2.1 (the "FPC modified LGPL License"), in which case the provisions
 * of this license are applicable instead of those above.
 * Please see the file LICENSE.txt for additional information concerning this
 * license.
 *
 * The Original Code is Graphics32
 *
 * The Initial Developer of the Original Code is
 * Alex A. Denisov
 *
 * Portions created by the Initial Developer are Copyright (C) 2000-2009
 * the Initial Developer. All Rights Reserved.
 *
 * ***** END LICENSE BLOCK ***** *)

interface

{$include GR32.inc}

uses
{$if defined(MSWINDOWS) and defined(DEBUG)}
  Windows, // For OutputDebugString. In interface section so we don't override TFixed
{$ifend}
  SysUtils,
  Classes,
  Types,
  GR32,
  GR32_VectorMaps,
  GR32_Rasterizers;

type
  ETransformError = class(Exception);
  ETransformNotImplemented = class(Exception);

type
  TFloatMatrix = array [0..2, 0..2] of TFloat;     // 3x3 TFloat precision
  TFixedMatrix = array [0..2, 0..2] of TFixed;     // 3x3 fixed precision

const
  IdentityMatrix: TFloatMatrix = (
    (1, 0, 0),
    (0, 1, 0),
    (0, 0, 1));

type
  TVector3f = array [0..2] of TFloat;
  TVector3i = array [0..2] of Integer;

// Matrix conversion routines
function FixedMatrix(const FloatMatrix: TFloatMatrix): TFixedMatrix; overload;
function FloatMatrix(const FixedMatrix: TFixedMatrix): TFloatMatrix; overload;

procedure Adjoint(var M: TFloatMatrix);
function Determinant(const M: TFloatMatrix): TFloat;
procedure Scale(var M: TFloatMatrix; Factor: TFloat);
procedure Invert(var M: TFloatMatrix);
function Mult(const M1, M2: TFloatMatrix): TFloatMatrix;
function VectorTransform(const M: TFloatMatrix; const V: TVector3f): TVector3f;


//------------------------------------------------------------------------------
//
//      TTransformation
//
//------------------------------------------------------------------------------
type
  TTransformation = class(TNotifiablePersistent)
  private
    FSrcRect: TFloatRect;
    procedure SetSrcRect(const Value: TFloatRect);
  protected
    TransformValid: Boolean;
    procedure PrepareTransform; virtual;
    procedure ReverseTransformInt(DstX, DstY: Integer; out SrcX, SrcY: Integer); virtual;
    procedure ReverseTransformFixed(DstX, DstY: TFixed; out SrcX, SrcY: TFixed); virtual;
    procedure ReverseTransformFloat(DstX, DstY: TFloat; out SrcX, SrcY: TFloat); virtual;
    procedure TransformInt(SrcX, SrcY: Integer; out DstX, DstY: Integer); virtual;
    procedure TransformFixed(SrcX, SrcY: TFixed; out DstX, DstY: TFixed); virtual;
    procedure TransformFloat(SrcX, SrcY: TFloat; out DstX, DstY: TFloat); virtual;
  public
    constructor Create; virtual;
    procedure Changed; override;
    function HasTransformedBounds: Boolean; virtual;
    function GetTransformedBounds: TFloatRect; overload;
    function GetTransformedBounds(const ASrcRect: TFloatRect): TFloatRect; overload; virtual;
    function ReverseTransform(const P: TPoint): TPoint; overload; virtual;
    function ReverseTransform(const P: TFixedPoint): TFixedPoint; overload; virtual;
    function ReverseTransform(const P: TFloatPoint): TFloatPoint; overload; virtual;
    function Transform(const P: TPoint): TPoint; overload; virtual;
    function Transform(const P: TFixedPoint): TFixedPoint; overload; virtual;
    function Transform(const P: TFloatPoint): TFloatPoint; overload; virtual;
    property SrcRect: TFloatRect read FSrcRect write SetSrcRect;
  end;

  TTransformationClass = class of TTransformation;


//------------------------------------------------------------------------------
//
//      TNestedTransformation
//
//------------------------------------------------------------------------------
type
  TNestedTransformation = class(TTransformation)
  private
    FItems: TList;
    FOwner: TPersistent;
    function GetCount: Integer;
    function GetItem(Index: Integer): TTransformation;
    procedure SetItem(Index: Integer; const Value: TTransformation);
  protected
    procedure PrepareTransform; override;
    procedure ReverseTransformFixed(DstX, DstY: TFixed; out SrcX, SrcY: TFixed); override;
    procedure ReverseTransformFloat(DstX, DstY: TFloat; out SrcX, SrcY: TFloat); override;
    procedure TransformFixed(SrcX, SrcY: TFixed; out DstX, DstY: TFixed); override;
    procedure TransformFloat(SrcX, SrcY: TFloat; out DstX, DstY: TFloat); override;
  public
    constructor Create; override;
    destructor Destroy; override;
    function Add(ItemClass: TTransformationClass): TTransformation;
    procedure Clear;
    procedure Delete(Index: Integer);
    function Insert(Index: Integer; ItemClass: TTransformationClass): TTransformation;

    property Owner: TPersistent read FOwner;
    property Count: Integer read GetCount;
    property Items[Index: Integer]: TTransformation read GetItem write SetItem; default;
  end;


//------------------------------------------------------------------------------
//
//      T3x3Transformation
//
//------------------------------------------------------------------------------
type
  T3x3Transformation = class(TTransformation)
  protected
    FMatrix, FInverseMatrix: TFloatMatrix;
    FFixedMatrix, FInverseFixedMatrix: TFixedMatrix;
    procedure PrepareTransform; override;
    procedure ReverseTransformFixed(DstX, DstY: TFixed; out SrcX, SrcY: TFixed); override;
    procedure ReverseTransformFloat(DstX, DstY: TFloat; out SrcX, SrcY: TFloat); override;
    procedure TransformFixed(SrcX, SrcY: TFixed; out DstX, DstY: TFixed); override;
    procedure TransformFloat(SrcX, SrcY: TFloat; out DstX, DstY: TFloat); override;
  public
    property Matrix: TFloatMatrix read FMatrix;
  end;


//------------------------------------------------------------------------------
//
//      TAffineTransformation
//
//------------------------------------------------------------------------------
type
  TAffineTransformation = class(T3x3Transformation)
  private
    FStack: ^TFloatMatrix;
    FStackLevel: Integer;
  public
    constructor Create; override;

    function GetTransformedBounds(const ASrcRect: TFloatRect): TFloatRect; override;
    procedure Push;
    procedure Pop;
    procedure Clear; overload;
    procedure Clear(BaseMatrix: TFloatMatrix); overload;
    procedure Rotate(Alpha: TFloat); overload; // degrees
    procedure Rotate(Cx, Cy, Alpha: TFloat); overload; // degrees
    procedure Skew(Fx, Fy: TFloat);
    procedure Scale(Sx, Sy: TFloat); overload;
    procedure Scale(Value: TFloat); overload;
    procedure Translate(Dx, Dy: TFloat);
  end;


//------------------------------------------------------------------------------
//
//      TProjectiveTransformation
//
//------------------------------------------------------------------------------
type
  TProjectiveTransformation = class(T3x3Transformation)
  private
    FQuadX: array [0..3] of TFloat;
    FQuadY: array [0..3] of TFloat;
    procedure SetX(Index: Integer; const Value: TFloat); {$IFDEF UseInlining} inline; {$ENDIF}
    procedure SetY(Index: Integer; const Value: TFloat); {$IFDEF UseInlining} inline; {$ENDIF}
    function GetX(Index: Integer): TFloat; {$IFDEF UseInlining} inline; {$ENDIF}
    function GetY(Index: Integer): TFloat; {$IFDEF UseInlining} inline; {$ENDIF}
  protected
    procedure PrepareTransform; override;
    procedure ReverseTransformFixed(DstX, DstY: TFixed; out SrcX, SrcY: TFixed); override;
    procedure ReverseTransformFloat(DstX, DstY: TFloat; out SrcX, SrcY: TFloat); override;
    procedure TransformFixed(SrcX, SrcY: TFixed; out DstX, DstY: TFixed); override;
    procedure TransformFloat(SrcX, SrcY: TFloat; out DstX, DstY: TFloat); override;
  public
    function GetTransformedBounds(const ASrcRect: TFloatRect): TFloatRect; override;
    property X[Index: Integer]: TFloat read GetX write SetX;
    property Y[index: Integer]: TFloat read GetX write SetY;
  published
    property X0: TFloat index 0 read GetX write SetX;
    property X1: TFloat index 1 read GetX write SetX;
    property X2: TFloat index 2 read GetX write SetX;
    property X3: TFloat index 3 read GetX write SetX;
    property Y0: TFloat index 0 read GetY write SetY;
    property Y1: TFloat index 1 read GetY write SetY;
    property Y2: TFloat index 2 read GetY write SetY;
    property Y3: TFloat index 3 read GetY write SetY;
  end;


//------------------------------------------------------------------------------
//
//      TProjectiveTransformationEx
//
//------------------------------------------------------------------------------
// Performs projective transformation between two convex quadrilaterals.
//------------------------------------------------------------------------------
// References:
//
// - "Fundamentals of Texture Mapping and Image Warping"
//   Paul S. Heckbert
//   www.cs.cmu.edu/~ph/texfund/texfund.pdf
//
// - "Projective Mappings for ImageWarping"
//   Paul S. Heckbert
//   http://graphics.cs.cmu.edu/courses/15-463/2008_fall/Papers/proj.pdf
//
// - "Geometric Tools for Computer Graphics"
//   David H. Eberly
//   https://www.amazon.com/Geometric-Computer-Graphics-Morgan-Kaufmann/dp/1558605940
//
// - "Perspective Mappings"
//   David H. Eberly
//   https://geometrictools.com/Documentation/PerspectiveMappings.pdf
//
//------------------------------------------------------------------------------
{$ifndef FPC} // FPC considers unicode identifiers a bug, so f*ck them and their silly toy compiler.
type
  TQuadrilateral = array[0..3] of TPoint;
  TFloatQuadrilateral = array[0..3] of TFloatPoint;

  TProjectiveTransformationEx = class(T3x3Transformation)
  private
    FExtrapolate: boolean;
    FSourceQuad: TFloatQuadrilateral;
    FDestQuad: TFloatQuadrilateral;
    procedure SetSourceQuad(const Value: TFloatQuadrilateral); {$IFDEF UseInlining} inline; {$ENDIF}
    procedure SetDestQuad(const Value: TFloatQuadrilateral); {$IFDEF UseInlining} inline; {$ENDIF}
    procedure SetSource(Index: Integer; const Value: TFloatPoint); {$IFDEF UseInlining} inline; {$ENDIF}
    procedure SetSourceX(Index: Integer; const Value: TFloat); {$IFDEF UseInlining} inline; {$ENDIF}
    procedure SetSourceY(Index: Integer; const Value: TFloat); {$IFDEF UseInlining} inline; {$ENDIF}
    procedure SetDest(Index: Integer; const Value: TFloatPoint); {$IFDEF UseInlining} inline; {$ENDIF}
    procedure SetDestX(Index: Integer; const Value: TFloat); {$IFDEF UseInlining} inline; {$ENDIF}
    procedure SetDestY(Index: Integer; const Value: TFloat); {$IFDEF UseInlining} inline; {$ENDIF}
    function GetSource(Index: Integer): TFloatPoint; {$IFDEF UseInlining} inline; {$ENDIF}
    function GetSourceX(Index: Integer): TFloat; {$IFDEF UseInlining} inline; {$ENDIF}
    function GetSourceY(Index: Integer): TFloat; {$IFDEF UseInlining} inline; {$ENDIF}
    function GetDest(Index: Integer): TFloatPoint; {$IFDEF UseInlining} inline; {$ENDIF}
    function GetDestX(Index: Integer): TFloat; {$IFDEF UseInlining} inline; {$ENDIF}
    function GetDestY(Index: Integer): TFloat; {$IFDEF UseInlining} inline; {$ENDIF}
  protected
    procedure PrepareTransform; override;
    procedure ReverseTransformFixed(DstX, DstY: TFixed; out SrcX, SrcY: TFixed); override;
    procedure ReverseTransformFloat(DstX, DstY: TFloat; out SrcX, SrcY: TFloat); override;
    procedure TransformFixed(SrcX, SrcY: TFixed; out DstX, DstY: TFixed); override;
    procedure TransformFloat(SrcX, SrcY: TFloat; out DstX, DstY: TFloat); override;
  public
    function GetTransformedBounds(const ASrcRect: TFloatRect): TFloatRect; override;

    property SourceQuad: TFloatQuadrilateral read FSourceQuad write SetSourceQuad;
    property Source[Index: Integer]: TFloatPoint read GetSource write SetSource;
    property SourceX[Index: Integer]: TFloat read GetSourceX write SetSourceX;
    property SourceY[Index: Integer]: TFloat read GetSourceY write SetSourceY;
    property DestQuad: TFloatQuadrilateral read FDestQuad write SetDestQuad;
    property Dest[Index: Integer]: TFloatPoint read GetDest write SetDest;
    property DestX[Index: Integer]: TFloat read GetDestX write SetDestX;
    property DestY[index: Integer]: TFloat read GetDestX write SetDestY;
  published
    // Set Extrapolate=True to have pixels beyond the destination be transformed.
    // This is done by having GetTransformedBounds return the passed source rect,
    // which in turn causes the rasterizer to process all pixels of the source
    // image instead of just the pixels covered by the target quad.
    property Extrapolate: boolean read FExtrapolate write FExtrapolate;

    property DestX0: TFloat index 0 read GetDestX write SetDestX;
    property DestX1: TFloat index 1 read GetDestX write SetDestX;
    property DestX2: TFloat index 2 read GetDestX write SetDestX;
    property DestX3: TFloat index 3 read GetDestX write SetDestX;
    property DestY0: TFloat index 0 read GetDestY write SetDestY;
    property DestY1: TFloat index 1 read GetDestY write SetDestY;
    property DestY2: TFloat index 2 read GetDestY write SetDestY;
    property DestY3: TFloat index 3 read GetDestY write SetDestY;

    property SourceX0: TFloat index 0 read GetSourceX write SetSourceX;
    property SourceX1: TFloat index 1 read GetSourceX write SetSourceX;
    property SourceX2: TFloat index 2 read GetSourceX write SetSourceX;
    property SourceX3: TFloat index 3 read GetSourceX write SetSourceX;
    property SourceY0: TFloat index 0 read GetSourceY write SetSourceY;
    property SourceY1: TFloat index 1 read GetSourceY write SetSourceY;
    property SourceY2: TFloat index 2 read GetSourceY write SetSourceY;
    property SourceY3: TFloat index 3 read GetSourceY write SetSourceY;
  end;
{$endif FPC}


//------------------------------------------------------------------------------
//
//      TTwirlTransformation
//
//------------------------------------------------------------------------------
type
  TTwirlTransformation = class(TTransformation)
  private
    Frx, Fry: TFloat;
    FTwirl: TFloat;
    procedure SetTwirl(const Value: TFloat);
  protected
    procedure PrepareTransform; override;
    procedure ReverseTransformFloat(DstX, DstY: TFloat; out SrcX, SrcY: TFloat); override;
  public
    constructor Create; override;
    function GetTransformedBounds(const ASrcRect: TFloatRect): TFloatRect; override;
  published
    property Twirl: TFloat read FTwirl write SetTwirl;
  end;


//------------------------------------------------------------------------------
//
//      TBloatTransformation
//
//------------------------------------------------------------------------------
type
  TBloatTransformation = class(TTransformation)
  private
    FBloatPower: TFloat;
    FBP: TFloat;
    FPiW, FPiH: TFloat;
    procedure SetBloatPower(const Value: TFloat);
  protected
    procedure PrepareTransform; override;
    procedure ReverseTransformFloat(DstX, DstY: TFloat; out SrcX, SrcY: TFloat); override;
    procedure TransformFloat(DstX, DstY: TFloat; out SrcX, SrcY: TFloat); override;
  public
    constructor Create; override;
  published
    property BloatPower: TFloat read FBloatPower write SetBloatPower;
  end;


//------------------------------------------------------------------------------
//
//      TDisturbanceTransformation
//
//------------------------------------------------------------------------------
type
  TDisturbanceTransformation = class(TTransformation)
  private
    FDisturbance: TFloat;
    procedure SetDisturbance(const Value: TFloat);
  protected
    procedure ReverseTransformFloat(DstX, DstY: TFloat; out SrcX, SrcY: TFloat); override;
  public
    function GetTransformedBounds(const ASrcRect: TFloatRect): TFloatRect; override;
  published
    property Disturbance: TFloat read FDisturbance write SetDisturbance;
  end;


//------------------------------------------------------------------------------
//
//      TFishEyeTransformation
//
//------------------------------------------------------------------------------
type
  TFishEyeTransformation = class(TTransformation)
  private
    Frx, Fry: TFloat;
    Faw, Fsr: TFloat;
    Sx, Sy: TFloat;
    FMinR: TFloat;
  protected
    procedure PrepareTransform; override;
    procedure ReverseTransformFloat(DstX, DstY: TFloat; out SrcX, SrcY: TFloat); override;
  end;


//------------------------------------------------------------------------------
//
//      TPolarTransformation
//
//------------------------------------------------------------------------------
type
  TPolarTransformation = class(TTransformation)
  private
    FDstRect: TFloatRect;
    FPhase: TFloat;
    Sx, Sy, Cx, Cy, Dx, Dy, Rt, Rt2, Rr, Rcx, Rcy: TFloat;
    procedure SetDstRect(const Value: TFloatRect);
    procedure SetPhase(const Value: TFloat);
  protected
    procedure PrepareTransform; override;
    procedure TransformFloat(SrcX, SrcY: TFloat; out DstX, DstY: TFloat); override;
    procedure ReverseTransformFloat(DstX, DstY: TFloat; out SrcX, SrcY: TFloat); override;
  public
    property DstRect: TFloatRect read FDstRect write SetDstRect;
    property Phase: TFloat read FPhase write SetPhase;
  end;


//------------------------------------------------------------------------------
//
//      TPathTransformation
//
//------------------------------------------------------------------------------
type
  TPathTransformation = class(TTransformation)
  private
    FTopLength: TFloat;
    FBottomLength: TFloat;
    FBottomCurve: TArrayOfFloatPoint;
    FTopCurve: TArrayOfFloatPoint;
    FTopHypot, FBottomHypot: array of record Dist, RecDist: TFloat end;
    procedure SetBottomCurve(const Value: TArrayOfFloatPoint);
    procedure SetTopCurve(const Value: TArrayOfFloatPoint);
  protected
    rdx, rdy: TFloat;
    procedure PrepareTransform; override;
    procedure TransformFloat(SrcX, SrcY: TFloat; out DstX, DstY: TFloat); override;
  public
    destructor Destroy; override;
    property TopCurve: TArrayOfFloatPoint read FTopCurve write SetTopCurve;
    property BottomCurve: TArrayOfFloatPoint read FBottomCurve write SetBottomCurve;
  end;


//------------------------------------------------------------------------------
//
//      TRadialDistortionTransformation
//
//------------------------------------------------------------------------------
type
  TRadialDistortionTransformation = class(TTransformation)
  protected
    FCoefficient1, FCoefficient2, FScale: TFloat;
    FFocalPoint: TFloatPoint;
    r_0, r_tgt_max, r_tgt_min: Single;
    FMapElements: Integer;
    Map: Array of TFloat;
    function LookUpReverseMap(const r_tgt: TFloat): TFloat;
    procedure SetCoefficient1(const Value: TFloat);
    procedure SetCoefficient2(const Value: TFloat);
    procedure SetScale(const Value: TFloat);
    procedure SetMapElements(const Value: Integer);
    procedure PrepareReverseMap;
    procedure PrepareTransform; override;
    procedure ReverseTransformFloat(DstX, DstY: TFloat; out SrcX, SrcY: TFloat); override;
    procedure TransformFloat(SrcX, SrcY: TFloat; out DstX, DstY: TFloat); override;
  public
    constructor Create; override;
    function HasTransformedBounds: Boolean; override;
  published
    property Coefficient1: TFloat read FCoefficient1 write SetCoefficient1;
    property Coefficient2: TFloat read FCoefficient2 write SetCoefficient2;
    property Scale: TFloat read FScale write SetScale;
    property MapElements: Integer read FMapElements write SetMapElements;
  end;


//------------------------------------------------------------------------------
//
//      TRemapTransformation
//
//------------------------------------------------------------------------------
type
  TRemapTransformation = class(TTransformation)
  private
    FVectorMap: TVectorMap;
    FScalingFixed: TFixedVector;
    FScalingFloat: TFloatVector;
    FCombinedScalingFixed: TFixedVector;
    FCombinedScalingFloat: TFloatVector;
    FSrcTranslationFixed: TFixedVector;
    FSrcScaleFixed: TFixedVector;
    FDstTranslationFixed: TFixedVector;
    FDstScaleFixed: TFixedVector;
    FSrcTranslationFloat: TFloatVector;
    FSrcScaleFloat: TFloatVector;
    FDstTranslationFloat: TFloatVector;
    FDstScaleFloat: TFloatVector;
    FOffsetFixed : TFixedVector;
    FOffsetInt : TPoint;
    FMappingRect: TFloatRect;
    FOffset: TFloatVector;
    procedure SetMappingRect(Rect: TFloatRect);
    procedure SetOffset(const Value: TFloatVector);
  protected
    procedure PrepareTransform; override;
    procedure ReverseTransformInt(DstX, DstY: Integer; out SrcX, SrcY: Integer); override;
    procedure ReverseTransformFloat(DstX, DstY: TFloat; out SrcX, SrcY: TFloat); override;
    procedure ReverseTransformFixed(DstX, DstY: TFixed; out SrcX, SrcY: TFixed); override;
  public
    constructor Create; override;
    destructor Destroy; override;
    function HasTransformedBounds: Boolean; override;
    function GetTransformedBounds(const ASrcRect: TFloatRect): TFloatRect; override;
    procedure Scale(Sx, Sy: TFloat);
    property MappingRect: TFloatRect read FMappingRect write SetMappingRect;
    property Offset: TFloatVector read FOffset write SetOffset;
    property VectorMap: TVectorMap read FVectorMap write FVectorMap;
  end;


//------------------------------------------------------------------------------
//
//      TSphereTransformation
//
//------------------------------------------------------------------------------
// Transform a map (planisphere) into a Spherical projection.
// By Marc LAFON (marc.lafon AT free.fr), 01 nov 2005
//------------------------------------------------------------------------------
type
  TSphereTransformation = class(TTransformation)
  private
    FMapWidth, FMapHeight: TFloat;
    FSquareRadius: TFloat;
    FCenter: TFloatPoint;
    FRadius: TFloat;
    FLongitude: TFloat;
    FLattitude: TFloat;
    FLattitudeSin: TFloat;
    FLattitudeCos: TFloat;
    FLattitudeSinInvRadius: TFloat;
    FLattitudeCosInvRadius: TFloat;
    FSrcRectTop: TFloat;
    FSrcRectLeft: TFloat;
    procedure SetCenter(const Value: TFloatPoint);
    procedure SetLattitude(const Value: TFloat);
    procedure SetLongitude(const Value: TFloat);
    procedure SetRadius(const Value: TFloat);
  protected
    procedure PrepareTransform; override;
    procedure ReverseTransformFloat(DstX, DstY: TFloat; out SrcX, SrcY: TFloat); override;
  public
    constructor Create; override;

    function HasTransformedBounds: Boolean; override;
    function GetTransformedBounds(const ASrcRect: TFloatRect): TFloatRect; override;

    // Return True if the (X,Y) point is in the Sphere projection
    function IsInSphere(CartesianX, CartesianY: TFloat):boolean;
    // Transform (X,Y) coordinate as Lattitude and Longitude coordinates in the Sphere
    function SphericalCoordinate(CartesianX, CartesianY: TFloat):TFloatPoint;
    // Transform Longitude and Lattitude coordinates (X,Y) into their screen projection.
    // Returns False if this point is on visible face.
    function ScreenCoordinate(var X, Y: TFloat):boolean;

    // Center of the Sphere in the Destination Bitmap
    property Center: TFloatPoint read FCenter write SetCenter;
    // Radius of the Sphere in the Destination Bitmap
    property Radius: TFloat read FRadius write SetRadius;
    // Rotation of the Sphere (Y-axe rotation angle)
    property Lattitude: TFloat read FLattitude write SetLattitude;
    // Rotation of the Sphere (X-axe rotation angle)
    property Longitude: TFloat read FLongitude write SetLongitude;
  end;


//------------------------------------------------------------------------------
//
//      Utilities
//
//------------------------------------------------------------------------------
function TransformPoints(Points: TArrayOfArrayOfFixedPoint; Transformation: TTransformation): TArrayOfArrayOfFixedPoint;

procedure Transform(Dst, Src: TCustomBitmap32; Transformation: TTransformation; Reverse: boolean = True); overload;
procedure Transform(Dst, Src: TCustomBitmap32; Transformation: TTransformation; const DstClip: TRect; Reverse: boolean = True); overload;
procedure Transform(Dst, Src: TCustomBitmap32; Transformation: TTransformation; Rasterizer: TRasterizer; Reverse: boolean = True); overload;
procedure Transform(Dst, Src: TCustomBitmap32; Transformation: TTransformation; Rasterizer: TRasterizer; const DstClip: TRect; Reverse: boolean = True); overload;

procedure RasterizeTransformation(Vectormap: TVectormap;
  Transformation: TTransformation; DstRect: TRect;
  CombineMode: TVectorCombineMode = vcmAdd;
  CombineCallback: TVectorCombineEvent = nil);

procedure SetBorderTransparent(ABitmap: TCustomBitmap32; ARect: TRect);

//------------------------------------------------------------------------------

{ FullEdge controls how the bitmap is resampled }
var
  FullEdge: Boolean = True;

resourcestring
  RCStrReverseTransformationNotImplemented = 'Reverse transformation is not implemented in %s.';
  RCStrForwardTransformationNotImplemented = 'Forward transformation is not implemented in %s.';
  RCStrTopBottomCurveNil = 'Top or bottom curve is nil';


//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------

implementation

uses
  Math,
  GR32_Blend,
  GR32_LowLevel,
  GR32_Math,
  GR32_Bindings,
  GR32_Resamplers,
  GR32_Geometry;

resourcestring
  RCStrSrcRectIsEmpty = 'SrcRect is empty!';
  RCStrMappingRectIsEmpty = 'MappingRect is empty!';
  RStrStackEmpty = 'Stack empty';

type
  { provides access to proctected members of TTransformation by typecasting }
  TTransformationAccess = class(TTransformation);

//------------------------------------------------------------------------------
//
//      A bit of linear algebra
//
//------------------------------------------------------------------------------

var
  DET_2x2_32: function(a1, a2, b1, b2: TFloat): TFloat;
  DET_3x3_32: function(a1, a2, a3, b1, b2, b3, c1, c2, c3: TFloat): TFloat;
  DET_2x2_64: function(a1, a2, b1, b2: Double): Double;


//------------------------------------------------------------------------------
// DET_2x2_32
//------------------------------------------------------------------------------
function DET_2x2_32_Pas(a1, a2, b1, b2: TFloat): TFloat;
begin
  Result := a1 * b2 - a2 * b1;
end;

{$IFNDEF PUREPASCAL}

{$if defined(TARGET_x86)}

function DET_2x2_32_ASM(a1, a2, b1, b2: TFloat): TFloat; {$IFDEF FPC}assembler; {$ENDIF}
asm
        FLD     A1.Single
        FMUL    B2.Single
        FLD     A2.Single
        FMUL    B1.Single
        FSUBP
end;

{$elseif defined(TARGET_x64)}

function DET_2x2_32_SSE(a1, a2, b1, b2: TFloat): TFloat; {$IFDEF FPC}assembler; nostackframe;{$ENDIF}
asm
  // XMM0: a1
  // XMM1: a2
  // XMM2: b1
  // XMM3: b2

        MULSS   XMM0, XMM3      // XMM0 <- a1 * b2
        MULSS   XMM1, XMM2      // XMM1 <- a2 * b1
        SUBSS   XMM0, XMM1      // Result <- (a1 * b2) - (a2 * b1)
end;

{$ifend}

{$ENDIF}


//------------------------------------------------------------------------------
// DET_2x2_64
//------------------------------------------------------------------------------
function DET_2x2_64_Pas(a1, a2, b1, b2: Double): Double; overload;
begin
  Result := a1 * b2 - a2 * b1;
end;

{$IFNDEF PUREPASCAL}

{$if defined(TARGET_x86)}

function DET_2x2_64_ASM(a1, a2, b1, b2: Double): Double; {$IFDEF FPC}assembler; {$ENDIF}
asm
        FLD     A1.Double
        FMUL    B2.Double
        FLD     A2.Double
        FMUL    B1.Double
        FSUBP
end;

{$elseif defined(TARGET_x64)}

function DET_2x2_64_SSE(a1, a2, b1, b2: Double): Double; {$IFDEF FPC}assembler; nostackframe;{$ENDIF}
asm
  // XMM0: a1
  // XMM1: a2
  // XMM2: b1
  // XMM3: b2

        MULSD   XMM0, XMM3      // XMM0 <- a1 * b2
        MULSD   XMM1, XMM2      // XMM1 <- a2 * b1
        SUBSD   XMM0, XMM1      // Result <- (a1 * b2) - (a2 * b1)
end;

{$ifend}

{$ENDIF}


//------------------------------------------------------------------------------
// DET_3x3_32
//------------------------------------------------------------------------------
function DET_3x3_32_Pas(a1, a2, a3, b1, b2, b3, c1, c2, c3: TFloat): TFloat; overload; {$IFDEF UseInlining} inline; {$ENDIF}
begin
  Result :=
    a1 * (b2 * c3 - b3 * c2) -
    b1 * (a2 * c3 - a3 * c2) +
    c1 * (a2 * b3 - a3 * b2);
end;


//------------------------------------------------------------------------------
//
//      Utilities
//
//------------------------------------------------------------------------------
procedure Adjoint(var M: TFloatMatrix);
var
  Tmp: TFloatMatrix;
begin
  Tmp := M;

  M[0,0] :=  DET_2x2_32(Tmp[1,1], Tmp[1,2], Tmp[2,1], Tmp[2,2]);
  M[0,1] := -DET_2x2_32(Tmp[0,1], Tmp[0,2], Tmp[2,1], Tmp[2,2]);
  M[0,2] :=  DET_2x2_32(Tmp[0,1], Tmp[0,2], Tmp[1,1], Tmp[1,2]);

  M[1,0] := -DET_2x2_32(Tmp[1,0], Tmp[1,2], Tmp[2,0], Tmp[2,2]);
  M[1,1] :=  DET_2x2_32(Tmp[0,0], Tmp[0,2], Tmp[2,0], Tmp[2,2]);
  M[1,2] := -DET_2x2_32(Tmp[0,0], Tmp[0,2], Tmp[1,0], Tmp[1,2]);

  M[2,0] :=  DET_2x2_32(Tmp[1,0], Tmp[1,1], Tmp[2,0], Tmp[2,1]);
  M[2,1] := -DET_2x2_32(Tmp[0,0], Tmp[0,1], Tmp[2,0], Tmp[2,1]);
  M[2,2] :=  DET_2x2_32(Tmp[0,0], Tmp[0,1], Tmp[1,0], Tmp[1,1]);
end;

//------------------------------------------------------------------------------

function Determinant(const M: TFloatMatrix): TFloat;
begin
  Result := DET_3x3_32_Pas(
    M[0,0], M[1,0], M[2,0],
    M[0,1], M[1,1], M[2,1],
    M[0,2], M[1,2], M[2,2]);
end;

//------------------------------------------------------------------------------

procedure Scale(var M: TFloatMatrix; Factor: TFloat);
var
  i, j: Integer;
begin
  for i := 0 to 2 do
    for j := 0 to 2 do
      M[i,j] := M[i,j] * Factor;
end;

//------------------------------------------------------------------------------

procedure Invert(var M: TFloatMatrix);
var
  Det: TFloat;
begin
  Det := Determinant(M);
  if Abs(Det) < 1E-5 then
    M := IdentityMatrix
  else
  begin
    Adjoint(M);
    Scale(M, 1 / Det);
  end;
end;

//------------------------------------------------------------------------------

function Mult(const M1, M2: TFloatMatrix): TFloatMatrix;
var
  i, j: Integer;
begin
  for i := 0 to 2 do
    for j := 0 to 2 do
      Result[i, j] :=
        M1[0, j] * M2[i, 0] +
        M1[1, j] * M2[i, 1] +
        M1[2, j] * M2[i, 2];
end;

//------------------------------------------------------------------------------

function VectorTransform(const M: TFloatMatrix; const V: TVector3f): TVector3f;
begin
  Result[0] := M[0,0] * V[0] + M[1,0] * V[1] + M[2,0] * V[2];
  Result[1] := M[0,1] * V[0] + M[1,1] * V[1] + M[2,1] * V[2];
  Result[2] := M[0,2] * V[0] + M[1,2] * V[1] + M[2,2] * V[2];
end;

//------------------------------------------------------------------------------

procedure SetBorderTransparent(ABitmap: TCustomBitmap32; ARect: TRect);
var
  I: Integer;
begin
  GR32.IntersectRect(ARect, ARect, ABitmap.BoundsRect);
  with ARect, ABitmap do
  if (Right > Left) and (Bottom > Top) and
    (Left < ClipRect.Right) and (Top < ClipRect.Bottom) and
    (Right > ClipRect.Left) and (Bottom > ClipRect.Top) then
  begin
    Dec(Right);
    Dec(Bottom);
    for I := Left to Right do
      begin
      ABitmap[I, Top] := ABitmap[I, Top] and $00FFFFFF;
      ABitmap[I, Bottom] := ABitmap[I, Bottom] and $00FFFFFF;
      end;
    for I := Top to Bottom do
    begin
      ABitmap[Left, I] := ABitmap[Left, I] and $00FFFFFF;
      ABitmap[Right, I] := ABitmap[Right, I] and $00FFFFFF;
    end;
    Changed;
  end;
end;

//------------------------------------------------------------------------------
// Transformation functions
//------------------------------------------------------------------------------
function TransformPoints(Points: TArrayOfArrayOfFixedPoint; Transformation: TTransformation): TArrayOfArrayOfFixedPoint;
var
  I, J: Integer;
begin
  if Points = nil then
    Result := nil
  else
  begin
    SetLength(Result, Length(Points));
    Transformation.PrepareTransform;

    for I := 0 to High(Result) do
    begin
      SetLength(Result[I], Length(Points[I]));
      if Length(Result[I]) > 0 then
        for J := 0 to High(Result[I]) do
          Transformation.TransformFixed(Points[I][J].X, Points[I][J].Y, Result[I][J].X, Result[I][J].Y);
    end;
  end;
end;

//------------------------------------------------------------------------------

procedure Transform(Dst, Src: TCustomBitmap32; Transformation: TTransformation; Reverse: boolean);
var
  Rasterizer: TRasterizer;
begin
  Rasterizer := DefaultRasterizerClass.Create;
  try
    Transform(Dst, Src, Transformation, Rasterizer, Reverse);
  finally
    Rasterizer.Free;
  end;
end;

//------------------------------------------------------------------------------

procedure Transform(Dst, Src: TCustomBitmap32; Transformation: TTransformation; const DstClip: TRect; Reverse: boolean);
var
  Rasterizer: TRasterizer;
begin
  Rasterizer := DefaultRasterizerClass.Create;
  try
    Transform(Dst, Src, Transformation, Rasterizer, DstClip, Reverse);
  finally
    Rasterizer.Free;
  end;
end;

//------------------------------------------------------------------------------

procedure Transform(Dst, Src: TCustomBitmap32; Transformation: TTransformation;
  Rasterizer: TRasterizer; Reverse: boolean);
begin
  Transform(Dst, Src, Transformation, Rasterizer, Dst.BoundsRect, Reverse);
end;

//------------------------------------------------------------------------------

procedure Transform(Dst, Src: TCustomBitmap32; Transformation: TTransformation;
  Rasterizer: TRasterizer; const DstClip: TRect; Reverse: boolean);
var
  DstRect: TRect;
  Transformer: TTransformer;
begin
  GR32.IntersectRect(DstRect, DstClip, Dst.ClipRect);

  if (DstRect.Right < DstRect.Left) or (DstRect.Bottom < DstRect.Top) then
    Exit;

  if not Dst.MeasuringMode then
  begin
    Transformer := TTransformer.Create(Src.Resampler, Transformation, Reverse);
    try
      Rasterizer.Sampler := Transformer;
      Rasterizer.Rasterize(Dst, DstRect, Src);
    finally
      Transformer.Free;
    end;
  end;
  Dst.Changed(DstRect);
end;

//------------------------------------------------------------------------------

procedure RasterizeTransformation(Vectormap: TVectormap;
  Transformation: TTransformation; DstRect: TRect;
  CombineMode: TVectorCombineMode = vcmAdd;
  CombineCallback: TVectorCombineEvent = nil);
var
  I, J: Integer;
  P, Q, Progression: TFixedVector;
  ProgressionX, ProgressionY: TFixed;
  MapPtr: PFixedPointArray;
begin
  GR32.IntersectRect(DstRect, VectorMap.BoundsRect, DstRect);
  if GR32.IsRectEmpty(DstRect) then
    Exit;

  if not TTransformationAccess(Transformation).TransformValid then
    TTransformationAccess(Transformation).PrepareTransform;

  case CombineMode of
    vcmAdd:
      begin
        with DstRect do
        for I := Top to Bottom - 1 do
        begin
          MapPtr := @VectorMap.Vectors[I * VectorMap.Width];
          for J := Left to Right - 1 do
          begin
            P := FixedPoint(Integer(J - Left), Integer(I - Top));
            Q := Transformation.ReverseTransform(P);
            Inc(MapPtr[J].X, Q.X - P.X);
            Inc(MapPtr[J].Y, Q.Y - P.Y);
          end;
        end;
      end;
    vcmReplace:
      begin
        with DstRect do
        for I := Top to Bottom - 1 do
        begin
          MapPtr := @VectorMap.Vectors[I * VectorMap.Width];
          for J := Left to Right - 1 do
          begin
            P := FixedPoint(Integer(J - Left), Integer(I - Top));
            Q := Transformation.ReverseTransform(P);
            MapPtr[J].X := Q.X - P.X;
            MapPtr[J].Y := Q.Y - P.Y;
          end;
        end;
      end;
  else // vcmCustom
    ProgressionX := Fixed(1 / (DstRect.Right - DstRect.Left - 1));
    ProgressionY := Fixed(1 / (DstRect.Bottom - DstRect.Top - 1));
    Progression.Y := 0;
    with DstRect do for I := Top to Bottom - 1 do
    begin
      Progression.X := 0;
      MapPtr := @VectorMap.Vectors[I * VectorMap.Width];
      for J := Left to Right - 1 do
      begin
        P := FixedPoint(Integer(J - Left), Integer(I - Top));
        Q := Transformation.ReverseTransform(P);
        Q.X := Q.X - P.X;
        Q.Y := Q.Y - P.Y;
        CombineCallback(Q, Progression, MapPtr[J]);

        Inc(Progression.X, ProgressionX);
      end;
     Inc(Progression.Y, ProgressionY);
    end;
  end;
end;

//------------------------------------------------------------------------------
// Matrix conversion routines
//------------------------------------------------------------------------------
function FixedMatrix(const FloatMatrix: TFloatMatrix): TFixedMatrix;
begin
  Result[0,0] := Round(FloatMatrix[0,0] * FixedOne);
  Result[0,1] := Round(FloatMatrix[0,1] * FixedOne);
  Result[0,2] := Round(FloatMatrix[0,2] * FixedOne);
  Result[1,0] := Round(FloatMatrix[1,0] * FixedOne);
  Result[1,1] := Round(FloatMatrix[1,1] * FixedOne);
  Result[1,2] := Round(FloatMatrix[1,2] * FixedOne);
  Result[2,0] := Round(FloatMatrix[2,0] * FixedOne);
  Result[2,1] := Round(FloatMatrix[2,1] * FixedOne);
  Result[2,2] := Round(FloatMatrix[2,2] * FixedOne);
end;

function FloatMatrix(const FixedMatrix: TFixedMatrix): TFloatMatrix;
begin
  Result[0,0] := FixedMatrix[0,0] * FixedToFloat;
  Result[0,1] := FixedMatrix[0,1] * FixedToFloat;
  Result[0,2] := FixedMatrix[0,2] * FixedToFloat;
  Result[1,0] := FixedMatrix[1,0] * FixedToFloat;
  Result[1,1] := FixedMatrix[1,1] * FixedToFloat;
  Result[1,2] := FixedMatrix[1,2] * FixedToFloat;
  Result[2,0] := FixedMatrix[2,0] * FixedToFloat;
  Result[2,1] := FixedMatrix[2,1] * FixedToFloat;
  Result[2,2] := FixedMatrix[2,2] * FixedToFloat;
end;


//------------------------------------------------------------------------------
//
//      TTransformation
//
//------------------------------------------------------------------------------
function TTransformation.GetTransformedBounds: TFloatRect;
begin
  Result := GetTransformedBounds(FSrcRect);
end;

procedure TTransformation.Changed;
begin
  TransformValid := False;
  inherited;
end;

constructor TTransformation.Create;
begin
  // virtual constructor to be overriden in derived classes
end;

function TTransformation.GetTransformedBounds(const ASrcRect: TFloatRect): TFloatRect;
begin
  Result := ASrcRect;
end;

function TTransformation.HasTransformedBounds: Boolean;
begin
  Result := True;
end;

procedure TTransformation.PrepareTransform;
begin
  // Dummy
end;

function TTransformation.ReverseTransform(const P: TFloatPoint): TFloatPoint;
begin
  if not TransformValid then
    PrepareTransform;
  ReverseTransformFloat(P.X, P.Y, Result.X, Result.Y);
end;

function TTransformation.ReverseTransform(const P: TFixedPoint): TFixedPoint;
begin
  if not TransformValid then
    PrepareTransform;
  ReverseTransformFixed(P.X, P.Y, Result.X, Result.Y);
end;

function TTransformation.ReverseTransform(const P: TPoint): TPoint;
begin
  if not TransformValid then
    PrepareTransform;
  ReverseTransformInt(P.X, P.Y, Result.X, Result.Y);
end;

procedure TTransformation.ReverseTransformFixed(DstX, DstY: TFixed; out SrcX, SrcY: TFixed);
var
  X, Y: TFloat;
begin
  ReverseTransformFloat(DstX * FixedToFloat, DstY * FixedToFloat, X, Y);
  SrcX := Fixed(X);
  SrcY := Fixed(Y);
end;

procedure TTransformation.ReverseTransformFloat(DstX, DstY: TFloat;
  out SrcX, SrcY: TFloat);
begin
  // ReverseTransformFloat is the top precisionlevel, all descendants must override at least this level!
  raise ETransformNotImplemented.CreateFmt(RCStrReverseTransformationNotImplemented, [Self.Classname]);
end;

procedure TTransformation.ReverseTransformInt(DstX, DstY: Integer;
  out SrcX, SrcY: Integer);
var
  X, Y: TFixed;
begin
  ReverseTransformFixed(DstX shl 16, DstY shl 16, X, Y);
  SrcX := FixedRound(X);
  SrcY := FixedRound(Y);
end;

procedure TTransformation.SetSrcRect(const Value: TFloatRect);
begin
  FSrcRect := Value;
  Changed;
end;

function TTransformation.Transform(const P: TFloatPoint): TFloatPoint;
begin
  if not TransformValid then
    PrepareTransform;
  TransformFloat(P.X, P.Y, Result.X, Result.Y);
end;

function TTransformation.Transform(const P: TFixedPoint): TFixedPoint;
begin
  if not TransformValid then
    PrepareTransform;
  TransformFixed(P.X, P.Y, Result.X, Result.Y);
end;

function TTransformation.Transform(const P: TPoint): TPoint;
begin
  if not TransformValid then
    PrepareTransform;
  TransformInt(P.X, P.Y, Result.X, Result.Y);
end;

procedure TTransformation.TransformFixed(SrcX, SrcY: TFixed; out DstX,
  DstY: TFixed);
var
  X, Y: TFloat;
begin
  TransformFloat(SrcX * FixedToFloat, SrcY * FixedToFloat, X, Y);
  DstX := Fixed(X);
  DstY := Fixed(Y);
end;

procedure TTransformation.TransformFloat(SrcX, SrcY: TFloat; out DstX, DstY: TFloat);
begin
  // TransformFloat is the top precisionlevel, all descendants must override at least this level!
  raise ETransformNotImplemented.CreateFmt(RCStrForwardTransformationNotImplemented, [Self.Classname]);
end;

procedure TTransformation.TransformInt(SrcX, SrcY: Integer; out DstX, DstY: Integer);
var
  X, Y: TFixed;
begin
  TransformFixed(SrcX shl 16, SrcY shl 16, X, Y);
  DstX := FixedRound(X);
  DstY := FixedRound(Y);
end;


//------------------------------------------------------------------------------
//
//      TNestedTransformation
//
//------------------------------------------------------------------------------
constructor TNestedTransformation.Create;
begin
  FItems := TList.Create;
end;

destructor TNestedTransformation.Destroy;
begin
  if Assigned(FItems) then
    Clear;
  FItems.Free;
  inherited;
end;

function TNestedTransformation.Add(
  ItemClass: TTransformationClass): TTransformation;
begin
  Result := ItemClass.Create;
  FItems.Add(Result);
end;

procedure TNestedTransformation.Clear;
begin
  BeginUpdate;
  try
    while FItems.Count > 0 do
      Delete(FItems.Count - 1);
  finally
    EndUpdate;
  end;
end;

procedure TNestedTransformation.Delete(Index: Integer);
begin
  TTransformation(FItems[Index]).Free;
  FItems.Delete(Index);
end;

function TNestedTransformation.GetCount: Integer;
begin
  Result := FItems.Count;
end;

function TNestedTransformation.GetItem(Index: Integer): TTransformation;
begin
  Result := FItems[Index];
end;

function TNestedTransformation.Insert(Index: Integer;
  ItemClass: TTransformationClass): TTransformation;
begin
  BeginUpdate;
  try
    Result := Add(ItemClass);
  finally
    EndUpdate;
  end;
end;

procedure TNestedTransformation.PrepareTransform;
var
  Index: Integer;
begin
  for Index := 0 to Count - 1 do
    TTransformation(FItems[Index]).PrepareTransform;
end;

procedure TNestedTransformation.ReverseTransformFixed(DstX, DstY: TFixed;
  out SrcX, SrcY: TFixed);
var
  Index: Integer;
begin
  for Index := 0 to Count - 1 do
  begin
    TTransformation(FItems[Index]).ReverseTransformFixed(DstX, DstY, SrcX,
      SrcY);
    DstX := SrcX;
    DstY := SrcY;
  end;
end;

procedure TNestedTransformation.ReverseTransformFloat(DstX, DstY: TFloat;
  out SrcX, SrcY: TFloat);
var
  Index: Integer;
begin
  for Index := 0 to Count - 1 do
  begin
    TTransformation(FItems[Index]).ReverseTransformFloat(DstX, DstY, SrcX,
      SrcY);
    DstX := SrcX;
    DstY := SrcY;
  end;
end;

procedure TNestedTransformation.SetItem(Index: Integer;
  const Value: TTransformation);
begin
  TCollectionItem(FItems[Index]).Assign(Value);
end;

procedure TNestedTransformation.TransformFixed(SrcX, SrcY: TFixed; out DstX,
  DstY: TFixed);
var
  Index: Integer;
begin
  for Index := 0 to Count - 1 do
  begin
    TTransformation(FItems[Index]).TransformFixed(SrcX, SrcY, DstX, DstY);
    SrcX := DstX;
    SrcY := DstY;
  end;
end;

procedure TNestedTransformation.TransformFloat(SrcX, SrcY: TFloat; out DstX,
  DstY: TFloat);
var
  Index: Integer;
begin
  for Index := 0 to Count - 1 do
  begin
    TTransformation(FItems[Index]).TransformFloat(SrcX, SrcY, DstX, DstY);
    SrcX := DstX;
    SrcY := DstY;
  end;
end;


//------------------------------------------------------------------------------
//
//      T3x3Transformation
//
//------------------------------------------------------------------------------
procedure T3x3Transformation.PrepareTransform;
begin
  FInverseMatrix := Matrix;
  Invert(FInverseMatrix);

  // calculate a fixed point (65536) factors
  FInverseFixedMatrix := FixedMatrix(FInverseMatrix);
  FFixedMatrix := FixedMatrix(Matrix);

  TransformValid := True;
end;

procedure T3x3Transformation.ReverseTransformFixed(DstX, DstY: TFixed; out SrcX,
  SrcY: TFixed);
begin
  SrcX := FixedMul(DstX, FInverseFixedMatrix[0, 0]) +
    FixedMul(DstY, FInverseFixedMatrix[1, 0]) + FInverseFixedMatrix[2, 0];
  SrcY := FixedMul(DstX, FInverseFixedMatrix[0, 1]) +
    FixedMul(DstY, FInverseFixedMatrix[1, 1]) + FInverseFixedMatrix[2, 1];
end;

procedure T3x3Transformation.ReverseTransformFloat(DstX, DstY: TFloat; out SrcX,
  SrcY: TFloat);
begin
  SrcX := DstX * FInverseMatrix[0, 0] + DstY * FInverseMatrix[1, 0] +
    FInverseMatrix[2, 0];
  SrcY := DstX * FInverseMatrix[0, 1] + DstY * FInverseMatrix[1, 1] +
    FInverseMatrix[2, 1];
end;

procedure T3x3Transformation.TransformFixed(SrcX, SrcY: TFixed; out DstX,
  DstY: TFixed);
begin
  DstX := FixedMul(SrcX, FFixedMatrix[0, 0]) +
    FixedMul(SrcY, FFixedMatrix[1, 0]) + FFixedMatrix[2, 0];
  DstY := FixedMul(SrcX, FFixedMatrix[0, 1]) +
    FixedMul(SrcY, FFixedMatrix[1, 1]) + FFixedMatrix[2, 1];
end;

procedure T3x3Transformation.TransformFloat(SrcX, SrcY: TFloat; out DstX,
  DstY: TFloat);
begin
  DstX := SrcX * Matrix[0, 0] + SrcY * Matrix[1, 0] + Matrix[2, 0];
  DstY := SrcX * Matrix[0, 1] + SrcY * Matrix[1, 1] + Matrix[2, 1];
end;


//------------------------------------------------------------------------------
//
//      TAffineTransformation
//
//------------------------------------------------------------------------------
constructor TAffineTransformation.Create;
begin
  FStackLevel := 0;
  FStack := nil;
  Clear;
end;

procedure TAffineTransformation.Clear;
begin
  FMatrix := IdentityMatrix;
  Changed;
end;

procedure TAffineTransformation.Clear(BaseMatrix: TFloatMatrix);
begin
  FMatrix := BaseMatrix;
  Changed;
end;

function TAffineTransformation.GetTransformedBounds(const ASrcRect: TFloatRect): TFloatRect;
var
  V1, V2, V3, V4: TVector3f;
begin
  V1[0] := ASrcRect.Left;  V1[1] := ASrcRect.Top;    V1[2] := 1;
  V2[0] := ASrcRect.Right; V2[1] := V1[1];           V2[2] := 1;
  V3[0] := V1[0];          V3[1] := ASrcRect.Bottom; V3[2] := 1;
  V4[0] := V2[0];          V4[1] := V3[1];           V4[2] := 1;
  V1 := VectorTransform(Matrix, V1);
  V2 := VectorTransform(Matrix, V2);
  V3 := VectorTransform(Matrix, V3);
  V4 := VectorTransform(Matrix, V4);
  Result.Left   := Min(Min(V1[0], V2[0]), Min(V3[0], V4[0]));
  Result.Right  := Max(Max(V1[0], V2[0]), Max(V3[0], V4[0]));
  Result.Top    := Min(Min(V1[1], V2[1]), Min(V3[1], V4[1]));
  Result.Bottom := Max(Max(V1[1], V2[1]), Max(V3[1], V4[1]));
end;

procedure TAffineTransformation.Push;
begin
  Inc(FStackLevel);
  ReallocMem(FStack, FStackLevel * SizeOf(TFloatMatrix));
  Move(FMatrix, FStack^[FStackLevel - 1], SizeOf(TFloatMatrix));
end;

procedure TAffineTransformation.Pop;
begin
  if FStackLevel <= 0 then
    raise Exception.Create(RStrStackEmpty);

  Move(FStack^[FStackLevel - 1], FMatrix, SizeOf(TFloatMatrix));
  Dec(FStackLevel);
  Changed;
end;

procedure TAffineTransformation.Rotate(Alpha: TFloat);
var
  S, C: TFloat;
  M: TFloatMatrix;
begin
  Alpha := DegToRad(Alpha);
  GR32_Math.SinCos(Alpha, S, C);
  M := IdentityMatrix;
  M[0, 0] := C;   M[1, 0] := S;
  M[0, 1] := -S;  M[1, 1] := C;
  FMatrix := Mult(M, Matrix);
  Changed;
end;

procedure TAffineTransformation.Rotate(Cx, Cy, Alpha: TFloat);
var
  S, C: TFloat;
  M: TFloatMatrix;
begin
  if (Cx <> 0) or (Cy <> 0) then
    Translate(-Cx, -Cy);
  Alpha := DegToRad(Alpha);
  GR32_Math.SinCos(Alpha, S, C);
  M := IdentityMatrix;
  M[0, 0] := C;   M[1, 0] := S;
  M[0, 1] := -S;  M[1, 1] := C;
  FMatrix := Mult(M, Matrix);
  if (Cx <> 0) or (Cy <> 0) then
    Translate(Cx, Cy);
  Changed;
end;

procedure TAffineTransformation.Scale(Sx, Sy: TFloat);
var
  M: TFloatMatrix;
begin
  M := IdentityMatrix;
  M[0, 0] := Sx;
  M[1, 1] := Sy;
  FMatrix := Mult(M, Matrix);
  Changed;
end;

procedure TAffineTransformation.Scale(Value: TFloat);
var
  M: TFloatMatrix;
begin
  M := IdentityMatrix;
  M[0, 0] := Value;
  M[1, 1] := Value;
  FMatrix := Mult(M, Matrix);
  Changed;
end;

procedure TAffineTransformation.Skew(Fx, Fy: TFloat);
var
  M: TFloatMatrix;
begin
  M := IdentityMatrix;
  M[1, 0] := Fx;
  M[0, 1] := Fy;
  FMatrix := Mult(M, Matrix);
  Changed;
end;

procedure TAffineTransformation.Translate(Dx, Dy: TFloat);
var
  M: TFloatMatrix;
begin
  M := IdentityMatrix;
  M[2, 0] := Dx;
  M[2, 1] := Dy;
  FMatrix := Mult(M, Matrix);

  Changed;
end;


//------------------------------------------------------------------------------
//
//      TProjectiveTransformation
//
//------------------------------------------------------------------------------
function TProjectiveTransformation.GetTransformedBounds(const ASrcRect: TFloatRect): TFloatRect;
begin
  Result.Left   := Min(Min(FQuadX[0], FQuadX[1]), Min(FQuadX[2], FQuadX[3]));
  Result.Right  := Max(Max(FQuadX[0], FQuadX[1]), Max(FQuadX[2], FQuadX[3]));
  Result.Top    := Min(Min(FQuadY[0], FQuadY[1]), Min(FQuadY[2], FQuadY[3]));
  Result.Bottom := Max(Max(FQuadY[0], FQuadY[1]), Max(FQuadY[2], FQuadY[3]));
end;

function TProjectiveTransformation.GetX(Index: Integer): TFloat;
begin
  Result := FQuadX[Index];
end;

function TProjectiveTransformation.GetY(Index: Integer): TFloat;
begin
  Result := FQuadY[Index];
end;

procedure TProjectiveTransformation.PrepareTransform;
var
  dx1, dx2, px, dy1, dy2, py: TFloat;
  g, h, k: TFloat;
  R: TFloatMatrix;
begin
  px  := FQuadX[0] - FQuadX[1] + FQuadX[2] - FQuadX[3];
  py  := FQuadY[0] - FQuadY[1] + FQuadY[2] - FQuadY[3];

  if (px = 0) and (py = 0) then
  begin
    // affine mapping
    FMatrix[0, 0] := FQuadX[1] - FQuadX[0];
    FMatrix[1, 0] := FQuadX[2] - FQuadX[1];
    FMatrix[2, 0] := FQuadX[0];

    FMatrix[0, 1] := FQuadY[1] - FQuadY[0];
    FMatrix[1, 1] := FQuadY[2] - FQuadY[1];
    FMatrix[2, 1] := FQuadY[0];

    FMatrix[0, 2] := 0;
    FMatrix[1, 2] := 0;
    FMatrix[2, 2] := 1;
  end
  else
  begin
    // projective mapping
    dx1 := FQuadX[1] - FQuadX[2];
    dx2 := FQuadX[3] - FQuadX[2];
    dy1 := FQuadY[1] - FQuadY[2];
    dy2 := FQuadY[3] - FQuadY[2];
    k := dx1 * dy2 - dx2 * dy1;
    if k <> 0 then
    begin
      k := 1 / k;
      g := (px * dy2 - py * dx2) * k;
      h := (dx1 * py - dy1 * px) * k;

      FMatrix[0, 0] := FQuadX[1] - FQuadX[0] + g * FQuadX[1];
      FMatrix[1, 0] := FQuadX[3] - FQuadX[0] + h * FQuadX[3];
      FMatrix[2, 0] := FQuadX[0];

      FMatrix[0, 1] := FQuadY[1] - FQuadY[0] + g * FQuadY[1];
      FMatrix[1, 1] := FQuadY[3] - FQuadY[0] + h * FQuadY[3];
      FMatrix[2, 1] := FQuadY[0];

      FMatrix[0, 2] := g;
      FMatrix[1, 2] := h;
      FMatrix[2, 2] := 1;
    end
    else
    begin
      FillChar(FMatrix, SizeOf(FMatrix), 0);
    end;
  end;

  // denormalize texture space (u, v)
  R := IdentityMatrix;
  R[0, 0] := 1 / (SrcRect.Right - SrcRect.Left);
  R[1, 1] := 1 / (SrcRect.Bottom - SrcRect.Top);
  FMatrix := Mult(FMatrix, R);

  R := IdentityMatrix;
  R[2, 0] := -SrcRect.Left;
  R[2, 1] := -SrcRect.Top;
  FMatrix := Mult(FMatrix, R);

  inherited;
end;

procedure TProjectiveTransformation.SetX(Index: Integer; const Value: TFloat);
begin
  FQuadX[Index] := Value;
  Changed;
end;

procedure TProjectiveTransformation.SetY(Index: Integer; const Value: TFloat);
begin
  FQuadY[Index] := Value;
  Changed;
end;

procedure TProjectiveTransformation.ReverseTransformFixed(DstX, DstY: TFixed;
  out SrcX, SrcY: TFixed);
var
  Z: TFixed;
  Zf: TFloat;
begin
  Z := FixedMul(FInverseFixedMatrix[0, 2], DstX) +
    FixedMul(FInverseFixedMatrix[1, 2], DstY) + FInverseFixedMatrix[2, 2];

  if Z = 0 then
    Exit;

{$IFDEF UseInlining}
  SrcX := FixedMul(DstX, FInverseFixedMatrix[0, 0]) +
    FixedMul(DstY, FInverseFixedMatrix[1, 0]) + FInverseFixedMatrix[2, 0];
  SrcY := FixedMul(DstX, FInverseFixedMatrix[0,1]) +
    FixedMul(DstY, FInverseFixedMatrix[1, 1]) + FInverseFixedMatrix[2, 1];
{$ELSE}
  inherited;
{$ENDIF}

  if Z <> FixedOne then
  begin
    Zf := FixedOne / Z;
    SrcX := Round(SrcX * Zf);
    SrcY := Round(SrcY * Zf);
  end;
end;

procedure TProjectiveTransformation.ReverseTransformFloat(
  DstX, DstY: TFloat;
  out SrcX, SrcY: TFloat);
var
  Z: TFloat;
begin
  Z := FInverseMatrix[0, 2] * DstX + FInverseMatrix[1, 2] * DstY +
    FInverseMatrix[2, 2];

  if Z = 0 then
    Exit;

{$IFDEF UseInlining}
  SrcX := DstX * FInverseMatrix[0, 0] + DstY * FInverseMatrix[1, 0] +
    FInverseMatrix[2, 0];
  SrcY := DstX * FInverseMatrix[0, 1] + DstY * FInverseMatrix[1, 1] +
    FInverseMatrix[2, 1];
{$ELSE}
  inherited;
{$ENDIF}

  if Z <> 1 then
  begin
    Z := 1 / Z;
    SrcX := SrcX * Z;
    SrcY := SrcY * Z;
  end;
end;

procedure TProjectiveTransformation.TransformFixed(SrcX, SrcY: TFixed;
  out DstX, DstY: TFixed);
var
  Z: TFixed;
  Zf: TFloat;
begin
  Z := FixedMul(FFixedMatrix[0, 2], SrcX) +
    FixedMul(FFixedMatrix[1, 2], SrcY) + FFixedMatrix[2, 2];

  if Z = 0 then
    Exit;

{$IFDEF UseInlining}
  DstX := FixedMul(SrcX, FFixedMatrix[0, 0]) +
    FixedMul(SrcY, FFixedMatrix[1, 0]) + FFixedMatrix[2, 0];
  DstY := FixedMul(SrcX, FFixedMatrix[0, 1]) +
    FixedMul(SrcY, FFixedMatrix[1, 1]) + FFixedMatrix[2, 1];
{$ELSE}
  inherited;
{$ENDIF}

  if Z <> FixedOne then
  begin
    Zf := FixedOne / Z;
    DstX := Round(DstX * Zf);
    DstY := Round(DstY * Zf);
  end;
end;

procedure TProjectiveTransformation.TransformFloat(SrcX, SrcY: TFloat;
  out DstX, DstY: TFloat);
var
  Z: TFloat;
begin
  Z := FMatrix[0, 2] * SrcX + FMatrix[1, 2] * SrcY + FMatrix[2, 2];

  if Z = 0 then Exit;

{$IFDEF UseInlining}
  DstX := SrcX * Matrix[0, 0] + SrcY * Matrix[1, 0] + Matrix[2, 0];
  DstY := SrcX * Matrix[0, 1] + SrcY * Matrix[1, 1] + Matrix[2, 1];
{$ELSE}
  inherited;
{$ENDIF}

  if Z <> 1 then
  begin
    Z := 1 / Z;
    DstX := DstX * Z;
    DstY := DstY * Z;
  end;
end;


//------------------------------------------------------------------------------
//
//      TProjectiveTransformationEx
//
//------------------------------------------------------------------------------
// Based on amBitmapEditorToolForwardProjectiveTransform by Anders Melander
//------------------------------------------------------------------------------
{$ifndef FPC}
function TProjectiveTransformationEx.GetTransformedBounds(const ASrcRect: TFloatRect): TFloatRect;
var
  i: integer;
  Bounds: TFloatQuadrilateral;
begin
  if (FExtrapolate) then
    Exit(ASrcRect);

  // Transform the coords of the source rect to find the coords of
  // the corresponding target quad. Then return the boinding box of
  // this quad.
  for i := 0 to High(Bounds) do
    ReverseTransformFloat(FSourceQuad[i].X, FSourceQuad[i].Y, Bounds[i].X, Bounds[i].Y);

  Result.Left   := Min(Min(Bounds[0].X, Bounds[1].X), Min(Bounds[2].X, Bounds[3].X));
  Result.Right  := Max(Max(Bounds[0].X, Bounds[1].X), Max(Bounds[2].X, Bounds[3].X));
  Result.Top    := Min(Min(Bounds[0].Y, Bounds[1].Y), Min(Bounds[2].Y, Bounds[3].Y));
  Result.Bottom := Max(Max(Bounds[0].Y, Bounds[1].Y), Max(Bounds[2].Y, Bounds[3].Y));

(* Naive; Does not take projection to DestQuad into account.
  Result.Left   := Min(Min(FDestQuad[0].X, FDestQuad[1].X), Min(FDestQuad[2].X, FDestQuad[3].X));
  Result.Right  := Max(Max(FDestQuad[0].X, FDestQuad[1].X), Max(FDestQuad[2].X, FDestQuad[3].X));
  Result.Top    := Min(Min(FDestQuad[0].Y, FDestQuad[1].Y), Min(FDestQuad[2].Y, FDestQuad[3].Y));
  Result.Bottom := Max(Max(FDestQuad[0].Y, FDestQuad[1].Y), Max(FDestQuad[2].Y, FDestQuad[3].Y));
*)
end;

function TProjectiveTransformationEx.GetSource(Index: Integer): TFloatPoint;
begin
  Result := FSourceQuad[Index];
end;

function TProjectiveTransformationEx.GetSourceX(Index: Integer): TFloat;
begin
  Result := FSourceQuad[Index].X;
end;

function TProjectiveTransformationEx.GetSourceY(Index: Integer): TFloat;
begin
  Result := FSourceQuad[Index].Y;
end;

function TProjectiveTransformationEx.GetDest(Index: Integer): TFloatPoint;
begin
  Result := FDestQuad[Index];
end;

function TProjectiveTransformationEx.GetDestX(Index: Integer): TFloat;
begin
  Result := FDestQuad[Index].X;
end;

function TProjectiveTransformationEx.GetDestY(Index: Integer): TFloat;
begin
  Result := FDestQuad[Index].Y;
end;

procedure TProjectiveTransformationEx.SetSource(Index: Integer; const Value: TFloatPoint);
begin
  FSourceQuad[Index] := Value;
  Changed;
end;

procedure TProjectiveTransformationEx.SetSourceQuad(const Value: TFloatQuadrilateral);
begin
  FSourceQuad := Value;
  Changed;
end;

procedure TProjectiveTransformationEx.SetSourceX(Index: Integer; const Value: TFloat);
begin
  FSourceQuad[Index].Y := Value;
  Changed;
end;

procedure TProjectiveTransformationEx.SetSourceY(Index: Integer; const Value: TFloat);
begin
  FSourceQuad[Index].X := Value;
  Changed;
end;

procedure TProjectiveTransformationEx.SetDest(Index: Integer; const Value: TFloatPoint);
begin
  FDestQuad[Index] := Value;
  Changed;
end;

procedure TProjectiveTransformationEx.SetDestQuad(const Value: TFloatQuadrilateral);
begin
  FDestQuad := Value;
  Changed;
end;

procedure TProjectiveTransformationEx.SetDestX(Index: Integer; const Value: TFloat);
begin
  FDestQuad[Index].X := Value;
  Changed;
end;

procedure TProjectiveTransformationEx.SetDestY(Index: Integer; const Value: TFloat);
begin
  FDestQuad[Index].Y := Value;
  Changed;
end;

procedure TProjectiveTransformationEx.PrepareTransform;
//------------------------------------------------------------------------------
// From "Fundamentals of Texture Mapping and Image Warping" by Paul S. Heckbert:
//------------------------------------------------------------------------------
//
// The general form of a projective mapping is a rational linear mapping:
//
//        au + bv + c            du + ev + f
//   x = -------------  ,   y = -------------                                                                   [1]
//        gu + hv + i            gu + hv + i
//
// Manipulation of projective mappings is much easier in the homogeneous matrix notation:
//
//   Pd = Ps * Msd
//
//        ┌    ┐   ┌         ┐┌    ┐
//        │ x' │   │ a  d  g ││ u' │
//      = │ y' │ = │ b  e  h ││ v' │
//        │ w  │   │ c  f  i ││ q  │
//        └    ┘   └         ┘└    ┘
//
//            T               T                     T              T
// where (x,y) = (x'/ w, y'/w) for w ≠ 0, and (u,v) = (u'/q, v'/q) for q ≠ 0.
//
// Although there are 9 coefficients in the matrix above, these mappings are homogeneous, so
// any nonzero scalar multiple of these matrices gives an equivalent mapping. Hence there are only
// 8 degrees of freedom in a 2-D projective mapping. We can assume without loss of generality that
// i=1 except in the special case that source point (0, 0)T maps to a point at infinity. A projective
// mapping is affine when g=h=0.
//
// ...
//
// Projective mappings may be composed by concatenating their matrices.
//
// Another remarkable property is that the inverse of a projective mapping is a projective mapping.
// This is intuitively explained by reversing the plane-to-plane mapping by which a projective mapping
// is defined. The matrix for the inverse mapping is the inverse or adjoint of the forward mapping. (The
// adjoint of a matrix is the transpose of the matrix of cofactors; M^-1 = adj(M)/det(M)).
// In homogeneous algebra, the adjoint matrix can be used in place of the inverse matrix whenever an
// inverse transform is needed, since the two are scalar multiples of each other, and the adjoint always
// exists, while the inverse does not if the matrix is singular. The inverse transformation is thus:
//
//   Ps = Msd * Pd
//
//        ┌    ┐   ┌         ┐┌    ┐
//        │ u' │   │ A  D  G ││ x' │
//      = │ v' │ = │ B  E  H ││ y' │
//        │ q  │   │ C  F  I ││ w  │
//        └    ┘   └         ┘└    ┘
//
//        ┌                     ┐┌    ┐
//        │ ei-fh  ch-bi  bf-ce ││ x' │
//      = │ fg-di  ai-cg  cd-af ││ y' │
//        │ dh-eg  bg-ah  ae-bd ││ w  │
//        └                     ┘└    ┘
//                                                                T           T
// When mapping a point by the inverse transform we compute (u, v) from (x, y). If w ≠ 0 and
// q ≠ 0 then we can choose w = 1 and calculate:
//
//        Ax + By + C            Dx + Ey + F
//   u = -------------  ,   y = -------------                                                                   [2]
//        Gx + Hy + I            Gx + Hy + I
//
// ...
//
// In an interactive image warper one might specify the four corners of source and destination quadrilaterals
// with a tablet or mouse, and wish to warp one area to the other. This sort of task is an ideal
// application of projective mappings, but how do we find the mapping matrix?
//
// A projective mapping has 8 degrees of freedom which can be determined from the source and
// destination coordinates of the four corners of a quadrilateral. Let the correspondence map (uk, vk)^T
// to (xk, yk)^T for vertices numbered cyclically k = 0,1,2,3. All coordinates are assumed to be real
// (finite). To compute the forward mapping matrix Msd, assuming that i= 1, we have eight equations
// in the eight unknowns a-h:
//
//         auk + bvk + c
//   xk = ---------------  =>  uka + vkb + c - ukxkg - vkxkh = xk
//         guk + hvk + 1
//
//         duk + evk + f
//   yk = ---------------  =>  ukd + vke + f - ukykg - vkykh = yk
//         guk + hvk + 1
//
// for k = 0,1,2,3. This can be rewritten as an 8 × 8 system:
//
//        ┌                                    ┐┌   ┐   ┌    ┐
//        │ u0  v0   1   0   0   0 -u0x0 -v0x0 ││ a │   │ x0 │
//        │ u1  v1   1   0   0   0 -u1x1 -v1x1 ││ b │   │ x1 │
//        │ u2  v2   1   0   0   0 -u2x2 -v2x2 ││ c │   │ x2 │
//        │ u3  v3   1   0   0   0 -u3x3 -v3x3 ││ d │ = │ x3 │
//        │  0   0   0  u0  v0   1 -u0y0 -v0y0 ││ e │   │ y0 │
//        │  0   0   0  u1  v1   1 -u1y1 -v1y1 ││ f │   │ y1 │
//        │  0   0   0  u2  v2   1 -u2y2 -v2y2 ││ g │   │ y2 │
//        │  0   0   0  u3  v3   1 -u3y3 -v3y3 ││ h │   │ y3 │
//        └                                    ┘└   ┘   └    ┘
//
// This linear system can be solved using Gaussian elimination or other methods for the forward mapping
// coefficients a-h. If the inverse mapping is desired instead, then either we compute the adjoint
// of Msd or we follow the same procedure, starting from equation [2] instead of [1], and solve an
// 8 × 8 system for coefficients A-H.
//
// In speed-critical special cases, there are more efficient formulas for computing the mapping
// matrix. The formula above handles the case where the polygon is a general quadrilateral in both
// source and destination spaces. We will consider three additional cases: square-to-quadrilateral,
// quadrilateral-to-square, and (again) the general quadrilateral-to-quadrilateral mapping.
//
// Case 1. The system is easily solved symbolically in the special case where the uv quadrilateral
// is a unit square. If the vertex correspondence is as follows:
//
//     x  y  u  v
//   -------------
//    x0 y0  0  0
//    x1 y1  1  0
//    x2 y2  1  1
//    x3 y3  0  1
//
// then the eight equations reduce to
//
//                       c = x0
//             a + c - gx1 = x1
//       a + c - gx2 - hx2 = x2
//             b + c - hx3 = x3
//                       f = y0
//             d + f - gy1 = y1
//   d + e + f - gy2 - hy2 = y2
//             e + f - hy3 = y3
//
// If we define
//
//   ∆x1 = x1 - x2,  ∆x2 = x3 - x2,  ∑x = x0 - x1 + x2 - x3                                                     [a.1]
//   ∆y1 = y1 - y2,  ∆y2 = y3 - y2,  ∑y = y0 - y1 + y2 - y3                                                     [a.2]
//
// then the solution splits into two sub-cases:
//
//   (a) ∑x = 0 and ∑y = 0. This implies that the xy polygon is a parallelogram, so the mapping is
//       affine, and a = x1 - x0,  b = x2 - x1,  c = x0,  d = y1 - y0,  e = y2 - y1,  f = y0,  g = 0,  h = 0.   [b]
//
//   (b) ∑x ≠ 0 or ∑y ≠ 0 gives a projective mapping:
//           │ ∑x  ∆x2 │   │ ∆x1 ∆x2 │
//       g = │         │ / │         │                                                                          [c.1]
//           │ ∑y  ∆y2 │   │ ∆y1 ∆y2 │
//
//
//           │ ∆x1  ∑x │   │ ∆x1 ∆x2 │
//       h = │         │ / │         │                                                                          [c.2]
//           │ ∆y1  ∑y │   │ ∆y1 ∆y2 │
//
//       a = x1 - x0 + gx1                                                                                      [c.3]
//       b = x3 - x0 + hx3
//       c = x0
//       d = y1 - y0 + gy1
//       e = y3 - y0 + hy3
//       f = y0
//
// This computation is much faster than a straightforward 8 × 8 system solver. The mapping above is
// easily generalized to map a rectangle to a quadrilateral by pre-multiplying with a scale and translation
// matrix.
//
// Case 2. The inverse mapping, a quadrilateral to a square, can also be optimized. It turns out
// that the most efficient algorithm for computing this is not purely symbolic, as in the previous case,
// but numerical. We use the square-to-quadrilateral formulas just described to find the inverse of the
// desired mapping, and then take its adjoint to compute the quadrilateral-to-square mapping.
//
// Case 3. Since we can compute quadrilateral-to-square and square-to-quadrilateral mappings
// quickly, the two mappings can easily be composed to yield a general quadrilateral-to-
// mapping. This solution method is faster than a general 8 × 8 system solver.
//
  procedure CreateProjectiveMapping(const Quad: TFloatQuadrilateral; var Matrix: TFloatMatrix);
  var
    ∑x, ∑y: TFloat;
    ∆x1, ∆x2, ∆y1, ∆y2: TFloat;
    g, h, k: TFloat;
  begin
    ∑x := Quad[0].X - Quad[1].X + Quad[2].X - Quad[3].X;                        // See [a]
    ∑y := Quad[0].Y - Quad[1].Y + Quad[2].Y - Quad[3].Y;

    if (IsZero(∑x)) and (IsZero(∑y)) then                                       // See [b]
    begin
      // Quadrilateral is a parallelogram - Mapping is affine
      //   ┌         ┐
      //   │ a  d  0 │
      //   │ b  e  0 │
      //   │ c  f  1 │
      //   └         ┘
      Matrix[0, 0] := Quad[1].X - Quad[0].X;    // a
      Matrix[1, 0] := Quad[2].X - Quad[1].X;    // b
      Matrix[2, 0] := Quad[0].X;                // c

      Matrix[0, 1] := Quad[1].Y - Quad[0].Y;    // d
      Matrix[1, 1] := Quad[2].Y - Quad[1].Y;    // e
      Matrix[2, 1] := Quad[0].Y;                // f

      Matrix[0, 2] := 0;                        // g
      Matrix[1, 2] := 0;                        // h
      Matrix[2, 2] := 1;                        // i
    end else
    begin
      // Projective mapping
      //   ┌         ┐
      //   │ a  d  g │
      //   │ b  e  h │
      //   │ c  f  1 │
      //   └         ┘
      ∆x1 := Quad[1].X - Quad[2].X;                                             // See [a]
      ∆x2 := Quad[3].X - Quad[2].X;
      ∆y1 := Quad[1].Y - Quad[2].Y;
      ∆y2 := Quad[3].Y - Quad[2].Y;

      k := ∆x1 * ∆y2 - ∆x2 * ∆y1;

      if (not IsZero(k)) then
      begin
        k := 1 / k; // Avoid (one) costly divisions below
        g := (∑x * ∆y2 - ∑y * ∆x2) * k;                                         // See [c]
        h := (∆x1 * ∑y - ∆y1 * ∑x) * k;

        Matrix[0, 0] := Quad[1].X - Quad[0].X + g * Quad[1].X;  // a
        Matrix[1, 0] := Quad[3].X - Quad[0].X + h * Quad[3].X;  // b
        Matrix[2, 0] := Quad[0].X;                              // c

        Matrix[0, 1] := Quad[1].Y - Quad[0].Y + g * Quad[1].Y;  // d
        Matrix[1, 1] := Quad[3].Y - Quad[0].Y + h * Quad[3].Y;  // e
        Matrix[2, 1] := Quad[0].Y;                              // f

        Matrix[0, 2] := g;                                      // g
        Matrix[1, 2] := h;                                      // h
        Matrix[2, 2] := 1;                                      // i
      end else
        Matrix := Default(TFloatMatrix);
    end;
  end;

var
  SourceMatrix: TFloatMatrix;
  DestMatrix: TFloatMatrix;
//  R: TFloatMatrix;
begin
  CreateProjectiveMapping(FSourceQuad, SourceMatrix);

  CreateProjectiveMapping(FDestQuad, DestMatrix);
  Invert(DestMatrix);

  FMatrix := Mult(SourceMatrix, DestMatrix);

(*
  // Denormalize texture space (u, v)
  // Scale
  R := IdentityMatrix;
  R[0, 0] := 1 / (SrcRect.Right - SrcRect.Left);
  R[1, 1] := 1 / (SrcRect.Bottom - SrcRect.Top);
  FMatrix := Mult(FMatrix, R);

  // Translate
  R := IdentityMatrix;
  R[2, 0] := -SrcRect.Left;
  R[2, 1] := -SrcRect.Top;
  FMatrix := Mult(FMatrix, R);
*)

  inherited;
end;

procedure TProjectiveTransformationEx.ReverseTransformFixed(DstX, DstY: TFixed; out SrcX, SrcY: TFixed);
var
  Z: TFixed;
  Zf: TFloat;
begin
  Z := FixedMul(FInverseFixedMatrix[0, 2], DstX) +
       FixedMul(FInverseFixedMatrix[1, 2], DstY) +
                FInverseFixedMatrix[2, 2];

  if Z = 0 then
    Exit;

{$IFDEF UseInlining}
  SrcX := FixedMul(FInverseFixedMatrix[0, 0], DstX) +
          FixedMul(FInverseFixedMatrix[1, 0], DstY) +
                   FInverseFixedMatrix[2, 0];

  SrcY := FixedMul(FInverseFixedMatrix[0, 1], DstX) +
          FixedMul(FInverseFixedMatrix[1, 1], DstY) +
                   FInverseFixedMatrix[2, 1];
{$ELSE}
  inherited;
{$ENDIF}

  if Z <> FixedOne then
  begin
    Zf := FixedOne / Z;
    SrcX := Round(SrcX * Zf);
    SrcY := Round(SrcY * Zf);
  end;
end;

procedure TProjectiveTransformationEx.ReverseTransformFloat(DstX, DstY: TFloat; out SrcX, SrcY: TFloat);
var
  Z: TFloat;
begin
  Z := FInverseMatrix[0, 2] * DstX +
       FInverseMatrix[1, 2] * DstY +
       FInverseMatrix[2, 2];

  if IsZero(Z) then
    Exit;

{$IFDEF UseInlining}
  SrcX := FInverseMatrix[0, 0] * DstX +
          FInverseMatrix[1, 0] * DstY +
          FInverseMatrix[2, 0];

  SrcY := FInverseMatrix[0, 1] * DstX +
          FInverseMatrix[1, 1] * DstY +
          FInverseMatrix[2, 1];
{$ELSE}
  inherited;
{$ENDIF}

  if Z <> 1 then
  begin
    Z := 1 / Z;
    SrcX := SrcX * Z;
    SrcY := SrcY * Z;
  end;
end;

procedure TProjectiveTransformationEx.TransformFixed(SrcX, SrcY: TFixed; out DstX, DstY: TFixed);
var
  Z: TFixed;
  Zf: TFloat;
begin
  Z := FixedMul(FFixedMatrix[0, 2], SrcX) +
       FixedMul(FFixedMatrix[1, 2], SrcY) +
                FFixedMatrix[2, 2];

  if Z = 0 then
    Exit;

{$IFDEF UseInlining}
  DstX := FixedMul(FFixedMatrix[0, 0], SrcX) +
          FixedMul(FFixedMatrix[1, 0], SrcY) +
                   FFixedMatrix[2, 0];
  DstY := FixedMul(FFixedMatrix[0, 1], SrcX) +
          FixedMul(FFixedMatrix[1, 1], SrcY) +
                   FFixedMatrix[2, 1];
{$ELSE}
  inherited;
{$ENDIF}

  if Z <> FixedOne then
  begin
    Zf := FixedOne / Z;
    DstX := Round(DstX * Zf);
    DstY := Round(DstY * Zf);
  end;
end;

procedure TProjectiveTransformationEx.TransformFloat(SrcX, SrcY: TFloat; out DstX, DstY: TFloat);
var
  Z: TFloat;
begin
  Z := FMatrix[0, 2] * SrcX +
       FMatrix[1, 2] * SrcY +
       FMatrix[2, 2];

  if IsZero(Z) then
    Exit;

{$IFDEF UseInlining}
  DstX := FMatrix[0, 0] * SrcX +
          FMatrix[1, 0] * SrcY +
          FMatrix[2, 0];
  DstY := FMatrix[0, 1] * SrcX +
          FMatrix[1, 1] * SrcY +
          FMatrix[2, 1];
{$ELSE}
  inherited;
{$ENDIF}

  if Z <> 1 then
  begin
    Z := 1 / Z;
    DstX := DstX * Z;
    DstY := DstY * Z;
  end;
end;
{$endif FPC}


//------------------------------------------------------------------------------
//
//      TTwirlTransformation
//
//------------------------------------------------------------------------------
constructor TTwirlTransformation.Create;
begin
  FTwirl := 0.03;
end;

function TTwirlTransformation.GetTransformedBounds(const ASrcRect: TFloatRect): TFloatRect;
var
  Cx, Cy, R: TFloat;
const
  CPiHalf: TFloat = 0.5 * Pi;
begin
  Cx := (ASrcRect.Left + ASrcRect.Right) * 0.5;
  Cy := (ASrcRect.Top + ASrcRect.Bottom) * 0.5;
  R := Max(Cx - ASrcRect.Left, Cy - ASrcRect.Top);
  Result.Left := Cx - R * CPiHalf;
  Result.Right := Cx + R * CPiHalf;
  Result.Top := Cy - R * CPiHalf;
  Result.Bottom := Cy + R * CPiHalf;
end;

procedure TTwirlTransformation.PrepareTransform;
begin
  with FSrcRect do
  begin
    Frx := (Right - Left) * 0.5;
    Fry := (Bottom - Top) * 0.5;
  end;
  TransformValid := True;
end;

procedure TTwirlTransformation.ReverseTransformFloat(DstX, DstY: TFloat;
  out SrcX, SrcY: TFloat);
var
  xf, yf, r, t: Single;
begin
  xf := DstX - Frx;
  yf := DstY - Fry;

  r := GR32_Math.Hypot(xf, yf);
  t := ArcTan2(yf, xf) + r * FTwirl;
  GR32_Math.SinCos(t, yf, xf);

  SrcX := Frx + r * xf;
  SrcY := Fry + r * yf;
end;

procedure TTwirlTransformation.SetTwirl(const Value: TFloat);
begin
  FTwirl := Value;
  Changed;
end;


//------------------------------------------------------------------------------
//
//      TBloatTransformation
//
//------------------------------------------------------------------------------
constructor TBloatTransformation.Create;
begin
  FBloatPower := 0.3;
end;

procedure TBloatTransformation.PrepareTransform;
begin
  FPiW := (Pi / (FSrcRect.Right - FSrcRect.Left));
  FPiH := (Pi / (FSrcRect.Bottom - FSrcRect.Top));
  FBP := FBloatPower * Max(FSrcRect.Right - FSrcRect.Left, FSrcRect.Bottom - FSrcRect.Top);
  TransformValid := True;
end;

procedure TBloatTransformation.ReverseTransformFloat(DstX, DstY: TFloat;
  out SrcX, SrcY: TFloat);
var
  SinY, CosY, SinX, CosX, t: Single;
begin
  GR32_Math.SinCos(FPiH * DstY, SinY, CosY);
  GR32_Math.SinCos(FPiW * DstX, SinX, CosX);
  t := FBP * SinY * SinX;
  SrcX := DstX + t * CosX;
  SrcY := DstY + t * CosY;
end;

procedure TBloatTransformation.TransformFloat(DstX, DstY: TFloat;
  out SrcX, SrcY: TFloat);
var
  SinY, CosY, SinX, CosX, t: Single;
begin
  GR32_Math.SinCos(-FPiH * DstY, SinY, CosY);
  GR32_Math.SinCos(-FPiW * DstX, SinX, CosX);
  t := FBP * SinY * SinX;
  SrcX := DstX + t * CosX;
  SrcY := DstY + t * CosY;
end;

procedure TBloatTransformation.SetBloatPower(const Value: TFloat);
begin
  FBloatPower := Value;
  Changed;
end;


//------------------------------------------------------------------------------
//
//      TFishEyeTransformation
//
//------------------------------------------------------------------------------
procedure TFishEyeTransformation.PrepareTransform;
begin
  with FSrcRect do
  begin
    Frx := (Right - Left) * 0.5;
    Fry := (Bottom - Top) * 0.5;
    if Frx <= Fry then
    begin
      FMinR := Frx;
      Sx := 1;
      Sy:= Frx / Fry;
    end
    else
    begin
      FMinR := Fry;
      Sx:= Fry / Frx;
      Sy := 1;
    end;
    Fsr := 1 / FMinR;
    Faw := ArcSin(Constrain(FMinR * Fsr, -1, 1));
    if Faw <> 0 then
      Faw := 1 / Faw;
    Faw := Faw * FMinR
  end;
  TransformValid := True;  
end;

procedure TFishEyeTransformation.ReverseTransformFloat(DstX, DstY: TFloat;
  out SrcX, SrcY: TFloat);
var
  d, Xrx, Yry: TFloat;
begin
  Yry := (DstY - Fry) * sy;
  Xrx := (DstX - Frx) * sx;
  d := GR32_Math.Hypot(Xrx, Yry);
  if (d < FMinR) and (d > 0) then
  begin
    d := ArcSin(d * Fsr) * Faw / d;
    SrcX := Frx + Xrx * d;
    SrcY := Fry + Yry * d;
  end
  else
  begin
    SrcX := DstX;
    SrcY := DstY;
  end;
end;


//------------------------------------------------------------------------------
//
//      TPolarTransformation
//
//------------------------------------------------------------------------------
procedure TPolarTransformation.PrepareTransform;
begin
  Sx := SrcRect.Right - SrcRect.Left;
  Sy := SrcRect.Bottom - SrcRect.Top;
  Cx := (DstRect.Left + DstRect.Right) * 0.5;
  Cy := (DstRect.Top + DstRect.Bottom) * 0.5;
  Dx := DstRect.Right - Cx;
  Dy := DstRect.Bottom - Cy;

  Rt := (1 / (PI * 2)) * Sx;

  Rt2 := Sx;
  if Rt2 <> 0 then
    Rt2 := 1 / Sx
  else
    Rt2 := 0.00000001;
  Rt2 := Rt2 * 2 * Pi;

  Rr := Sy;
  if Rr <> 0 then
    Rr := 1 / Rr
  else
    Rr := 0.00000001;

  Rcx := Cx;
  if Rcx <> 0 then
    Rcx := 1 / Rcx
  else
    Rcx := 0.00000001;

  Rcy := Cy;
  if Rcy <> 0 then
    Rcy := 1 / Rcy
   else
    Rcy := 0.00000001;

  TransformValid := True;
end;

procedure TPolarTransformation.SetDstRect(const Value: TFloatRect);
begin
  FDstRect := Value;
  Changed;
end;

procedure TPolarTransformation.TransformFloat(SrcX, SrcY: TFloat; out DstX,
  DstY: TFloat);
var
  R, Theta, S, C: TFloat;
begin
  Theta := (SrcX - SrcRect.Left) * Rt2 + Phase;
  R := (SrcY - SrcRect.Bottom) * Rr;
  GR32_Math.SinCos(Theta, S, C);

  DstX := Dx * R * C + Cx;
  DstY := Dy * R * S + Cy;
end;

procedure TPolarTransformation.ReverseTransformFloat(DstX, DstY: TFloat;
  out SrcX, SrcY: TFloat);
const
  PI2 = 2 * PI;
var
  Dcx, Dcy, Theta: TFloat;
begin
  Dcx := (DstX - Cx) * Rcx;
  Dcy := (DstY - Cy) * Rcy;

  Theta := ArcTan2(Dcy, Dcx) + Pi - Phase;
  if Theta < 0 then
    Theta := Theta + PI2;

  SrcX := SrcRect.Left + Theta * Rt;
  SrcY := SrcRect.Bottom - GR32_Math.Hypot(Dcx, Dcy) * Sy;
end;


procedure TPolarTransformation.SetPhase(const Value: TFloat);
begin
  FPhase := Value;
  Changed;
end;


//------------------------------------------------------------------------------
//
//      TPathTransformation
//
//------------------------------------------------------------------------------
destructor TPathTransformation.Destroy;
begin
  FTopHypot := nil;
  FBottomHypot := nil;
  inherited;
end;

procedure TPathTransformation.PrepareTransform;
var
  I: Integer;
  L, DDist: TFloat;
begin
  if not (Assigned(FTopCurve) and Assigned(FBottomCurve)) then
    raise ETransformError.Create(RCStrTopBottomCurveNil);

  SetLength(FTopHypot, Length(FTopCurve));
  SetLength(FBottomHypot, Length(FBottomCurve));

  L := 0;
  for I := 0 to High(FTopCurve) - 1 do
  begin
    FTopHypot[I].Dist := L;
    with FTopCurve[I + 1] do
      L := L + GR32_Math.Hypot(FTopCurve[I].X - X, FTopCurve[I].Y - Y);
  end;
  FTopLength := L;

  for I := 1 to High(FTopCurve) do
    with FTopHypot[I] do
    begin
      DDist := Dist - FTopHypot[I - 1].Dist;
      if DDist <> 0 then
        RecDist := 1 / DDist
      else
      if I > 1 then
        RecDist := FTopHypot[I - 1].RecDist
      else
        RecDist := 0;
    end;

  L := 0;
  for I := 0 to High(FBottomCurve) - 1 do
  begin
    FBottomHypot[I].Dist := L;
    with FBottomCurve[I + 1] do
      L := L + GR32_Math.Hypot(FBottomCurve[I].X - X, FBottomCurve[I].Y - Y);
  end;
  FBottomLength := L;

  for I := 1 to High(FBottomCurve) do
    with FBottomHypot[I] do
    begin
      DDist := Dist - FBottomHypot[I - 1].Dist;
      if DDist <> 0 then
        RecDist := 1 / DDist
      else
      if I > 1 then
        RecDist := FBottomHypot[I - 1].RecDist
      else
        RecDist := 0;
    end;

  rdx := 1 / (SrcRect.Right - SrcRect.Left);
  rdy := 1 / (SrcRect.Bottom - SrcRect.Top);

  TransformValid := True;
end;

procedure TPathTransformation.SetBottomCurve(const Value: TArrayOfFloatPoint);
begin
  FBottomCurve := Value;
  Changed;
end;

procedure TPathTransformation.SetTopCurve(const Value: TArrayOfFloatPoint);
begin
  FTopCurve := Value;
  Changed;
end;

procedure TPathTransformation.TransformFloat(SrcX, SrcY: TFloat; out DstX,
  DstY: TFloat);
var
  I, H: Integer;
  X, Y, fx, dx, dy, r, Tx, Ty, Bx, By: TFloat;
begin
  X := (SrcX - SrcRect.Left) * rdx;
  Y := (SrcY - SrcRect.Top) * rdy;

  fx := X * FTopLength;
  I := 1;
  H := High(FTopHypot);
  while (FTopHypot[I].Dist < fx) and (I < H) do
    Inc(I);


  with FTopHypot[I] do
    r := (Dist - fx) * RecDist;

  dx := (FTopCurve[I - 1].X - FTopCurve[I].X);
  dy := (FTopCurve[I - 1].Y - FTopCurve[I].Y);
  Tx := FTopCurve[I].X + r * dx;
  Ty := FTopCurve[I].Y + r * dy;

  fx := X * FBottomLength;
  I := 1;
  H := High(FBottomHypot);
  while (FBottomHypot[I].Dist < fx) and (I < H) do
    Inc(I);

  with FBottomHypot[I] do
    r := (Dist - fx) * RecDist;

  dx := (FBottomCurve[I - 1].X - FBottomCurve[I].X);
  dy := (FBottomCurve[I - 1].Y - FBottomCurve[I].Y);
  Bx := FBottomCurve[I].X + r * dx;
  By := FBottomCurve[I].Y + r * dy;

  DstX := Tx + Y * (Bx - Tx);
  DstY := Ty + Y * (By - Ty);
end;


//------------------------------------------------------------------------------
//
//      TDisturbanceTransformation
//
//------------------------------------------------------------------------------
function TDisturbanceTransformation.GetTransformedBounds(
  const ASrcRect: TFloatRect): TFloatRect;
begin
  Result := ASrcRect;
  GR32.InflateRect(Result, 0.5 * FDisturbance, 0.5 * FDisturbance);
end;

procedure TDisturbanceTransformation.ReverseTransformFloat(DstX,
  DstY: TFloat; out SrcX, SrcY: TFloat);
begin
  SrcX := DstX + (Random - 0.5) * FDisturbance;
  SrcY := DstY + (Random - 0.5) * FDisturbance;
end;

procedure TDisturbanceTransformation.SetDisturbance(const Value: TFloat);
begin
  FDisturbance := Value;
  Changed;  
end;


//------------------------------------------------------------------------------
//
//      TRadialDistortionTransformation
//
//------------------------------------------------------------------------------
constructor TRadialDistortionTransformation.Create;
begin
  FCoefficient1 := 0;
  FCoefficient2 := 0;
  FScale := 1;
  FMapElements := 0;
end;

function TRadialDistortionTransformation.HasTransformedBounds: Boolean;
begin
  Result := False;
end;

procedure TRadialDistortionTransformation.PrepareReverseMap;
var
  i, j, LowerI, UpperI, jmax: Integer;
  r_src, r_tgt, LowerValue, UpperValue: TFloat;
{$if defined(MSWINDOWS) and defined(DEBUG)}
  // some counters to evaluate the mapping
  interpolated, unset, mapToSameIndex, IndexOutOfRange: Integer;
{$ifend}
begin
  if MapElements <= 1 then
    MapElements := Trunc(r_0);

  r_tgt_max := 2;
  r_tgt_min := -0.5;

  SetLength(Map, MapElements);
  for i := 0 to High(Map) do
    Map[i] := -1;

  jmax := 1000;
{$if defined(MSWINDOWS) and defined(DEBUG)}
  mapToSameIndex := 0;
  IndexOutOfRange := 0;
{$ifend}
  for j := 0 to jmax do
  begin
    r_src := j/jmax*2;
    r_tgt := Scale*(1 + FCoefficient1 * Sqr(r_src) + FCoefficient2 * Power(r_src, 4));
    Assert(InRange(r_tgt, r_tgt_min, r_tgt_max));
    i := Trunc((r_tgt*r_src-r_tgt_min)/(r_tgt_max-r_tgt_min)*(High(Map)-1));
    if not InRange(i, 0, High(Map)) then
    begin
{$if defined(MSWINDOWS) and defined(DEBUG)}
      Inc(IndexOutOfRange);
      // OutputDebugString(PChar(Format('PrepareReverseMap: i=%d out of range (0, MapElements=%d), r_tgt=%f', [ i, MapElements, r_tgt ])))
{$ifend}
    end
    else
    if Map[i]<>-1 then
    begin
{$if defined(MSWINDOWS) and defined(DEBUG)}
      Inc(mapToSameIndex);
      // OutputDebugString(PChar(Format('PrepareReverseMap: Map[i=%d] already has value %f (wanted to put %f there)', [ i, Map[i], r_tgt ])))
{$ifend}
    end
    else
      Map[i] := r_tgt;
  end;

{$if defined(MSWINDOWS) and defined(DEBUG)}
  unset := 0;
  for i := 0 to High(Map) do
  begin
    if Map[i] = -1 then
      Inc(unset);
  end;
{$ifend}

  // linear interpolation where Map[i] == -1 (but no extrapolation)
  i := 0;
  LowerI := -1;
  LowerValue := -1;
{$if defined(MSWINDOWS) and defined(DEBUG)}
  interpolated := 0;
{$ifend}
  repeat
    if Map[i] = -1 then
    begin
      if LowerI <> -1 then
      begin
        UpperI := i+1;
        while (UpperI<=High(Map)) and (Map[UpperI] = -1) do
          Inc(UpperI);
        if UpperI<=High(Map) then
        begin
          UpperValue := Map[UpperI];
          for j := LowerI+1 to UpperI-1 do
          begin
            Map[j] := LowerValue + (UpperValue-LowerValue) * (j-LowerI) / (UpperI - LowerI);
{$if defined(MSWINDOWS) and defined(DEBUG)}
            Inc(interpolated);
{$ifend}
          end;
        end;
      end;
    end
    else
    begin
      LowerI := i;
      LowerValue := Map[i];
    end;
    Inc(i);
  until i > High(Map);
{$if defined(MSWINDOWS) and defined(DEBUG)}
  OutputDebugString(PChar(Format(
    'TRadialDistortionTransformation.PrepareReverseMap: mapToSameIndex=%d. IndexOutOfRange=%d. %d out of %d map elements were uninitialized, %d of these were interpolated',
    [ mapToSameIndex, IndexOutOfRange, unset, High(Map), interpolated ])));
{$ifend}

  for i := 0 to High(Map) do
  begin
    if Map[i] = -1 then
      Map[i] := 1;
  end;
{$if defined(MSWINDOWS) and defined(DEBUG)}
  OutputDebugString(PChar(Format('TRadialDistortionTransformation.PrepareReverseMap: MinValue(Map)=%f MaxValue(Map)=%f', [ MinValue(Map), MaxValue(Map) ])));
{$ifend}
end;

procedure TRadialDistortionTransformation.PrepareTransform;
var
  r: TRect;
begin
  if GR32.IsRectEmpty(SrcRect) then
    raise Exception.Create(RCStrSrcRectIsEmpty);
  TransformValid := not GR32.IsRectEmpty(SrcRect);
  if Not TransformValid then
    Exit;

  // center / focal point relative to which all (un)distortions are calculated
  FFocalPoint.x := (SrcRect.Right + SrcRect.Left) / 2;
  FFocalPoint.y := (SrcRect.Bottom + SrcRect.Top) / 2;

  r := MakeRect(SrcRect);
  r_0 := Sqrt(2*Sqr(Min(r.Right - r.Left, r.Bottom - r.Top)))/2;

  PrepareReverseMap;
end;

function TRadialDistortionTransformation.LookUpReverseMap(const r_tgt: TFloat): TFloat;
var
  index: Integer;
begin
  index := Trunc((r_tgt-r_tgt_min) / (r_tgt_max-r_tgt_min) * High(Map));
  if not InRange(index, 0, High(Map)) then
    raise Exception.Create(Format('TRadialDistortionTransformation.LookUpReverseMap: Index %d out of range (0..%d)', [ index, MapElements ]));
  Result := Map[index];
end;

procedure TRadialDistortionTransformation.ReverseTransformFloat(DstX, DstY: TFloat;
  out SrcX, SrcY: TFloat);
var
  r_tgt, r_src: Single;
  d: TFloatPoint;
begin
  d.x := DstX;
  d.y := DstY;
  r_tgt := Distance(FFocalPoint, d)/r_0;

  r_src := LookUpReverseMap(r_tgt);

  SrcX := FFocalPoint.X + (d.X-FFocalPoint.X) / r_src;
  SrcY := FFocalPoint.Y + (d.Y-FFocalPoint.Y) / r_src;
end;

procedure TRadialDistortionTransformation.SetCoefficient1(const Value: TFloat);
begin
  FCoefficient1 := Value;
  Changed;
end;

procedure TRadialDistortionTransformation.SetCoefficient2(const Value: TFloat);
begin
  FCoefficient2 := Value;
  Changed;
end;

procedure TRadialDistortionTransformation.SetScale(const Value: TFloat);
begin
  FScale := Value;
  Changed;
end;

procedure TRadialDistortionTransformation.SetMapElements(const Value: Integer);
begin
  FMapElements := Value;
  Changed;
end;

procedure TRadialDistortionTransformation.TransformFloat(SrcX, SrcY: TFloat; out DstX, DstY: TFloat);
var
  r_tgt, r_src: Single;
  d: TFloatPoint;
begin
  d.x := SrcX;
  d.y := SrcY;
  r_src := Distance(FFocalPoint, d)/r_0;
  r_tgt := Scale*(1 + FCoefficient1 * Sqr(r_src) + FCoefficient2 * Power(r_src, 4));
  DstX := FFocalPoint.X + (d.X-FFocalPoint.X) * r_tgt;
  DstY := FFocalPoint.Y + (d.Y-FFocalPoint.Y) * r_tgt;
end;


//------------------------------------------------------------------------------
//
//      TRemapTransformation
//
//------------------------------------------------------------------------------
constructor TRemapTransformation.Create;
begin
  inherited;
  FScalingFixed := FixedPoint(1, 1);
  FScalingFloat := FloatPoint(1, 1);
  FOffset := FloatPoint(0,0);
  FVectorMap := TVectorMap.Create;
  // Ensuring initial setup to avoid exceptions
  FVectorMap.SetSize(1, 1);
end;

destructor TRemapTransformation.Destroy;
begin
  FVectorMap.Free;
  inherited;
end;

//------------------------------------------------------------------------------

function TRemapTransformation.GetTransformedBounds(const ASrcRect: TFloatRect): TFloatRect;
const
  InfRect: TFloatRect = (Left: -Infinity; Top: -Infinity; Right: Infinity; Bottom: Infinity);
begin
  // We can't predict the ultimate bounds without transforming each vector in
  // the vector map, return the absolute biggest possible transformation bounds
  Result := InfRect;
end;

function TRemapTransformation.HasTransformedBounds: Boolean;
begin
  Result := False;
end;

procedure TRemapTransformation.PrepareTransform;
begin
  if GR32.IsRectEmpty(SrcRect) then
    raise Exception.Create(RCStrSrcRectIsEmpty);
  if GR32.IsRectEmpty(FMappingRect) then
    raise Exception.Create(RCStrMappingRectIsEmpty);
  with SrcRect do
  begin
    FSrcTranslationFloat.X := Left;
    FSrcTranslationFloat.Y := Top;
    FSrcScaleFloat.X := (Right - Left) / (FVectorMap.Width - 1);
    FSrcScaleFloat.Y := (Bottom - Top) / (FVectorMap.Height - 1);
    FSrcTranslationFixed := FixedPoint(FSrcTranslationFloat);
    FSrcScaleFixed := FixedPoint(FSrcScaleFloat);
  end;

  with FMappingRect do
  begin
    FDstTranslationFloat.X := Left;
    FDstTranslationFloat.Y := Top;
    FDstScaleFloat.X := (FVectorMap.Width - 1) / (Right - Left);
    FDstScaleFloat.Y := (FVectorMap.Height - 1) / (Bottom - Top);
    FCombinedScalingFloat.X := FDstScaleFloat.X * FScalingFloat.X;
    FCombinedScalingFloat.Y := FDstScaleFloat.Y * FScalingFloat.Y;
    FCombinedScalingFixed := FixedPoint(FCombinedScalingFloat);
    FDstTranslationFixed := FixedPoint(FDstTranslationFloat);
    FDstScaleFixed := FixedPoint(FDstScaleFloat);
  end;
  TransformValid := True;
end;

procedure TRemapTransformation.ReverseTransformFixed(DstX, DstY: TFixed;
  out SrcX, SrcY: TFixed);
begin
  with FVectorMap.FixedVectorX[DstX - FOffsetFixed.X, DstY - FOffsetFixed.Y] do
  begin
    DstX := DstX - FDstTranslationFixed.X;
    DstX := FixedMul(DstX , FDstScaleFixed.X);
    DstX := DstX + FixedMul(X, FCombinedScalingFixed.X);
    DstX := FixedMul(DstX, FSrcScaleFixed.X);
    SrcX := DstX + FSrcTranslationFixed.X;

    DstY := DstY - FDstTranslationFixed.Y;
    DstY := FixedMul(DstY, FDstScaleFixed.Y);
    DstY := DstY + FixedMul(Y, FCombinedScalingFixed.Y);
    DstY := FixedMul(DstY, FSrcScaleFixed.Y);
    SrcY := DstY + FSrcTranslationFixed.Y;
  end;
end;

procedure TRemapTransformation.ReverseTransformFloat(DstX, DstY: TFloat;
  out SrcX, SrcY: TFloat);
begin
  with FVectorMap.FloatVectorF[DstX - FOffset.X, DstY - FOffset.Y] do
  begin
    DstX := DstX - FDstTranslationFloat.X;
    DstY := DstY - FDstTranslationFloat.Y;
    DstX := DstX * FDstScaleFloat.X;
    DstY := DstY * FDstScaleFloat.Y;

    DstX := DstX + X * FCombinedScalingFloat.X;
    DstY := DstY + Y * FCombinedScalingFloat.Y;

    DstX := DstX * FSrcScaleFloat.X;
    DstY := DstY * FSrcScaleFloat.Y;
    SrcX := DstX + FSrcTranslationFloat.X;
    SrcY := DstY + FSrcTranslationFloat.Y;
  end;
end;

procedure TRemapTransformation.ReverseTransformInt(DstX, DstY: Integer;
  out SrcX, SrcY: Integer);
begin
  with FVectorMap.FixedVector[DstX - FOffsetInt.X, DstY - FOffsetInt.Y] do
  begin
    DstX := DstX * FixedOne - FDstTranslationFixed.X;
    DstY := DstY * FixedOne - FDstTranslationFixed.Y;
    DstX := FixedMul(DstX, FDstScaleFixed.X);
    DstY := FixedMul(DstY, FDstScaleFixed.Y);

    DstX := DstX + FixedMul(X, FCombinedScalingFixed.X);
    DstY := DstY + FixedMul(Y, FCombinedScalingFixed.Y);

    DstX := FixedMul(DstX, FSrcScaleFixed.X);
    DstY := FixedMul(DstY, FSrcScaleFixed.Y);
    SrcX := FixedRound(DstX + FSrcTranslationFixed.X);
    SrcY := FixedRound(DstY + FSrcTranslationFixed.Y);
  end;
end;

procedure TRemapTransformation.Scale(Sx, Sy: TFloat);
begin
  FScalingFixed.X := Fixed(Sx);
  FScalingFixed.Y := Fixed(Sy);
  FScalingFloat.X := Sx;
  FScalingFloat.Y := Sy;
  Changed;  
end;

procedure TRemapTransformation.SetMappingRect(Rect: TFloatRect);
begin
  FMappingRect := Rect;
  Changed;
end;

procedure TRemapTransformation.SetOffset(const Value: TFloatVector);
begin
  FOffset := Value;
  FOffsetInt := Point(Value);
  FOffsetFixed := FixedPoint(Value);
  Changed;
end;


//------------------------------------------------------------------------------
//
//      TSphereTransformation
//
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
// Utilities
//------------------------------------------------------------------------------
procedure Modulo2Pi(var Angle: TFloat); {Result is between 0 and 2PI }
{$if defined(PUREPASCAL) or (not defined(TARGET_x86))}
begin
  Angle := GR32_Math.FloatMod(Angle, PI*2);
end;
{$else}
asm
        FLDPI
        FADD    ST,ST               // 2PI
        FLD     DWORD ptr [Angle]
        FPREM                       // calc Modulo
        FLDZ
        FCOMIP  ST,ST(1)            // Compare 0 and Modulo (+pop 0)
        JNB     @@1                 // if Modulo >= 0 then
        FSTP    DWORD ptr [Angle]   //  return Modulo...
        FSTP    ST(0)               //  POP the rest (2PI)
        JMP     @@2
@@1:    FADDP                       // add Modulo and Rest (2PI)
        FSTP    DWORD ptr [Angle]   // Modulo+2PI...
@@2:    FWAIT
end;
{$ifend}


constructor TSphereTransformation.Create;
begin
  inherited;
  FRadius := 1;
end;

//------------------------------------------------------------------------------

function TSphereTransformation.GetTransformedBounds(const ASrcRect: TFloatRect): TFloatRect;
begin
  // There is not direct relation between SourceRect and DestRect !
  // During transformation process this rect will be clipped.
  Result.Left := FCenter.X - FRadius;
  Result.Top := FCenter.Y - FRadius;
  Result.Bottom := FCenter.Y + FRadius;
  Result.Right := FCenter.X + FRadius;
end;

function TSphereTransformation.HasTransformedBounds: Boolean;
begin
  Result := False;
end;

function TSphereTransformation.IsInSphere(CartesianX, CartesianY: TFloat): boolean;
begin
  if not TransformValid then
    PrepareTransform;

  CartesianX := CartesianX - FCenter.X;
  CartesianY := CartesianY - FCenter.Y;

  Result := (FSquareRadius >= (CartesianX * CartesianX + CartesianY * CartesianY));
end;

procedure TSphereTransformation.PrepareTransform;
begin
  // Invariants during transformation
  FMapWidth := (SrcRect.Width - 1) / (2 * PI);
  FMapHeight := (SrcRect.Height - 1) / PI;
  FSquareRadius := Sqr(FRadius);
  GR32_Math.SinCos(FLattitude, FLattitudeSin, FLattitudeCos);
  FLattitudeSinInvRadius := -FLattitudeSin / FRadius;
  FLattitudeCosInvRadius := FLattitudeCos / FRadius;
  FSrcRectTop := SrcRect.Top;
  FSrcRectLeft := SrcRect.Left;

  TransformValid := True;
end;

procedure TSphereTransformation.ReverseTransformFloat(DstX, DstY: TFloat; out SrcX, SrcY: TFloat);
// FPC currently refuses to compile the ASM version.
// Consider deprecating it as it's not really worth the effort - or replace it with a SSE version
{$if defined(PUREPASCAL) or (not defined(TARGET_x86)) or (defined(FPC))}
var
  Dist: TFloat;
begin
  // Screen projection on sphere
  DstX := DstX - FCenter.X; // = Y
  DstY := FCenter.Y - DstY; // = Z
  Dist := DstX * DstX + DstY * DstY;

  if (Dist > FSquareRadius) then // Not projectable on the sphere
  begin
    SrcX := -1;
    SrcY := -1;
    Exit;
  end;

  Dist := Sqrt(FSquareRadius - Dist);

  // Apply rotations
  DstX := Arctan2(DstX, Dist * FLattitudeCos + DstY * FLattitudeSin) + FLongitude;
  Modulo2Pi(DstX);

  // Map projection
  SrcX := SrcRect.Left + DstX * FMapWidth;
  SrcY := SrcRect.Top + ArcCos(DstY * FLattitudeCosInvRadius + Dist * FLattitudeSinInvRadius) * FMapHeight;
end;
{$else}
{Assembler version (FPU) ... 4% faster on a P4 }
asm
// screen projection on sphere
//  DstX := DstX - FCenterX; // = Y
    fld   DstX               // DstX
    fsub  [eax].FCenter.X // DstX'
//  DstY := FCenterY - DstY; // = Z
    fld   [eax].FCenter.Y // FCenterY | DstX'
    fsub  DstY               // DstY'    | DstX'
//  x := DstX * DstX + DstY * DstY;
    fld   st(0)                    // Z    | Z    | Y
    fmul  st(0),st(1)              // ZZ   | Z    | Y
    fld   st(2)                    // Y    | ZZ   | Z   | Y
    fmul  st(0),st(3)
    faddp                          // X'   | Z    | Y
//  if (FSquareRadius < x) then // not projetable in the sphere.
    fld [eax].FSquareRadius
    fcomi st(0),st(1) // st(0) < st(1)
    jnbe   @@1
    fstp  st(0)
    fstp  st(0)
    fstp  st(0)
    fstp  st(0)
//    SrcX := -1;
    mov   [SrcX],$bf800000
//    SrcY := -1;
    mov   [SrcY],$bf800000
//    Exit;
    jmp @@fin
@@1:
//  x := sqrt(FSquareRadius - x);
    fsubrp
    fsqrt                          // X    | Z    | Y
// apply rotations
//  DstX := Arctan2(Y,X * FLattitudeCos + Z * FLattitudeSin) + FLongitude; // Lon
    fxch  st(2)                   // Y     | Z    | Y
    fld   st(2)                   // X     | Y    | Z    | X
    fmul  [eax].FLattitudeCos
    fld   st(2)                   // Z     | Xx.  | Y    | Z    | X
    fmul  [eax].FLattitudeSin     // Zx.   | Xx.  | Y    | Z    | X
    faddp                         // Xx+Zx | Y    | Z    | X
    fpatan
    fadd  [eax].FLongitude  // DstX  | Z    | X
//  if DstX > PI2S then
    fldpi
    fadd  st(0),st(0)       // 2PI   | DstX | Z    | X
    fcomi st(0),st(1) // st(0) < st(1)
    jnb   @@test2
//    DstX := DstX - PI2S
    fsubp st(1),st(0)
    jmp   @@testfin
//  else if DstX < 0 then
@@test2:
    fldz
    fcomip st(0),st(2) // st(0) < st(2)
    jb @@test3
//    DstX := DstX + PI2S;
    faddp
    jmp   @@testfin
@@test3:
    fstp  st(0)
@@testfin:
// Map projection
//  SrcX := DstX * FMapWidth;
    fmul  [eax].FMapWidth
    FADD  [eax].FSrcRectLeft
    fstp  dword ptr [SrcX]// Z    | X
//  SrcY := ArcCos(Z * FLattitudeCosInvRadius + x * FLattitudeSinInvRadius) * FMapHeight;
    fmul  [eax].FLattitudeCosInvRadius
    fxch
    fmul  [eax].FLattitudeSinInvRadius
    faddp
    FLD1                 // 1      | X
    FLD   ST(1)          // X      | 1      | X
    FMUL  ST(0),ST(0)    // X²     | 1      | X
    FSUBP ST(1),ST(0)    // 1 - X² | X
    FABS                 //<- avoid rounding errors...
    FSQRT                // sqrt(.)| X
    FXCH  st(1)
    FPATAN               // result |
    fmul  [eax].FMapHeight
    FADD  [eax].FSrcRectTop
    fstp  dword ptr [SrcY]
@@fin:
    fwait
end;
{$ifend}

function TSphereTransformation.ScreenCoordinate(var X, Y: TFloat): boolean;
var
  SinLong, CosLong, SinLat, CosLat: TFloat;
begin
  if not TransformValid then
    PrepareTransform;

  GR32_Math.SinCos(X - FLongitude, SinLong, CosLong);
  GR32_Math.SinCos(Y, SinLat, CosLat);

  Result := (SinLat * CosLong * FLattitudeCos >= CosLat * FLattitudeSin);

  if Result then
  begin
    X := FCenter.X + FRadius * SinLat * SinLong;
    Y := FCenter.Y - FRadius * (SinLat * CosLong * FLattitudeSin + CosLat * FLattitudeCos);
  end;
end;

procedure TSphereTransformation.SetCenter(const Value: TFloatPoint);
begin
  if FCenter <> Value then
  begin
    FCenter := Value;
    Changed;
  end;
end;

procedure TSphereTransformation.SetLattitude(const Value: TFloat);
begin
  if FLattitude <> Value then
  begin
    FLattitude := Value;
    Modulo2Pi(FLattitude);
    Changed;
  end;
end;

procedure TSphereTransformation.SetLongitude(const Value: TFloat);
begin
  if FLongitude <> Value then
  begin
    FLongitude := Value;
    Modulo2Pi(FLongitude);
    Changed;
  end;
end;

procedure TSphereTransformation.SetRadius(const Value: TFloat);
begin
  if (Value > 0) and (FRadius <> Value) then
  begin
    FRadius := Value;
    Changed;
  end;
end;

function TSphereTransformation.SphericalCoordinate(CartesianX, CartesianY: TFloat): TFloatPoint;
var
  Dist: TFloat;
begin
  if not TransformValid then
    PrepareTransform;

  // Screen projection on sphere
  CartesianX := CartesianX - FCenter.X; // = Y
  CartesianY := FCenter.Y - CartesianY; // = Z
  Dist := CartesianX * CartesianX + CartesianY * CartesianY;

  if (Dist > FSquareRadius) then // Not projectable in the sphere.
  begin
    Result.X := 0;
    Result.Y := 0;
    Exit;
  end;

  Dist := Sqrt(FSquareRadius - Dist);

  // Apply rotations
  Result.X := Arctan2(CartesianX, Dist * FLattitudeCos + CartesianY * FLattitudeSin) + FLongitude;
  Modulo2Pi(Result.X);
  Result.Y := ArcCos(CartesianY * FLattitudeCosInvRadius + Dist * FLattitudeSinInvRadius) - (PI / 2);
end;


//------------------------------------------------------------------------------
//
//      Bindings
//
//------------------------------------------------------------------------------

var
  TransformsRegistry: TFunctionRegistry;

procedure RegisterBindings;
begin
  TransformsRegistry := NewRegistry('GR32_Transforms bindings');
  TransformsRegistry.RegisterBinding(@@DET_2x2_32, 'DET_2x2_32');
  TransformsRegistry.RegisterBinding(@@DET_3x3_32, 'DET_3x3_32');
  TransformsRegistry.RegisterBinding(@@DET_2x2_64, 'DET_2x2_64');

  // DET_2x2_32
  TransformsRegistry[@@DET_2x2_32].Add(@DET_2x2_32_Pas, [isPascal]).Name := 'DET_2x2_32_Pas';
{$IFNDEF PUREPASCAL}
{$if defined(TARGET_x86)}
  TransformsRegistry[@@DET_2x2_32].Add(@DET_2x2_32_ASM, [isAssembler]).Name := 'DET_2x2_32_ASM';
{$elseif defined(TARGET_x64) and (not defined(OMIT_SSE2))}
  TransformsRegistry[@@DET_2x2_32].Add(@DET_2x2_32_SSE, [isSSE]).Name := 'DET_2x2_32_SSE';
{$ifend}
{$ENDIF}

  // DET_2x2_64
  TransformsRegistry[@@DET_2x2_64].Add(@DET_2x2_64_Pas, [isPascal]).Name := 'DET_2x2_64_Pas';
{$IFNDEF PUREPASCAL}
{$if defined(TARGET_x86)}
  TransformsRegistry[@@DET_2x2_64].Add(@DET_2x2_64_ASM, [isAssembler]).Name := 'DET_2x2_64_ASM';
{$elseif defined(TARGET_x64) and (not defined(OMIT_SSE2))}
  TransformsRegistry[@@DET_2x2_64].Add(@DET_2x2_64_SSE, [isSSE]).Name := 'DET_2x2_64_SSE';
{$ifend}
{$ENDIF}

  // DET_3x3_32
  TransformsRegistry[@@DET_3x3_32].Add(@DET_3x3_32_Pas, [isPascal]).Name := 'DET_3x3_32_Pas';

  TransformsRegistry.RebindAll;
end;

//------------------------------------------------------------------------------

initialization
  RegisterBindings;

end.