{
     Apophysis Copyright (C) 2001-2004 Mark Townsend
     Apophysis Copyright (C) 2005-2006 Ronald Hordijk, Piotr Borys, Peter Sdobnov
     Apophysis Copyright (C) 2007-2008 Piotr Borys, Peter Sdobnov
     
     Apophysis "3D hack" Copyright (C) 2007-2008 Peter Sdobnov
     Apophysis "7X" Copyright (C) 2009-2010 Georg Kiehne
     Apophysis AV "Phoenix Edition" Copyright (C) 2021-2022 Alice V. Koryagina

     This program is free software; you can redistribute it and/or modify
     it under the terms of the GNU General Public License as published by
     the Free Software Foundation; either version 2 of the License, or
     (at your option) any later version.

     This program is distributed in the hope that it will be useful,
     but WITHOUT ANY WARRANTY; without even the implied warranty of
     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
     GNU General Public License for more details.

     You should have received a copy of the GNU General Public License
     along with this program; if not, write to the Free Software
     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
}

unit ControlPoint;

interface

//{$define VAR_STR}
//{$define GAUSSIAN_DOF}

uses
  Classes, Windows, Cmap, XForm, XFormMan, Binary, SysUtils,
  Math, Bezier (* {$ifdef CPUX86}, AsmRandom {$endif} *);

const
  SUB_BATCH_SIZE = 10000;
  PROP_TABLE_SIZE = 1024;
  PREFILTER_WHITE = (1 shl 26);
  FILTER_CUTOFF = 1.8;  // AV: maybe move it to ImageMaker?
  BRIGHT_ADJUST = 2.3;
  //FUSE = 15; // AV: moved to Global since it became a variable

  vRandom = -1;  // AV: index of randomly chosen variation for random flames

type

  TTriangle = record
    x: array[0..2] of double;
    y: array[0..2] of double;
  end;
  TTriangles = array[-1..NXFORMS] of TTriangle;

  TSPoint = record
    x: double;
    y: double;
  end;
  TSRect = record
    Left, Top, Right, Bottom: double;
  end;

{
  TMapPalette = record
    Red:   array[0..255] of byte;
    Green: array[0..255] of byte;
    Blue:  array[0..255] of byte;
  end;

  TColorMaps = record
    Identifier: string;
    UGRFile: string;
  end;
}

type

  TPointsArray = array of TCPpoint;
  TPointsXYArray = array of TXYpoint;
  //T2CPointsArray = array of T2Cpoint;

  TControlPoint = class
  public
    finalXform: TXForm;
    finalXformEnabled: boolean;
    useFinalXform: boolean;
    soloXform: integer;
    curvePoints: array [0..3] of BezierPoints;
    curveWeights: array [0..3] of BezierWeights;

    Transparency: boolean;

    cameraPitch, cameraYaw, cameraPersp, cameraDOF: double;
    cameraZpos, cameraRoll: double;
    ProjectionFunc: procedure(pPoint: PCPPoint) of object;

    xform: array[0..NXFORMS] of TXForm;

    noLinearFix: boolean;
    cmap: TColorMap;
    cmapindex: integer;
    time: double;  // AV: is used for interpolation
    Fbrightness: double; // 1.0 = normal
    contrast: double; // 1.0 = normal
    gamma: double;
    Width: integer;
    Height: integer;
    spatial_oversample: integer;
    name: string;
    center: array[0..1] of double; // camera center
    vibrancy: double; // blend between color algs (0=old,1=new)
    hue_rotation: double; // applies to cmap, 0-1
    background: array[0..3] of smallint; // Changed to integers so no conversion needed - mt
    zoom: double; // effects ppu and sample density
    pixels_per_unit: double; // and scale
    spatial_filter_radius: double; // variance of gaussian
    sample_density: extended; // samples per pixel (not bucket)
    (* in order to motion blur more accurately we compute the logs of the
      sample density many times and average the results.  we interplate
      only this many times. *)
    actual_density: extended; // for incomplete renders
    nbatches: integer; // this much color resolution.  but making it too high induces clipping
    white_level: integer;

    FAngle: Double;
    symmetry: integer; // color speed

    estimator, estimator_min, estimator_curve: double; // density estimator.
    jitters: integer;  // <-- AV: deprecated
    gamma_threshold: double;
    enable_de : boolean;
    used_plugins : TStringList;
    comment: string; // AV: holds user's comment on flame

    {
    pulse: array[0..1, 0..1] of double; // [i][0]=magnitute [i][1]=frequency */
    wiggle: array[0..1, 0..1] of double; // frequency is /minute, assuming 30 frames/s */
    cmap_inter: integer; // if this is true, then color map interpolates one entry
                      // at a time with a bright edge
    variation:  TVariation;  // <-- AV: deprecated
    nick, url: string;  // <-- AV: deprecated
    PropTable: array of TXForm; // <-- AV: declared in TXform unit
    FTwoColorDimensions: Boolean;
    xdata : string;
    }

    procedure FillUsedPlugins;

  private
    invalidXform: TXForm;

    CameraMatrix: array[0..2, 0..2] of double;
    DofCoef: double;
    {$ifdef GAUSSIAN_DOF}
    gauss_rnd: array [0..3] of double;
    gauss_N: integer;
    {$endif}

    sinPitch, cosPitch, sinRoll, cosRoll,
    sinYaw, cosYaw: double;

    procedure ProjectNone(pPoint: PCPPoint);
    procedure ProjectPitch(pPoint: PCPPoint);
    procedure ProjectPitchYaw(pPoint: PCPPoint);
    procedure ProjectPitchDOF(pPoint: PCPPoint);
    procedure ProjectPitchYawDOF(pPoint: PCPPoint);
    // AV: for rotation around Y-axis
    procedure ProjectPitchRoll(pPoint: PCPPoint);
    procedure ProjectPitchRollYaw(pPoint: PCPPoint);
    procedure ProjectPitchRollYawDOF(pPoint: PCPPoint);

    function getppux: double;
    function getppuy: double;

    procedure PreCalcBounds(var maxx, minx, maxy, miny: double); // AV
    procedure DirectCopy_AV(const cp1: TControlPoint);

    function GetBrightness: double;
    procedure SetBrightness(br: double);
    function GetRelativeGammaThreshold: double;
    procedure SetRelativeGammaThreshold(gtr: double);

  public

    procedure SaveToStringlist(sl: TStringlist);
    procedure SaveToFile(Filename: string);
    procedure SaveToBinary(const handle: File);

    procedure ParseString(aString: string);
    procedure ParseStringList(sl: TStringlist);
    procedure CalcBoundbox;
    function BlowsUp(NrPoints: integer): boolean;
	
    procedure Clear;

    procedure InterpolateX(cp1, cp2: TControlPoint; Tm: double);
    procedure InterpolateAll(cp1, cp2: TControlPoint; Tm: double; it: integer);

    procedure IterateXY(NrPoints: integer; var Points: TPointsXYArray);
    procedure IterateXYC(NrPoints: integer; var Points: TPointsArray);

    procedure Prepare;

    function Clone: TControlPoint;  // AV: rewritten
    procedure Copy(cp1: TControlPoint; KeepSizes: boolean = false); // AV: rewritten

    function HasFinalXForm: boolean;

    // CP-specific functions moved from unit Main
    function NumXForms: integer;
    function TrianglesFromCP(var Triangles: TTriangles): integer;
    procedure GetFromTriangles(const Triangles: TTriangles; const t: integer);

    procedure GetTriangle(var Triangle: TTriangle; const n: integer);
    procedure GetPostTriangle(var Triangle: TTriangle; const n: integer);

    procedure NormalizeProbabilities;  // AV
    procedure CalculateWeights;  // AV
    procedure CalculateColorSpeed;  // AV
    procedure AdjustScale(w, h: integer);

    constructor Create;
    destructor Destroy; override;

    procedure ZoomtoRect(R: TSRect);
    procedure ZoomOuttoRect(R: TSRect);
    procedure MoveRect(R: TSRect);
    procedure ZoomIn(Factor: double);
    procedure Rotate(Angle: double);

    property ppux: double read getppux;
    property ppuy: double read getppuy;

    property brightness: double
      read GetBrightness
      write SetBrightness;
    property gammaThreshRelative: double
      read GetRelativeGammaThreshold
      write SetRelativeGammaThreshold;
  end;

function add_symmetry_to_control_point(var cp: TControlPoint; sym: integer): integer;
procedure tile_control_point(var cp: TControlPoint; sym: integer);

function CalcUPRMagn(const cp: TControlPoint): double;
function CalcBinaryFlameSize(cp: TControlPoint): integer;
procedure PrepareToInterpolation(var SourceCp, TargetCp: TControlPoint); // AV

{ Math operations }
function line_dist(x, y, x1, y1, x2, y2: double): double;
function dist(x1, y1, x2, y2: double): double;
function det(a, b, c, d: double): double;
function solve3(x1, x2, x1h, y1, y2, y1h, z1, z2, z1h: double;
                var a, b, e: double): double;

// AV: moved specific variables from Global here due to often name-space conflicts
var
  MainTriangles:  TTriangles;
  Variation: integer; // Current variation

implementation

uses Global;

{ TControlPoint }

procedure TControlPoint.FillUsedPlugins;
var
  i, j, f : integer;
  v : double;
  s : string;
begin
  used_plugins.Clear;

  if self.finalXformEnabled then  // AV
    f := Min(NumXForms, NXFORMS)
  else
    f := Min(NumXForms - 1, NXFORMS);

  for i := 0 to f do
    with xform[i] do begin
      for j := 0 to NRVAR - 1 do begin
        v := self.xform[i].GetVariation(j);
        s := Varnames(j);
        if (v <> 0) and // uses variation
          (used_plugins.IndexOf(s) < 0)  // not listed yet
         then used_plugins.Add(s);
      end;
    end;

end;

constructor TControlPoint.Create;
var
  i: word;
begin
  for i := 0 to NXFORMS do
    xform[i] := TXForm.Create;

  invalidXform := TXForm.Create;
  soloXform := -1;

  background[0] := 0;
  background[1] := 0;
  background[2] := 0;

  for i := 0 to 3 do
  begin
    curvePoints[i][0].x := 0.00; curvePoints[i][0].y := 0.00; curveWeights[i][0] := 1;
    curvePoints[i][1].x := 0.00; curvePoints[i][1].y := 0.00; curveWeights[i][1] := 1;
    curvePoints[i][2].x := 1.00; curvePoints[i][2].y := 1.00; curveWeights[i][2] := 1;
    curvePoints[i][3].x := 1.00; curvePoints[i][3].y := 1.00; curveWeights[i][3] := 1;
  end;

  center[0] := 0;
  center[1] := 0;

  pixels_per_unit := 50;

  width := 100;
  Height := 100;

  spatial_oversample := 1;
  spatial_filter_radius := 0.5;

  FAngle := 0;
  gamma := 1;
  vibrancy := 1;
  contrast := 1;
  Fbrightness := 1;
  gamma_threshold := 0.01;  

  sample_density := 50;
  zoom := 0;
  
  estimator := 9.0;
  estimator_min := 0.0;
  estimator_curve := 0.4;
  enable_de := false;
  jitters := 1;
  nbatches := 1;
  
  white_level := 200;
  hue_rotation := 1; // AV
  
  finalXformEnabled := false;
  Transparency := false;

  cameraPitch := 0;
  cameraYaw := 0;
  cameraRoll := 0;
  cameraPersp := 0;
  cameraZpos := 0;
  cameraDOF := 0;

  used_plugins := TStringList.Create;
  comment := '';
  
  {
  pulse[0][0] := 0;
  pulse[0][1] := 60;
  pulse[1][0] := 0;
  pulse[1][1] := 60;

  wiggle[0][0] := 0;
  wiggle[0][1] := 60;
  wiggle[1][0] := 0;
  wiggle[1][1] := 60;

  xdata := '';
  FTwoColorDimensions := False;
  }
end;

destructor TControlPoint.Destroy;
var
  i: word;
begin
  for i := 0 to NXFORMS do
    xform[i].Free;
  invalidXform.Free;

  used_plugins.Free;  // <<< -X- commenting out = hack - fixme!
  // AV: uncommented to save user's PC memory from leaking - fixme!
  inherited;
end;

procedure TControlPoint.Prepare;
var
  i, n: Integer;
  propsum: double;
  LoopValue: double;
  j: integer;
  TotValue: double;

  k: integer;
  tp: array[0..NXFORMS] of double;
begin
//  SetLength(PropTable, PROP_TABLE_SIZE);

  //totValue := 0;
  n := NumXforms;
  assert(n > 0);

  finalXform := xform[n];
  finalXform.Prepare;
  useFinalXform := FinalXformEnabled and HasFinalXform;
  for i := 0 to n - 1 do begin
    xform[i].Prepare;
    //totValue := totValue + xform[i].density;
  end;
  invalidXform.PrepareInvalidXForm;

  if soloXform >= 0 then begin
    for i := 0 to n - 1 do begin
      xform[i].transOpacity := 0;
    end;
    xform[soloXform].transOpacity := 1;
  end;

  for k := 0 to n - 1 do begin
    totValue := 0;
    SetLength(xform[k].PropTable, PROP_TABLE_SIZE);

    for i := 0 to n - 1 do begin
      tp[i] := xform[i].density * xform[k].modWeights[i];
      totValue := totValue + tp[i];
    end;

    if totValue > 0 then begin
     LoopValue := 0;
      for i := 0 to PROP_TABLE_SIZE-1 do begin
        propsum := 0;
        j := -1;
        repeat
          inc(j);
          propsum := propsum + tp[j];//xform[j].density;
        until (propsum > LoopValue) or (j = n - 1);

      //assert(tp[j]<>0);

        xform[k].PropTable[i] := xform[j];
        LoopValue := LoopValue + TotValue / PROP_TABLE_SIZE;
      end;
    end
    else begin
      for i := 0 to PROP_TABLE_SIZE-1 do
        xform[k].PropTable[i] := invalidXform;
    end;
  end;

  (*
  CameraMatrix[0, 0] := cos(-CameraYaw);
  CameraMatrix[1, 0] := -sin(-CameraYaw);
  CameraMatrix[2, 0] := 0;
  CameraMatrix[0, 1] := cos(CameraPitch) * sin(-CameraYaw);
  CameraMatrix[1, 1] := cos(CameraPitch) * cos(-CameraYaw);
  CameraMatrix[2, 1] := -sin(CameraPitch);
  CameraMatrix[0, 2] := sin(CameraPitch) * sin(-CameraYaw);
  CameraMatrix[1, 2] := sin(CameraPitch) * cos(-CameraYaw);
  CameraMatrix[2, 2] := cos(CameraPitch);
   *)
   // AV: added 3D camera precalc
  SinCos(cameraPitch, sinPitch, cosPitch);
  SinCos(cameraRoll, sinRoll, cosRoll);
  SinCos(-cameraYaw, sinYaw, cosYaw);
   // AV: recalculate the matrix in order to support Y-rotation
  CameraMatrix[0, 0] := cosRoll * cosYaw;
  CameraMatrix[1, 0] := -cosRoll * sinYaw;
  CameraMatrix[2, 0] := sinRoll;
  CameraMatrix[0, 1] := sinPitch * sinRoll * cosYaw + cosPitch * sinYaw;
  CameraMatrix[1, 1] := -sinPitch * sinRoll * sinYaw + cosPitch * cosYaw;
  CameraMatrix[2, 1] := -sinPitch * cosRoll;
  CameraMatrix[0, 2] := -cosPitch * sinRoll * cosYaw + sinPitch * sinYaw;
  CameraMatrix[1, 2] := cosPitch * sinRoll * sinYaw + sinPitch * cosYaw;
  CameraMatrix[2, 2] := cosPitch * cosRoll;
  DofCoef := 0.1 * CameraDOF;
  {$ifdef GAUSSIAN_DOF}
  gauss_rnd[0] := random;
  gauss_rnd[1] := random;
  gauss_rnd[2] := random;
  gauss_rnd[3] := random;
  gauss_N := 0;
  {$endif}

  if (CameraDOF <> 0) then begin
    if (CameraRoll = 0) then
    begin
      if (CameraYaw <> 0) then
        ProjectionFunc := ProjectPitchYawDOF
      else
        ProjectionFunc := ProjectPitchDOF;
    end else  // Roll and DOF  
        ProjectionFunc := ProjectPitchRollYawDOF;
  end
  else if (CameraPitch <> 0) or (CameraRoll <> 0) or (CameraYaw <> 0) then begin
    if (CameraRoll <> 0) then
      ProjectionFunc := ProjectPitchRollYaw
    else if (CameraYaw <> 0) then  // Roll  = 0
      ProjectionFunc := ProjectPitchYaw
    else  // Roll = Yaw = 0
      ProjectionFunc := ProjectPitch;
  end
  else ProjectionFunc := ProjectNone;
end;

procedure TControlPoint.IterateXY(NrPoints: integer; var Points: TPointsXYArray);
var
  i: Integer;
  px, py: double;
  pPoint: PXYPoint;

  xf: TXform;
begin
  px := 2 * random - 1;
  py := 2 * random - 1;

  try
    xf := xform[0]; //random(NumXForms)];
    for i := 0 to FUSE do begin
      xf := xf.PropTable[Random(PROP_TABLE_SIZE)];
      xf.NextPointXY(px,py);
    end;

    pPoint := @Points[0];

    if UseFinalXform then
      for i := 0 to NrPoints - 1 do begin
        xf := xf.PropTable[Random(PROP_TABLE_SIZE)];
        xf.NextPointXY(px,py);
		    if (xf.transOpacity = 0) or
          ((xf.transOpacity < 1) and (random > xf.transOpacity)) then
			    pPoint^.x := 1e300 // MaxDouble // hack
		    else begin
            pPoint^.X := px;
            pPoint^.Y := py;
        end;
        finalXform.NextPointXY(pPoint^.X, pPoint^.y);
        Inc(pPoint);
      end
      else
        for i := 0 to NrPoints - 1 do begin
          xf := xf.PropTable[Random(PROP_TABLE_SIZE)];
          xf.NextPointXY(px,py);
          if (xf.transOpacity = 0) or   // AV
            ((xf.transOpacity < 1) and (random > xf.transOpacity)) then  // AV
             pPoint^.x := 1e300 // MaxDouble // hack
          else begin
            pPoint.X := px;
            pPoint.Y := py;
          end;
          Inc(pPoint);
        end;
  except
    on EMathError do begin
      exit;
    end;
  end;
end;

procedure TControlPoint.IterateXYC(NrPoints: integer; var Points: TPointsArray);
var
  i: Integer;
  p: TCPPoint;
  pPoint: PCPPoint;

  xf: TXform;
begin
  p.x := 2 * random - 1;
  p.y := 2 * random - 1;
  p.z := 0; // AV: fixed - someone forgot to initialize this variable
  p.c := random;

  try
    xf := xform[0]; //random(NumXForms)];
    for i := 0 to FUSE do begin
      xf := xf.PropTable[Random(PROP_TABLE_SIZE)];
      xf.NextPoint(p);
    end;

    pPoint := @Points[0];

    if UseFinalXform then
      for i := 0 to NrPoints - 1 do begin
        xf := xf.PropTable[Random(PROP_TABLE_SIZE)];
        xf.NextPoint(p);

        if (xf.transOpacity = 0) or  // AV: moved here from rendering unit
          ((xf.transOpacity < 1) and (random > xf.transOpacity)) then
          pPoint^.x := 1e300 // MaxDouble // hack
        else
		    begin
			    finalXform.NextPointTo(p, pPoint^);
		      ProjectionFunc(pPoint);
		    end;
        Inc(pPoint);
      end
    else
      for i := 0 to NrPoints - 1 do begin
        xf := xf.PropTable[Random(PROP_TABLE_SIZE)];
        xf.NextPoint(p);

        if (xf.transOpacity = 0) or   // AV: moved here from rendering unit
          ((xf.transOpacity < 1) and (random > xf.transOpacity)) then
          pPoint^.x := 1e300 // MaxDouble // hack
        else 
		    begin
          pPoint^ := p;
          ProjectionFunc(pPoint);
        end;
        Inc(pPoint);
      end;
  except
    on EMathError do begin
      exit;
    end;
    on EInvalidPointer do exit; // AV  
  end;
end;

procedure TControlPoint.ProjectNone(pPoint: PCPPoint);
var
  zr: double;
begin
  zr := 1 - cameraPersp * (pPoint^.z - CameraZpos);

  pPoint^.x := pPoint^.x / zr;
  pPoint^.y := pPoint^.y / zr;
  //pPoint^.z := pPoint^.z - CameraZpos;	// AV: ?
end;

procedure TControlPoint.ProjectPitch(pPoint: PCPPoint);
var
  y, z, zr: double;
begin
  z := pPoint^.z - CameraZpos;
  y := CameraMatrix[1,1] * pPoint^.y  - sinPitch * z;
  zr := 1 - cameraPersp *
    (CameraMatrix[1,2] * pPoint^.y + cosPitch * z);

  pPoint^.x := pPoint^.x / zr;
  pPoint^.y := y / zr;
  //pPoint^.z := pPoint^.z - CameraZpos; 
end;

procedure TControlPoint.ProjectPitchYaw(pPoint: PCPPoint);
var
  x, y, z, zr: double;
begin
  z := pPoint^.z - CameraZpos;
  x := cosYaw * pPoint^.x - sinYaw * pPoint^.y;
  y := CameraMatrix[0,1] * pPoint^.x + CameraMatrix[1,1] * pPoint^.y - sinPitch * z;
  zr := 1 - cameraPersp *
    (CameraMatrix[0,2]* pPoint^.x + CameraMatrix[1,2]*pPoint^.y + cosPitch * z);

  pPoint^.x := x / zr;
  pPoint^.y := y / zr;
  //pPoint^.z := pPoint^.z - CameraZpos;  
end;

procedure TControlPoint.ProjectPitchDOF(pPoint: PCPPoint);
var
  x, y, z, zr, dr: double;
  dsin, dcos: double;
  t: double;
begin
  z := pPoint^.z - CameraZpos;
  y := CameraMatrix[1,1] * pPoint^.y   - sinPitch * z;
  z := CameraMatrix[1,2] * pPoint^.y + cosPitch * z;
  zr := 1 - cameraPersp * z;

{$ifdef GAUSSIAN_DOF}
  asm
    fld     qword ptr [eax + gauss_rnd]
    fadd    qword ptr [eax + gauss_rnd+8]
    fadd    qword ptr [eax + gauss_rnd+16]
    fadd    qword ptr [eax + gauss_rnd+24]
    fld1
    fadd    st, st
    fsubp   st(1),st
    fmul    qword ptr [eax + dofCoef]
    fmul    qword ptr [z]
    fstp    qword ptr [dr]
    call    AsmRandExt
    mov     edx, [eax + gauss_N]
    fst     qword ptr [eax + gauss_rnd + edx*8]
    inc     edx
    and     edx,$03
    mov     [eax + gauss_N], edx
    fadd    st, st
    fldpi
    fmulp
    fsincos
    fstp    qword ptr [dcos]
    fstp    qword ptr [dsin]
  end;
{$else}
  t := random * 2 * pi;
  SinCos(t, dsin, dcos);
  dr := random * dofCoef * z;
{
  asm
    fld     qword ptr [z]
    fmul    st, st
    fmul    qword ptr [eax + dofCoef]
    fldpi
    fadd    st, st
    call    AsmRandExt
    fmulp
    fsincos
    fstp    qword ptr [dcos]
    fstp    qword ptr [dsin]
    call    AsmRandExt
    fmulp
    fstp    qword ptr [dr]
  end;
}
{$endif}

  pPoint^.x := (pPoint^.x + dr*dcos) / zr;
  pPoint^.y := (y + dr*dsin) / zr;
  //pPoint^.z := pPoint^.z - CameraZpos;
end;

procedure TControlPoint.ProjectPitchYawDOF(pPoint: PCPPoint);
var
  x, y, z, zr, dr: double;
  dsin, dcos: double;
  t : double;
begin
  z := pPoint^.z - CameraZpos;
  x := cosYaw * pPoint^.x - sinYaw * pPoint^.y;
  y := CameraMatrix[0,1] * pPoint^.x + CameraMatrix[1,1] * pPoint^.y - sinPitch * z;
  z := CameraMatrix[0,2] * pPoint^.x + CameraMatrix[1,2] * pPoint^.y + cosPitch * z;
  zr := 1 - cameraPersp * z;

{$ifdef GAUSSIAN_DOF}
  asm
    fld     qword ptr [eax + gauss_rnd]
    fadd    qword ptr [eax + gauss_rnd+8]
    fadd    qword ptr [eax + gauss_rnd+16]
    fadd    qword ptr [eax + gauss_rnd+24]
    fld1
    fadd    st, st
    fsubp   st(1),st
    fmul    qword ptr [eax + dofCoef]
    fmul    qword ptr [z]
    fstp    qword ptr [dr]
    call    AsmRandExt
    mov     edx, [eax + gauss_N]
    fst     qword ptr [eax + gauss_rnd + edx*8]
    inc     edx
    and     edx,$03
    mov     [eax + gauss_N], edx
    fadd    st, st
    fldpi
    fmulp
    fsincos
    fstp    qword ptr [dcos]
    fstp    qword ptr [dsin]
  end;
{$else}
  t := random * 2 * pi;
  SinCos(t, dsin, dcos);
  dr := random * dofCoef * z;
{
  asm
    fld     qword ptr [z]
    fmul    st, st
    fmul    qword ptr [eax + dofCoef]
    fldpi
    fadd    st, st
    call    AsmRandExt
    fmulp
    fsincos
    fstp    qword ptr [dcos]
    fstp    qword ptr [dsin]
    call    AsmRandExt
    fmulp
    fstp    qword ptr [dr]
  end;
}
{$endif}

  pPoint^.x := (x + dr*dcos) / zr;
  pPoint^.y := (y + dr*dsin) / zr;
 // pPoint^.z := pPoint^.z - CameraZpos;
end;

procedure TControlPoint.ProjectPitchRoll(pPoint: PCPPoint);
var  // AV: added Y-rotation
  x, y, z, zr: double;
begin
  z := pPoint^.z - CameraZpos;
  x := cosRoll * pPoint^.x + CameraMatrix[1,0] * pPoint^.y + CameraMatrix[1,0] * z;
  y := cosPitch * pPoint^.y - sinPitch * z;
  zr := 1 - cameraPersp *
    (sinRoll * pPoint^.x + CameraMatrix[1,2] * pPoint^.y + CameraMatrix[2,2] * z);

  pPoint^.x := x / zr;
  pPoint^.y := y / zr;
end;

procedure TControlPoint.ProjectPitchRollYaw(pPoint: PCPPoint);
var  // AV: added Y-rotation
  x, y, z, zr: double;
begin
  z := pPoint^.z - CameraZpos;
  x := CameraMatrix[0,0] * pPoint^.x + CameraMatrix[1,0] * pPoint^.y + CameraMatrix[2,0]  * z;
  y := CameraMatrix[0,1] * pPoint^.x + CameraMatrix[1,1] * pPoint^.y + CameraMatrix[2,1]  * z;
  zr := 1 - cameraPersp *
    (CameraMatrix[0,2]* pPoint^.x + CameraMatrix[1,2]*pPoint^.y + CameraMatrix[2,2] * z);

  pPoint^.x := x / zr;
  pPoint^.y := y / zr;
end;

procedure TControlPoint.ProjectPitchRollYawDOF(pPoint: PCPPoint);
var   // AV: added Y-rotation
  x, y, z, zr, dr: double;
  dsin, dcos: double;
  t : double;
begin
  z := pPoint^.z - CameraZpos;
  x := CameraMatrix[0,0] * pPoint^.x + CameraMatrix[1,0] * pPoint^.y + CameraMatrix[2,0]  * z;
  y := CameraMatrix[0,1] * pPoint^.x + CameraMatrix[1,1] * pPoint^.y + CameraMatrix[2,1]  * z;
  z := CameraMatrix[0,2]* pPoint^.x + CameraMatrix[1,2]*pPoint^.y + CameraMatrix[2,2] * z;
  zr := 1 - cameraPersp * z;

  t := random * 2 * pi;
  SinCos(t, dsin, dcos);
  dr := random * dofCoef * z;

  pPoint^.x := (x + dr * dcos) / zr;
  pPoint^.y := (y + dr * dsin) / zr;
end;

function TControlPoint.BlowsUp(NrPoints: integer): boolean;
const
  limit = 1E10;  // AV
var
  i, n: Integer;
  px, py: double;
  minx, maxx, miny, maxy: double;
  Points: TPointsXYArray;
  CurrentPoint: PXYPoint;

  xf: TXForm;
begin
  Result := false;

  n := min(SUB_BATCH_SIZE, NrPoints);
  SetLength(Points, n);

  px := 2 * random - 1;
  py := 2 * random - 1;

  Prepare;

  try
    xf := xform[random(NumXForms)];
    for i := 0 to FUSE do begin
      xf := xf.PropTable[Random(PROP_TABLE_SIZE)];
      xf.NextPointXY(px,py);
    end;

    CurrentPoint := @Points[0];
    for i := 0 to n-1 do begin
      xf := xf.PropTable[Random(PROP_TABLE_SIZE)];
      xf.NextPointXY(px,py);
      CurrentPoint.X := px;
      CurrentPoint.Y := py;
      Inc(CurrentPoint);
      // random CPs don't use finalXform...
    end;

    // It is possible that the transformation will grow very large but remain below the overflow line
    minx := limit; //1E10;
    maxx := -limit; //-1E10;
    miny := limit; //1E10;
    maxy := -limit; //-1E10;
    for i := 0 to n-1 do begin
      // AV: rewrote the block to fix the crash
      if InRange(Points[i].x, -limit, limit) then
      begin
        minx := min(minx, Points[i].x);
        maxx := max(maxx, Points[i].x);
      end
      else
      begin
        minx := -limit;
        maxx := limit;
      end;
      if InRange(Points[i].y, -limit, limit) then
      begin
        miny := min(miny, Points[i].y);
        maxy := max(maxy, Points[i].y);
      end
      else
      begin
        miny := -limit;
        maxy := limit;
      end;
    end;

    if ((Maxx - MinX) > 1000) or ((Maxy - Miny) > 1000) then
      Result := True;

  except
    on EMathError do begin
      Result := True;
      Exit;
    end;
  end;

end;

procedure TControlPoint.ParseString(aString: string);
var
  ParseValues: TStringList;
  ParsePos: integer;
  CurrentToken: string;
  CurrentXForm: integer;
  i: integer;
  OldDecimalSperator: Char;
  v: double;
begin
  ParseValues := TStringList.Create;
  ParseValues.CommaText := AString;

  OldDecimalSperator := FormatSettings.DecimalSeparator;
  FormatSettings.DecimalSeparator := '.';

  CurrentXForm := 0;

  ParsePos := 0;
  while (ParsePos < ParseValues.Count) do begin
    CurrentToken := ParseValues[ParsePos];
    if AnsiCompareText(CurrentToken, 'xform') = 0 then begin
      Inc(ParsePos);
      CurrentXForm := StrToInt(ParseValues[ParsePos]);
    (*end else if AnsiCompareText(CurrentToken, 'plugins') = 0 then begin
      used_plugins.Clear;
      i := 0;
      while true do begin
        if (ParsePos + 1) >= ParseValues.Count then
          break;

        Inc(ParsePos);
        used_plugins.Add(ParseValues[ParsePos]);
        Inc(i);
      end; *)
    end else if AnsiCompareText(CurrentToken, 'finalxformenabled') = 0 then begin
      Inc(ParsePos);
      finalxformenabled := StrToInt(ParseValues[ParsePos]) <> 0;
    end else if AnsiCompareText(CurrentToken, 'soloxform') = 0 then begin
      Inc(ParsePos);
      soloxform := StrToInt(ParseValues[ParsePos]);
    end else if AnsiCompareText(CurrentToken, 'time') = 0 then begin
      Inc(ParsePos);
      time := StrToFloat(ParseValues[ParsePos]);
    end else if AnsiCompareText(CurrentToken, 'brightness') = 0 then begin
      Inc(ParsePos);
      brightness := StrToFloat(ParseValues[ParsePos]) / BRIGHT_ADJUST;
    end else if AnsiCompareText(CurrentToken, 'zoom') = 0 then begin // mt
      Inc(ParsePos); // mt
      zoom := StrToFloat(ParseValues[ParsePos]); // mt
    end else if AnsiCompareText(CurrentToken, 'angle') = 0 then begin
      Inc(ParsePos);
      FAngle := StrToFloat(ParseValues[ParsePos]);
// 3d camera stuff
    end else if AnsiCompareText(CurrentToken, 'cam_pitch') = 0 then begin
      Inc(ParsePos);
      cameraPitch := StrToFloat(ParseValues[ParsePos]);
    end else if AnsiCompareText(CurrentToken, 'cam_yaw') = 0 then begin
      Inc(ParsePos);
      cameraYaw := StrToFloat(ParseValues[ParsePos]);
    end else if AnsiCompareText(CurrentToken, 'cam_roll') = 0 then begin
      Inc(ParsePos); // AV: add third angle
      cameraRoll := StrToFloat(ParseValues[ParsePos]);
    end else if AnsiCompareText(CurrentToken, 'cam_persp') = 0 then begin
      Inc(ParsePos);
      cameraPersp := StrToFloat(ParseValues[ParsePos]);
    end else if AnsiCompareText(CurrentToken, 'cam_zpos') = 0 then begin
      Inc(ParsePos);
      cameraZpos := StrToFloat(ParseValues[ParsePos]);
    end else if AnsiCompareText(CurrentToken, 'cam_dof') = 0 then begin
      Inc(ParsePos);
      cameraDOF := abs(StrToFloat(ParseValues[ParsePos]));
// end 3d
    end else if AnsiCompareText(CurrentToken, 'contrast') = 0 then begin
      Inc(ParsePos);
      contrast := StrToFloat(ParseValues[ParsePos]);
    end else if AnsiCompareText(CurrentToken, 'gamma') = 0 then begin
      Inc(ParsePos);
      gamma := StrToFloat(ParseValues[ParsePos]);
    end else if AnsiCompareText(CurrentToken, 'vibrancy') = 0 then begin
      Inc(ParsePos);
      vibrancy := StrToFloat(ParseValues[ParsePos]);
    end else if AnsiCompareText(CurrentToken, 'gamma_threshold') = 0 then begin
      Inc(ParsePos);
      gamma_threshold := StrToFloat(ParseValues[ParsePos]);
    end else if AnsiCompareText(CurrentToken, 'hue_rotation') = 0 then begin
      Inc(ParsePos);
      hue_rotation := StrToFloat(ParseValues[ParsePos]);
    end else if AnsiCompareText(CurrentToken, 'zoom') = 0 then begin
      Inc(ParsePos);
      zoom := StrToFloat(ParseValues[ParsePos]);
    end else if AnsiCompareText(CurrentToken, 'imagesize') = 0 then begin
      Inc(ParsePos);
      Width := StrToInt(ParseValues[ParsePos]);
      Inc(ParsePos);
      Height := StrToInt(ParseValues[ParsePos]);
    end else if AnsiCompareText(CurrentToken, 'image_size') = 0 then begin
      Inc(ParsePos);
      Width := StrToInt(ParseValues[ParsePos]);
      Inc(ParsePos);
      Height := StrToInt(ParseValues[ParsePos]);
    end else if AnsiCompareText(CurrentToken, 'center') = 0 then begin
      Inc(ParsePos);
      center[0] := StrToFloat(ParseValues[ParsePos]);
      Inc(ParsePos);
      center[1] := StrToFloat(ParseValues[ParsePos]);
    end else if AnsiCompareText(CurrentToken, 'background') = 0 then begin
      Inc(ParsePos);
      // Trap conversion errors for older parameters
      try
        background[0] := StrToInt(ParseValues[ParsePos]);
      except on EConvertError do
          background[0] := 0;
      end;
      Inc(ParsePos);
      try
        background[1] := StrToInt(ParseValues[ParsePos]);
      except on EConvertError do
          background[1] := 0;
      end;
      Inc(ParsePos);
      try
        background[2] := StrToInt(ParseValues[ParsePos]);
      except on EConvertError do
          background[2] := 0;
      end;
    end else if AnsiCompareText(CurrentToken, 'curves') = 0 then begin
      for i := 0 to 3 do
      begin
        Inc(ParsePos);curvePoints[i][0].x := StrToFloat(ParseValues[ParsePos]);
        Inc(ParsePos);curvePoints[i][0].y := StrToFloat(ParseValues[ParsePos]);
        Inc(ParsePos);curveWeights[i][0] := StrToFloat(ParseValues[ParsePos]);
        Inc(ParsePos);curvePoints[i][1].x := StrToFloat(ParseValues[ParsePos]);
        Inc(ParsePos);curvePoints[i][1].y := StrToFloat(ParseValues[ParsePos]);
        Inc(ParsePos);curveWeights[i][1] := StrToFloat(ParseValues[ParsePos]);
        Inc(ParsePos);curvePoints[i][2].x := StrToFloat(ParseValues[ParsePos]);
        Inc(ParsePos);curvePoints[i][2].y := StrToFloat(ParseValues[ParsePos]);
        Inc(ParsePos);curveWeights[i][2] := StrToFloat(ParseValues[ParsePos]);
        Inc(ParsePos);curvePoints[i][3].x := StrToFloat(ParseValues[ParsePos]);
        Inc(ParsePos);curvePoints[i][3].y := StrToFloat(ParseValues[ParsePos]);
        Inc(ParsePos);curveWeights[i][3] := StrToFloat(ParseValues[ParsePos]);
      end;
    end {else if AnsiCompareText(CurrentToken, 'pulse') = 0 then begin
      Inc(ParsePos);
      pulse[0, 0] := StrToFloat(ParseValues[ParsePos]);
      Inc(ParsePos);
      pulse[0, 1] := StrToFloat(ParseValues[ParsePos]);
      Inc(ParsePos);
      pulse[1, 0] := StrToFloat(ParseValues[ParsePos]);
      Inc(ParsePos);
      pulse[1, 1] := StrToFloat(ParseValues[ParsePos]);
    end else if AnsiCompareText(CurrentToken, 'wiggle') = 0 then begin
      Inc(ParsePos);
      wiggle[0, 0] := StrToFloat(ParseValues[ParsePos]);
      Inc(ParsePos);
      wiggle[0, 1] := StrToFloat(ParseValues[ParsePos]);
      Inc(ParsePos);
      wiggle[1, 0] := StrToFloat(ParseValues[ParsePos]);
      Inc(ParsePos);
      wiggle[1, 1] := StrToFloat(ParseValues[ParsePos]);
    end} else if AnsiCompareText(CurrentToken, 'pixels_per_unit') = 0 then begin
      Inc(ParsePos);
      pixels_per_unit := StrToFloat(ParseValues[ParsePos]);
    end else if AnsiCompareText(CurrentToken, 'spatial_filter_radius') = 0 then begin
      Inc(ParsePos);
      spatial_filter_radius := StrToFloat(ParseValues[ParsePos]);
    end else if AnsiCompareText(CurrentToken, 'spatial_oversample') = 0 then begin
      Inc(ParsePos);
      spatial_oversample := StrToInt(ParseValues[ParsePos]);
    end else if AnsiCompareText(CurrentToken, 'sample_density') = 0 then begin
      Inc(ParsePos);
      sample_density := StrToFloat(ParseValues[ParsePos]);
    end else if AnsiCompareText(CurrentToken, 'nbatches') = 0 then begin
      Inc(ParsePos);
      nbatches := StrToInt(ParseValues[ParsePos]);
    end else if AnsiCompareText(CurrentToken, 'white_level') = 0 then begin
      Inc(ParsePos);
      white_level := StrToInt(ParseValues[ParsePos]);
    end else if AnsiCompareText(CurrentToken, 'cmap') = 0 then begin
      Inc(ParsePos);
      cmapindex := StrToInt(ParseValues[ParsePos]);
    end {else if AnsiCompareText(CurrentToken, 'cmap_inter') = 0 then begin
      Inc(ParsePos);
      cmap_inter := StrToInt(ParseValues[ParsePos]);
    end else if AnsiCompareText(CurrentToken, 'palette') = 0 then begin
//      Inc(ParsePos);
//      cmapindex := StrToInt(ParseValues[ParsePos]);
      OutputDebugString(Pchar('NYI import Palette'));
    end} else if AnsiCompareText(CurrentToken, 'density') = 0 then begin
      Inc(ParsePos);
      xform[CurrentXForm].Density := StrToFloat(ParseValues[ParsePos]);
    end else if AnsiCompareText(CurrentToken, 'color') = 0 then begin
      Inc(ParsePos);
      xform[CurrentXForm].color := StrToFloat(ParseValues[ParsePos]);
    end else if AnsiCompareText(CurrentToken, 'symmetry') = 0 then begin
      Inc(ParsePos);
      xform[CurrentXForm].symmetry := StrToFloat(ParseValues[ParsePos]);
    end else if AnsiCompareText(CurrentToken, 'coefs') = 0 then begin
      Inc(ParsePos);
      xform[CurrentXForm].c[0, 0] := StrToFloat(ParseValues[ParsePos]);
      Inc(ParsePos);
      xform[CurrentXForm].c[0, 1] := StrToFloat(ParseValues[ParsePos]);
      Inc(ParsePos);
      xform[CurrentXForm].c[1, 0] := StrToFloat(ParseValues[ParsePos]);
      Inc(ParsePos);
      xform[CurrentXForm].c[1, 1] := StrToFloat(ParseValues[ParsePos]);
      Inc(ParsePos);
      xform[CurrentXForm].c[2, 0] := StrToFloat(ParseValues[ParsePos]);
      Inc(ParsePos);
      xform[CurrentXForm].c[2, 1] := StrToFloat(ParseValues[ParsePos]);
    end else if AnsiCompareText(CurrentToken, 'post') = 0 then begin
      Inc(ParsePos);
      xform[CurrentXForm].p[0, 0] := StrToFloat(ParseValues[ParsePos]);
      Inc(ParsePos);
      xform[CurrentXForm].p[0, 1] := StrToFloat(ParseValues[ParsePos]);
      Inc(ParsePos);
      xform[CurrentXForm].p[1, 0] := StrToFloat(ParseValues[ParsePos]);
      Inc(ParsePos);
      xform[CurrentXForm].p[1, 1] := StrToFloat(ParseValues[ParsePos]);
      Inc(ParsePos);
      xform[CurrentXForm].p[2, 0] := StrToFloat(ParseValues[ParsePos]);
      Inc(ParsePos);
      xform[CurrentXForm].p[2, 1] := StrToFloat(ParseValues[ParsePos]);
    end else if AnsiCompareText(CurrentToken, 'postxswap') = 0 then begin
      Inc(ParsePos);
      xform[CurrentXForm].postXswap := (ParseValues[ParsePos] = '1');
    end else if AnsiCompareText(CurrentToken, 'autozscale') = 0 then begin
      Inc(ParsePos);
      xform[CurrentXForm].autoZscale := (ParseValues[ParsePos] = '1');
	end
  ////////////////////
	else if AnsiCompareText(CurrentToken, 'var_order') = 0 then begin
      for i := 0 to NRVAR - 1 do 
        xform[CurrentXForm].ifs[i] := Varnames(i);

		i := 0;
    while true do begin
			if (ParsePos + 1) >= ParseValues.Count then
				break;
			if CharInSet(ParseValues[ParsePos + 1][1], ['a'..'z', 'A'..'Z', '.']) then
				break;

			Inc(ParsePos);
			xform[CurrentXForm].ifs[i] := Varnames(StrToInt(ParseValues[ParsePos]));
			Inc(i);
    end;
	end
  //////////////////////////
	else if AnsiCompareText(CurrentToken, 'vars') = 0 then begin
      for i := 0 to NRVAR - 1 do begin
        xform[CurrentXForm].SetVariation(i, 0.0);
      end;

      i := 0;
      while true do begin
        if (ParsePos + 1) >= ParseValues.Count then
          break;
        if CharInSet(ParseValues[ParsePos + 1][1], ['a'..'z', 'A'..'Z']) then
          break;

        Inc(ParsePos);
        xform[CurrentXForm].SetVariation(i, StrToFloat(ParseValues[ParsePos]));
        Inc(i);
      end;
    end else if AnsiCompareText(CurrentToken, 'variables') = 0 then begin
{
      v := 0;
      for i:= 0 to GetNrVariableNames-1 do begin
        xform[CurrentXForm].SetVariable(GetVariableNameAt(i), v);
      end;
}
      for i:= 0 to GetNrVariableNames-1 do begin
        xform[CurrentXForm].ResetVariable(GetVariableNameAt(i));
      end;

      i := 0;
      while true do begin
        if (ParsePos + 1) >= ParseValues.Count then
          break;
        if CharInSet(ParseValues[ParsePos + 1][1], ['a'..'z', 'A'..'Z']) then
          break;

        Inc(ParsePos);
        v := StrToFloat(ParseValues[ParsePos]);
        xform[CurrentXForm].SetVariable(GetVariableNameAt(i), v);
        Inc(i);
      end;

    end else if AnsiCompareText(CurrentToken, 'chaos') = 0 then begin
      i := 0;
      while true do begin
        if (ParsePos + 1) >= ParseValues.Count then
          break;
        if CharInSet(ParseValues[ParsePos + 1][1], ['a'..'z', 'A'..'Z']) then
          break;

        Inc(ParsePos);
        v := StrToFloat(ParseValues[ParsePos]);
        xform[CurrentXForm].modWeights[i] := v;
        Inc(i);
      end;

    end else if AnsiCompareText(CurrentToken, 'plotmode') = 0 then begin
      Inc(ParsePos);
      if((StrToInt(ParseValues[ParsePos]) = 1)) then
        xform[CurrentXForm].transOpacity := 0;
    end else if AnsiCompareText(CurrentToken, 'opacity') = 0 then begin
      Inc(ParsePos);
      xform[CurrentXForm].transOpacity := (StrToFloat(ParseValues[ParsePos]));
    end else if AnsiCompareText(CurrentToken, 'var_color') = 0 then begin
      Inc(ParsePos);
      xform[CurrentXForm].pluginColor := (StrToFloat(ParseValues[ParsePos]));
    end else begin
      OutputDebugString(Pchar('Unknown Token: ' + CurrentToken));
    end;

    Inc(ParsePos);
  end;
  GetCmap(cmapindex, hue_rotation, Cmap);

  ParseValues.Free;

  FormatSettings.DecimalSeparator := OldDecimalSperator;
end;

(*
procedure TControlPoint.SetVariation(vari: TVariation);
var
  i, j, v: integer;
  rv: integer;
  VarPossible: boolean;
begin
  FillVarDisturb;
  VarPossible := false;
  for j := 0 to NRVAR - 1 do begin
    VarPossible := VarPossible or Variations[j];
  end;

  if VarPossible then begin
    repeat
      rv := var_distrib[random(Length(var_distrib))];
    until Variations[rv];
  end else begin
    rv := 0;
  end;

  for i := 0 to NXFORMS - 1 do begin
    for j := 0 to NRVAR - 1 do begin
      xform[i].SetVariation(j, 0.0);
    end;

    if vari = vRandom then
    begin
      if rv < 0 then
      begin
        if VarPossible then begin
          repeat
            v := Mixed_var_distrib[random(Length(mixed_var_distrib))];
          until Variations[v]; // Use only Variations set in options
        end else begin
          v := 0;
        end;
        xform[i].SetVariation(v, 1.0);
      end
      else
        xform[i].SetVariation(rv, 1.0);
    end
    else
      xform[i].SetVariation(integer(vari), 1);
  end;
end;

procedure TControlPoint.RandomCP(min: integer = 2; max: integer = NXFORMS; calc: boolean = true);
var
  nrXforms: integer;
  i, j: integer;
  v, rv: integer;
  VarPossible: boolean;
begin
//hue_rotation := random;
  hue_rotation := 1;
  cmapindex := RANDOMCMAP;
  GetCmap(cmapindex, hue_rotation, cmap);
  time := 0.0;

  nrXforms := random(Max - (Min - 1)) + Min;

  FillVarDisturb;
  VarPossible := false;
  for j := 0 to NRVAR - 1 do begin
    VarPossible := VarPossible or Variations[j];
  end;

  if VarPossible then begin
    repeat
      rv := var_distrib[random(Length(var_distrib))];
    until Variations[rv];
  end else begin
    rv := 0;
  end;

  for i := 0 to NXFORMS - 1 do begin
    xform[i].density := 0;
  end;

  for i := 0 to nrXforms - 1 do begin
    xform[i].density := 1.0 / nrXforms;
    xform[i].color := i / (nrXforms - 1);

    xform[i].c[0][0] := 2 * random - 1;
    xform[i].c[0][1] := 2 * random - 1;
    xform[i].c[1][0] := 2 * random - 1;
    xform[i].c[1][1] := 2 * random - 1;
    xform[i].c[2][0] := 4 * random - 2;
    xform[i].c[2][1] := 4 * random - 2;

    for j := 0 to NRVAR - 1 do begin
      xform[i].SetVariation(j, 0);
    end;

    if rv < 0 then begin
      if VarPossible then begin
        repeat
          v := Mixed_var_distrib[random(Length(mixed_var_distrib))];
        until Variations[v]; // use only variations set in options
      end else begin
        v := 0;
      end;

      xform[i].SetVariation(v, 1);
    end else
      xform[i].SetVariation(rv, 1);

  end;
  if calc then CalcBoundbox;
end;

procedure TControlPoint.RandomCP1;
var
  i, j: Integer;
begin
  RandomCP;
  for i := 0 to NXFORMS - 1 do begin
    for j := 1 to NRVAR - 1 do
      xform[i].SetVariation(j, 0);
    xform[i].SetVariation(0, 1);
  end;

  CalcBoundbox;
end;
*)

procedure TControlPoint.PreCalcBounds(var maxx, minx, maxy, miny: double);
const
  limit = 50;  // AV
var
  Points: TPointsArray; // AV: fixed - was TPointsXYArray;
  i, j: integer;
  deltax, deltay: double;
  cntminx, cntmaxx: integer;
  cntminy, cntmaxy: integer;
  LimitOutSidePoints: integer;
begin
    SetLength(Points, SUB_BATCH_SIZE);
    IterateXYC(SUB_BATCH_SIZE, points);  // AV: fixed - was IterateXY

    LimitOutSidePoints := Round(0.05 * SUB_BATCH_SIZE);
    // who've been coding this nonsence?!
    minx := limit; // 1e99;
    maxx := -limit; // -1e99;
    miny := limit; // 1e99;
    maxy := -limit;  // -1e99;
    for i := 0 to SUB_BATCH_SIZE - 1 do begin
    //if Points[i].x > 1e200 then continue;
    // AV: rewrote the block to fix the crash
      if InRange(Points[i].x, -limit, limit) then
      begin
        minx := min(minx, Points[i].x);
        maxx := max(maxx, Points[i].x);
      end
      else
      begin
        minx := -limit;
        maxx := limit;
      end;
      if InRange(Points[i].y, -limit, limit) then
      begin
        miny := min(miny, Points[i].y);
        maxy := max(maxy, Points[i].y);
      end
      else
      begin
        miny := -limit;
        maxy := limit;
      end;
    end;

    deltax := (maxx - minx) * 0.25;
    maxx := (maxx + minx) / 2;
    minx := maxx;

    deltay := (maxy - miny) * 0.25;
    maxy := (maxy + miny) / 2;
    miny := maxy;

    for j := 0 to 10 do begin
      cntminx := 0;
      cntmaxx := 0;
      cntminy := 0;
      cntmaxy := 0;
      for i := 0 to SUB_BATCH_SIZE - 1 do begin
       // AV: rewrote the block to fix the crash
        if InRange(Points[i].x, -limit, limit) then
        begin
          if (Points[i].x < minx) then Inc(cntminx);
          if (Points[i].x > maxx) then Inc(cntmaxx);
        end;
        if InRange(Points[i].y, -limit, limit) then
        begin
          if (Points[i].y < miny) then Inc(cntminy);
          if (Points[i].y > maxy) then Inc(cntmaxy);
        end;
      end;

      if (cntMinx < LimitOutSidePoints) then begin
        minx := minx + deltax;
      end else begin
        minx := minx - deltax;
      end;

      if (cntMaxx < LimitOutSidePoints) then begin
        maxx := maxx - deltax;
      end else begin
        maxx := maxx + deltax;
      end;

      deltax := deltax / 2;

      if (cntMiny < LimitOutSidePoints) then begin
        miny := miny + deltay;
      end else begin
        miny := miny - deltay;
      end;

      if (cntMaxy < LimitOutSidePoints) then begin
        maxy := maxy - deltay;
      end else begin
        maxy := maxy + deltay;
      end;

      deltay := deltay / 2;
    end;

    if (abs(maxx - minx) > 1000) or
       (abs(maxy - miny) > 1000) then
      raise EMathError.Create('Flame area too large');

    center[0] := (minx + maxx) / 2;
    center[1] := (miny + maxy) / 2;
    if (abs(maxx - minx) > 0.001) and (abs(maxy - miny) > 0.001) then
      pixels_per_unit := 0.65 * Min(Width / abs(maxx - minx),
                        Height / abs(maxy - miny))
    else
      pixels_per_unit := 10;
end;

procedure TControlPoint.CalcBoundbox;
var
  minx, maxx, miny, maxy: double;
begin
  try
    Prepare;
    PreCalcBounds(maxx, minx, maxy, miny);  // AV
  except on EMathError do
    begin	// default
      center[0] := 0;
      center[1] := 0;
      pixels_per_unit := 10;
    end;
  end;
end;

function CalcUPRMagn(const cp: TControlPoint): double;
var
  minx, maxx, miny, maxy: double;
  xLength, yLength: double;
begin
  try
    cp.PreCalcBounds(maxx, minx, maxy, miny);  // AV

   // Calculate magn for UPRs
    xLength := maxx - minx;
    yLength := maxy - miny;
    if xLength >= yLength then
    begin
      result := 1 / xLength * 2;
    end
    else
    begin
      result := 1 / yLength * 2;
    end;

  except on E: EMathError do
    begin// default
      cp.center[0] := 0;
      cp.center[1] := 0;
      cp.pixels_per_unit := 10;
      raise Exception.Create('CalcUPRMagn: ' + e.Message);
    end;
  end;
end;

procedure PrepareToInterpolation(var SourceCp, TargetCp: TControlPoint); // AV
var
    i, j: integer;
    t, ns, nt, maxt: smallint;
    vs, vt: double;
    vn: string;
begin
  { reset linear variation on (temporary) invisible xforms }
  ns := SourceCp.NumXForms;
  nt := TargetCp.NumXForms;
  if ns > nt then begin
    for j := nt + 1 to ns - 1 do
      TargetCp.xform[j].SetVariation(0, 0);
    maxt := ns;  // max number of actually used xforms
  end
  else if nt > ns then begin
    for j := ns + 1 to nt - 1 do
      SourceCp.xform[j].SetVariation(0, 0);
    maxt := nt;
  end
  else // if nt = ns then
    maxt := nt;

  { adjust (temporary unused) variable values }
  for t := 0 to maxt do
    for i := 0 to GetNrVariableNames - 1 do
    begin
      j := GetVariationIndexFromVariableNameIndex(i);
      vs := SourceCp.xform[t].GetVariation(j);
      vt := TargetCp.xform[t].GetVariation(j);
      if (vs <> 0) and (vt = 0) then
      begin
        vn := GetVariableNameAt(i);
        SourceCp.xform[t].GetVariable(vn, vs);
        TargetCp.xform[t].SetVariable(vn, vs);
      end
      else if (vt <> 0) and (vs = 0) then
      begin
        vn := GetVariableNameAt(i);
        TargetCp.xform[t].GetVariable(vn, vt);
        SourceCp.xform[t].SetVariable(vn, vt);
      end;
    end;

  { adjust final transforms }
  if (SourceCp.finalXformEnabled or TargetCp.finalXformEnabled)
  and (nt <> ns) then  // otherwise it's already done
  begin
    for i := 0 to GetNrVariableNames - 1 do
    begin
      j := GetVariationIndexFromVariableNameIndex(i);
      vs := SourceCp.xform[ns].GetVariation(j);
      vt := TargetCp.xform[nt].GetVariation(j);
      if (vs <> 0) and (vt = 0) then
      begin
        vn := GetVariableNameAt(i);
        SourceCp.xform[ns].GetVariable(vn, vs);
        TargetCp.xform[nt].SetVariable(vn, vs);
      end
      else if (vt <> 0) and (vs = 0) then
      begin
        vn := GetVariableNameAt(i);
        TargetCp.xform[nt].GetVariable(vn, vt);
        SourceCp.xform[ns].SetVariable(vn, vt);
      end;
    end;
  end;

  { AV: improve visibility of mid-frames }
  SourceCp.NormalizeProbabilities;
  TargetCp.NormalizeProbabilities;
end;

// AV: this one from mine is used for scripting, the original method was left for mutants only
procedure TControlPoint.InterpolateAll(cp1, cp2: TControlPoint; Tm: double; it: integer);
var
  c0, c1, k: double;
  i, j: integer;
  numTX, numTX1, numTX2: word;
  r, s, t: array[0..2] of double;
  v1, v2, t1, t2: double;
  f1, f2: boolean;

  function AngOpt(ang0, ang1: double): double;
  var delta: double;
  begin
    delta := ang1 - ang0;
    if (delta > pi) then Result := ang1 - 2 * pi
    else if (delta < -pi)then Result := ang1 + 2 * pi
    else Result := ang1;
  end;

begin
  if (it > 3) or (it < 0) then exit;
  
  t1 := cp1.time;
  t2 := cp2.time;
  
  if (tm = t1) then
  begin
	  self.Copy(cp1);
	  exit;
  end
  else if (tm = t2) then
  begin
	  self.Copy(cp2);
    width := cp1.width;
    height := cp1.height;
	  exit;
  end;  
  
  if (t2 - t1) > 1E-6 then
  begin
    if (it < 2) then begin  // 0 or 1: cosine interpolation
      k := (t2 - tm) / (t2 - t1);
      c0 := (1 - cos(pi * k)) *  0.5;
    end
    else  // 2 or 3: linear interpolation
      c0 := (t2 - tm) / (t2 - t1);
    c1 := 1 - c0;
  end else begin
    c0 := 1;
    c1 := 0;
  end;

  self.time := Tm;
  
  if Odd(it) then  // 1 or 3: hsv interpolation
      for i := 0 to 255 do
      begin
      r[0] := cp1.cmap[i][0] / 255;
      r[1] := cp1.cmap[i][1] / 255;
      r[2] := cp1.cmap[i][2] / 255;
      rgb2hsv(r, s);
      r[0] := cp2.cmap[i][0] / 255;
      r[1] := cp2.cmap[i][1] / 255;
      r[2] := cp2.cmap[i][2] / 255;
      rgb2hsv(r, t);
      t[0] := c0 * s[0] + c1 * t[0];
      t[1] := c0 * s[1] + c1 * t[1];
      t[2] := c0 * s[2] + c1 * t[2];
      hsv2rgb(t, r);
      cmap[i][0] := Round(255 * r[0]);
      cmap[i][1] := Round(255 * r[1]);
      cmap[i][2] := Round(255 * r[2]);
    end
  else  // 0 or 2: rgb interpolation
  for i := 0 to 255 do
  begin
    cmap[i][0] := Round(c0 * cp1.cmap[i][0] + c1 * cp2.cmap[i][0]);
    cmap[i][1] := Round(c0 * cp1.cmap[i][1] + c1 * cp2.cmap[i][1]);
    cmap[i][2] := Round(c0 * cp1.cmap[i][2] + c1 * cp2.cmap[i][2]);
  end;

  self.cmapindex := -1;

  Fbrightness := c0 * cp1.Fbrightness + c1 * cp2.Fbrightness;
  contrast := c0 * cp1.contrast + c1 * cp2.contrast;
  gamma := c0 * cp1.gamma + c1 * cp2.gamma;
  vibrancy := c0 * cp1.vibrancy + c1 * cp2.vibrancy;
  gamma_threshold := c0 * cp1.gamma_threshold + c1 * cp2.gamma_threshold;
  width := cp1.width;
  height := cp1.height;
  spatial_oversample := Round(c0 * cp1.spatial_oversample + c1 * cp2.spatial_oversample);
  center[0] := c0 * cp1.center[0] + c1 * cp2.center[0];
  center[1] := c0 * cp1.center[1] + c1 * cp2.center[1];
  FAngle := c0 * cp1.FAngle + c1 * AngOpt(cp1.FAngle, cp2.FAngle);
  pixels_per_unit := c0 * cp1.pixels_per_unit + c1 * cp2.pixels_per_unit;
  background[0] := Round(c0 * cp1.background[0] + c1 * cp2.background[0]);
  background[1] := Round(c0 * cp1.background[1] + c1 * cp2.background[1]);
  background[2] := Round(c0 * cp1.background[2] + c1 * cp2.background[2]);
  spatial_filter_radius := c0 * cp1.spatial_filter_radius + c1 * cp2.spatial_filter_radius;
  sample_density := c0 * cp1.sample_density + c1 * cp2.sample_density;
  zoom := c0 * cp1.zoom + c1 * cp2.zoom;
  nbatches := Round(c0 * cp1.nbatches + c1 * cp2.nbatches);
  // white_level := Round(c0 * cp1.white_level + c1 * cp2.white_level);

  { AV: global 3D-parameters}
  cameraPitch := c0 * cp1.cameraPitch + c1 * AngOpt(cp1.cameraPitch, cp2.cameraPitch);
  cameraYaw := c0 * cp1.cameraYaw + c1 * AngOpt(cp1.cameraYaw, cp2.cameraYaw);
  cameraRoll := c0 * cp1.cameraRoll + c1 * AngOpt(cp1.cameraRoll, cp2.cameraRoll);
  cameraPersp := c0 * cp1.cameraPersp + c1 * cp2.cameraPersp;
  cameraZpos := c0 * cp1.cameraZpos + c1 * cp2.cameraZpos;
  cameraDOF := c0 * cp1.cameraDOF + c1 * cp2.cameraDOF;

  for i := 0 to 3 do begin
    curvePoints[i][0].x := c0 * cp1.curvePoints[i][0].x + c1 * cp2.curvePoints[i][0].x;
	  curvePoints[i][0].y := c0 * cp1.curvePoints[i][0].y + c1 * cp2.curvePoints[i][0].y;
	  curveWeights[i][0] := c0 * cp1.curveWeights[i][0] + c1 * cp2.curveWeights[i][0];
    curvePoints[i][1].x := c0 * cp1.curvePoints[i][1].x + c1 * cp2.curvePoints[i][1].x;
	  curvePoints[i][1].y := c0 * cp1.curvePoints[i][1].y + c1 * cp2.curvePoints[i][1].y;
	  curveWeights[i][1] := c0 * cp1.curveWeights[i][1] + c1 * cp2.curveWeights[i][1];
    curvePoints[i][2].x := c0 * cp1.curvePoints[i][2].x + c1 * cp2.curvePoints[i][2].x;
	  curvePoints[i][2].y := c0 * cp1.curvePoints[i][2].y + c1 * cp2.curvePoints[i][2].y;
	  curveWeights[i][2] := c0 * cp1.curveWeights[i][2] + c1 * cp2.curveWeights[i][2];
    curvePoints[i][3].x := c0 * cp1.curvePoints[i][2].x + c1 * cp2.curvePoints[i][2].x;
	  curvePoints[i][3].y := c0 * cp1.curvePoints[i][3].y + c1 * cp2.curvePoints[i][3].y;
	  curveWeights[i][3] := c0 * cp1.curveWeights[i][3] + c1 * cp2.curveWeights[i][3];
  end;

  numTX1 := cp1.NumXForms;
  numTX2 := cp2.NumXForms;
  numTX := max(numTX1, numTX2); // actual xforms including final
  
  f1 := cp1.FinalXformEnabled; // cp1.hasFinalXform;
  f2 := cp2.FinalXformEnabled; // cp2.hasFinalXform;

  if f1 or f2 then
  begin
    if f1 then
    begin
      if numTX1 < numTX then
      begin
        cp1.xform[numTX].Assign(cp1.xform[numTX1]);
        cp1.xform[numTX1].Clear;
        cp1.xform[numTX1].SetVariation(0,0);
      end;
    end
    else
      cp1.xform[numTX].symmetry := 1;

    if f2 then
    begin
      if numTX2 < numTX then
      begin
        cp2.xform[numTX].Assign(cp2.xform[numTX2]);
        cp2.xform[numTX2].Clear;
        cp2.xform[numTX2].SetVariation(0,0);
      end;
    end
    else
      cp2.xform[numTX].symmetry := 1;

    // final XForm
    finalXformEnabled := True;
    xform[numTX].color := c0 * cp1.xform[numTX].color + c1 * cp2.xform[numTX].color;
    xform[numTX].symmetry := c0 * cp1.xform[numTX].symmetry + c1 * cp2.xform[numTX].symmetry;
    xform[numTX].pluginColor := c0 * cp1.xform[numTX].pluginColor + c1 * cp2.xform[numTX].pluginColor;
    xform[numTX].ifs.Assign(cp2.xform[numTX].ifs); // AV: target variation order
    for j := 0 to NrVar-1 do begin
      v1 := cp1.xform[numTX].GetVariation(j);
      v2 := cp2.xform[numTX].GetVariation(j);
      if v1 = v2 then
        xform[numTX].SetVariation(j, v1)
      else
        xform[numTX].SetVariation(j, c0 * v1 + c1 * v2);
    end;
    for j:= 0 to GetNrVariableNames-1 do begin
      cp1.xform[numTX].GetVariable(GetVariableNameAt(j), v1);
      cp2.xform[numTX].GetVariable(GetVariableNameAt(j), v2);
      if v1 <> v2 then v1 := c0 * v1 + c1 * v2;
      xform[numTX].SetVariable(GetVariableNameAt(j), v1);
    end;
    for j := 0 to 2 do begin
      xform[numTX].c[j, 0] := c0 * cp1.xform[numTX].c[j, 0] + c1 * cp2.xform[numTX].c[j, 0];
      xform[numTX].c[j, 1] := c0 * cp1.xform[numTX].c[j, 1] + c1 * cp2.xform[numTX].c[j, 1];
      xform[numTX].p[j, 0] := c0 * cp1.xform[numTX].p[j, 0] + c1 * cp2.xform[numTX].p[j, 0];
      xform[numTX].p[j, 1] := c0 * cp1.xform[numTX].p[j, 1] + c1 * cp2.xform[numTX].p[j, 1];
    end;
  end
  else begin
    finalXformEnabled := False;
	  xform[numTX].symmetry := 1;
  end;

  // regular xforms
  for i := 0 to numTX-1 do begin
    xform[i].density := c0 * cp1.xform[i].density + c1 * cp2.xform[i].density;
    xform[i].color := c0 * cp1.xform[i].color + c1 * cp2.xform[i].color;
    xform[i].symmetry := c0 * cp1.xform[i].symmetry + c1 * cp2.xform[i].symmetry;
	  xform[i].transOpacity := c0 * cp1.xform[i].transOpacity + c1 * cp2.xform[i].transOpacity;
	  xform[i].pluginColor := c0 * cp1.xform[i].pluginColor + c1 * cp2.xform[i].pluginColor;
	  for j := 0 to numTX - 1 do
	    self.xform[i].modWeights[j] := c0 * cp1.xform[i].modWeights[j] + c1 * cp2.xform[i].modWeights[j];

    xform[i].ifs.Assign(cp2.xform[i].ifs); // AV: target variation order
    for j := 0 to NrVar-1 do begin
      v1 := cp1.xform[i].GetVariation(j);
      v2 := cp2.xform[i].GetVariation(j);
      if v1 = v2 then
        xform[i].SetVariation(j, v1)
      else
        xform[i].SetVariation(j, c0 * v1 + c1 * v2);
    end;
    for j:= 0 to GetNrVariableNames-1 do begin
      cp1.xform[i].GetVariable(GetVariableNameAt(j), v1);
      cp2.xform[i].GetVariable(GetVariableNameAt(j), v2);
      if v1 <> v2 then v1 := c0 * v1 + c1 * v2;
      xform[i].SetVariable(GetVariableNameAt(j), v1);
    end;
    // interpol matrices
    for j := 0 to 2 do begin
      xform[i].c[j, 0] := c0 * cp1.xform[i].c[j, 0] + c1 * cp2.xform[i].c[j, 0];
      xform[i].c[j, 1] := c0 * cp1.xform[i].c[j, 1] + c1 * cp2.xform[i].c[j, 1];
	    xform[i].p[j, 0] := c0 * cp1.xform[i].p[j, 0] + c1 * cp2.xform[i].p[j, 0];
      xform[i].p[j, 1] := c0 * cp1.xform[i].p[j, 1] + c1 * cp2.xform[i].p[j, 1];
    end;
  end;

  if numTX1 < numTX then
  begin
    cp1.xform[numTX1].Assign(cp1.xform[numTX]);
    cp1.xform[numTX].Clear;
  end;
  if numTX2 < numTX then
  begin
    cp2.xform[numTX2].Assign(cp2.xform[numTX]);
    cp2.xform[numTX].Clear;
  end;
end;

procedure TControlPoint.InterpolateX(cp1, cp2: TControlPoint; Tm: double);
var
  c0, c1: double;
  i, j: integer;
  v1, v2: double;

  nXforms1, nXmax: integer;
begin
  if (cp2.time - cp1.time) > 1E-6 then begin
    c0 := (cp2.time - tm) / (cp2.time - cp1.time);
    c1 := 1 - c0;
  end else begin
    c0 := 1;
    c1 := 0;
  end;

  // AV: replaced all useless copying into temporary Result variable
  // by direct assignments to self ControlPoint instance
  self.time := Tm;

  {  // AV: since from now this is used only for mutants, we can skip it
  if cp1.cmap_inter = 0 then
    for i := 0 to 255 do begin
      r[0] := cp1.cmap[i][0] / 255;
      r[1] := cp1.cmap[i][1] / 255;
      r[2] := cp1.cmap[i][2] / 255;
      rgb2hsv(r, s);
      r[0] := cp2.cmap[i][0] / 255;
      r[1] := cp2.cmap[i][1] / 255;
      r[2] := cp2.cmap[i][2] / 255;
      rgb2hsv(r, t);
      t[0] := c0 * s[0] + c1 * t[0];
      t[1] := c0 * s[1] + c1 * t[1];
      t[2] := c0 * s[2] + c1 * t[2];
      hsv2rgb(t, r);
      Result.cmap[i][0] := Round(255 * r[0]);
      Result.cmap[i][1] := Round(255 * r[1]);
      Result.cmap[i][2] := Round(255 * r[2]);
    end;
  }
  cmap := cp1.cmap; // AV
  cmapindex := -1;

  Fbrightness := c0 * cp1.Fbrightness + c1 * cp2.Fbrightness;
  contrast := c0 * cp1.contrast + c1 * cp2.contrast;
  gamma := c0 * cp1.gamma + c1 * cp2.gamma;
  vibrancy := c0 * cp1.vibrancy + c1 * cp2.vibrancy;
  gamma_threshold := c0 * cp1.gamma_threshold + c1 * cp2.gamma_threshold;
  width := cp1.width;
  height := cp1.height;
  spatial_oversample := Round(c0 * cp1.spatial_oversample + c1 * cp2.spatial_oversample);
  center[0] := c0 * cp1.center[0] + c1 * cp2.center[0];
  center[1] := c0 * cp1.center[1] + c1 * cp2.center[1];
  FAngle := c0 * cp1.FAngle + c1 * cp2.FAngle;
  pixels_per_unit := c0 * cp1.pixels_per_unit + c1 * cp2.pixels_per_unit;
  {  // AV: since from now this is used only for mutants, we can skip it
  Result.background[0] := c0 * cp1.background[0] + c1 * cp2.background[0];
  Result.background[1] := c0 * cp1.background[1] + c1 * cp2.background[1];
  Result.background[2] := c0 * cp1.background[2] + c1 * cp2.background[2];
  }
  spatial_filter_radius := c0 * cp1.spatial_filter_radius + c1 * cp2.spatial_filter_radius;
  sample_density := c0 * cp1.sample_density + c1 * cp2.sample_density;
  zoom := c0 * cp1.zoom + c1 * cp2.zoom;
  nbatches := Round(c0 * cp1.nbatches + c1 * cp2.nbatches);
  {
  Result.white_level := Round(c0 * cp1.white_level + c1 * cp2.white_level);
  for i := 0 to 3 do begin
    Result.pulse[i div 2][i mod 2] := c0 * cp1.pulse[i div 2][i mod 2] + c1 * cp2.pulse[i div 2][i mod 2];
    Result.wiggle[i div 2][i mod 2] := c0 * cp1.wiggle[i div 2][i mod 2] + c1 * cp2.wiggle[i div 2][i mod 2];
  end;
  }
  cameraPitch := c0 * cp1.cameraPitch;  // AV
  cameraYaw := c0 * cp1.cameraYaw;    // AV
  cameraRoll := c0 * cp1.cameraRoll;  // AV

  // save finalxform from mut(il)ation ;)
  {
  nXforms1 := cp1.NumXForms;
  if cp1.HasFinalXForm then
  begin
    if nXforms1 < NXFORMS then
    begin
      cp1.xform[NXFORMS].Assign(cp1.xform[nXforms1]);
      cp1.xform[nXforms1].Clear;
    end;
  end
  else begin
    cp1.xform[NXFORMS].Clear;
    cp1.xform[NXFORMS].symmetry := 1;
  end;

  nXforms2 := cp2.NumXForms;
  if cp2.HasFinalXForm then
  begin
    if nXforms2 < NXFORMS then
    begin
      cp2.xform[NXFORMS].Assign(cp2.xform[nXforms2]);
      cp2.xform[nXforms2].Clear;
    end;
  end
  else begin
    cp2.xform[NXFORMS].Clear;
    cp2.xform[NXFORMS].symmetry := 1;
  end;
  }
  nXforms1 := cp1.NumXForms;
  nXmax := cp2.NumXForms;
  nXmax := max(nxForms1, nXmax);
  // AV: since we don't change final xforms, we can speed up the calculations
  for i := 0 to nXmax - 1 do
  begin
    xform[i].density := c0 * cp1.xform[i].density + c1 * cp2.xform[i].density;
    xform[i].color := c0 * cp1.xform[i].color + c1 * cp2.xform[i].color;
    xform[i].symmetry := c0 * cp1.xform[i].symmetry + c1 * cp2.xform[i].symmetry;
    xform[i].transOpacity := c0 * cp1.xform[i].transOpacity + c1;  // AV

    for j := 0 to NrVar-1 do begin
      //xform[i].SetVariation(j, c0 * cp1.xform[i].GetVariation(j) + c1 * cp2.xform[i].GetVariation(j));	
      v1 := cp1.xform[i].GetVariation(j);
      v2 := cp2.xform[i].GetVariation(j);
      if v1 = v2 then
        xform[i].SetVariation(j, v1)
      else
        xform[i].SetVariation(j, c0 * v1 + c1 * v2);
    end;
    for j:= 0 to GetNrVariableNames-1 do begin
      cp1.xform[i].GetVariable(GetVariableNameAt(j), v1);
      cp2.xform[i].GetVariable(GetVariableNameAt(j), v2);
      if v1 <> v2 then v1 := c0 * v1 + c1 * v2;  // AV
      xform[i].SetVariable(GetVariableNameAt(j), v1);
    end;
    // interpol matrix
    for j := 0 to 2 do begin
      xform[i].c[j, 0] := c0 * cp1.xform[i].c[j, 0] + c1 * cp2.xform[i].c[j, 0];
      xform[i].c[j, 1] := c0 * cp1.xform[i].c[j, 1] + c1 * cp2.xform[i].c[j, 1];
      // AV
      xform[i].p[j, 0] := c0 * cp1.xform[i].p[j, 0] + c1 * cp2.xform[i].p[j, 0];
      xform[i].p[j, 1] := c0 * cp1.xform[i].p[j, 1] + c1 * cp2.xform[i].p[j, 1];
    end;
  end;

  finalXformEnabled := cp1.finalXformEnabled;
  // finalxform was supposed to be mutate-able too, but somehow it's always
  // getting confused by random-generated mutants :-\

  if cp1.HasFinalXForm then begin
    xform[nXmax].Assign(cp1.xform[nXforms1]);
    if nXforms1 < nXmax then  // reset color speed of cp1's final transform
      xform[nXforms1].symmetry := 0;
  end
  else
    xform[nXmax].Symmetry := 1;
  {
  if NumXForms < NXFORMS then
  begin
    xform[NumXForms].Assign(cp1.xform[NXFORMS]);
    xform[NXFORMS].Clear;
  end;

  // restore finalxforms in source CPs
  if nXforms1 < NXFORMS then
  begin
    cp1.xform[nXforms1].Assign(cp1.xform[NXFORMS]);
    cp1.xform[NXFORMS].Clear;
  end;
  if nXforms2 < NXFORMS then
  begin
    cp2.xform[nXforms2].Assign(cp2.xform[NXFORMS]);
    cp2.xform[NXFORMS].Clear;
  end;
  }
end;

procedure TControlPoint.SaveToFile(Filename: string);
var
  sl: TStringlist;
begin
  sl := TStringlist.Create;

  SaveToStringlist(sl);

  sl.SaveToFile(filename);
  sl.Free;
end;

procedure TControlPoint.SaveToStringlist(sl: TStringlist);
var
  i, j, k: Integer;
  s: string;
  OldDecimalSperator: Char;
  v: double;
  str: string;
  curves: string;
begin
  OldDecimalSperator := FormatSettings.DecimalSeparator;
  FormatSettings.DecimalSeparator := '.';

  sl.add(format('time %f', [time]));
  if cmapindex >= 0 then
    sl.add(format('cmap %d', [cmapindex]));
  sl.add(format('zoom %g', [zoom])); // mt
  sl.add(format('angle %g', [FAngle]));

  sl.add(format('cam_pitch %g', [cameraPitch]));
  sl.add(format('cam_yaw %g', [cameraYaw]));
  sl.add(format('cam_roll %g', [cameraRoll]));  // AV
  sl.add(format('cam_persp %g', [cameraPersp]));
  sl.add(format('cam_zpos %g', [cameraZpos]));
  sl.add(format('cam_dof %g', [cameraDOF]));

  for i := 0 to 3 do
  begin
    curves := curves + FloatToStr(curvePoints[i][0].x) + ' ';
    curves := curves + FloatToStr(curvePoints[i][0].y) + ' ';
    curves := curves + FloatToStr(curveWeights[i][0]) + ' ';

    curves := curves + FloatToStr(curvePoints[i][1].x) + ' ';
    curves := curves + FloatToStr(curvePoints[i][1].y) + ' ';
    curves := curves + FloatToStr(curveWeights[i][1]) + ' ';

    curves := curves + FloatToStr(curvePoints[i][2].x) + ' ';
    curves := curves + FloatToStr(curvePoints[i][2].y) + ' ';
    curves := curves + FloatToStr(curveWeights[i][2]) + ' ';

    curves := curves + FloatToStr(curvePoints[i][3].x) + ' ';
    curves := curves + FloatToStr(curvePoints[i][3].y) + ' ';
    curves := curves + FloatToStr(curveWeights[i][3]) + ' ';
  end;

  curves := trim(curves);
  sl.Add(Format('curves %s', [curves]));

  sl.add(format('image_size %d %d center %g %g pixels_per_unit %f',
    [Width, Height, center[0], center[1], pixels_per_unit]));
  sl.add(format('spatial_oversample %d spatial_filter_radius %f',
    [spatial_oversample, spatial_filter_radius]));
  sl.add(format('sample_density %g', [sample_density]));
//  sl.add(format('nbatches %d white_level %d background %f %f %f', - changed to integers - mt
  sl.add(format('nbatches %d white_level %d background %d %d %d',
    [nbatches, white_level, background[0], background[1], background[2]]));
  sl.add(format('brightness %f gamma %f contrast %f vibrancy %f gamma_threshold %f hue_rotation %f', {cmap_inter %d',}
    [Fbrightness * BRIGHT_ADJUST, gamma, contrast, vibrancy, gamma_threshold, hue_rotation {, cmap_inter}]));  // AV: added contrast
  sl.add(format('finalxformenabled %d', [ifthen(finalxformenabled, 1, 0)]));
  sl.add(format('soloxform %d', [soloXform]));

  (*str := '';
  for i := 0 to used_plugins.Count-1 do begin
    str := str + used_plugins[i];
    if (i = used_plugins.Count-1) then break;
    str := str + ' ';
  end;
  sl.Add(format('plugins %s', [str]));   *)

  for i := 0 to Min(NumXForms+1, NXFORMS) do
    with xform[i] do begin
      //if density = 0 then continue; - FinalXform has weight=0

      sl.add(format('xform %d density %g color %g symmetry %g', [i, density, color, symmetry]));
	  
	 /////// AV: try to write variation order
	  s := 'var_order'; 
	  for j := 0 to NRVAR - 1 do begin
        s := format('%s %d', [s, GetVariationIndex(ifs[j])]);
      end;
      sl.add(s);
	//////////////////////////
	
      s := 'vars';
      for j := 0 to NRVAR - 1 do begin
        s := format('%s %g', [s, GetVariation(j)]);
      end;
      sl.add(s);
      s := 'variables';
      for j:= 0 to GetNrVariableNames-1 do begin
{$ifndef VAR_STR}
        GetVariable(GetVariableNameAt(j), v);
        s := format('%s %g', [s, v]);
{$else}
        s := s + ' ' + GetVariableStr(GetVariableNameAt(j));
{$endif}
      end;
      sl.add(s);
      sl.Add(format('coefs %.6f %.6f %.6f %.6f %.6f %.6f',
        [c[0][0], c[0][1], c[1][0], c[1][1], c[2][0], c[2][1]]));
      sl.Add(format('post %.6f %.6f %.6f %.6f %.6f %.6f',
        [p[0][0], p[0][1], p[1][0], p[1][1], p[2][0], p[2][1]]));
      if postXswap then
        sl.Add('postxswap 1')
      else
        sl.Add('postxswap 0');
      if autoZscale then
        sl.Add('autozscale 1')
      else
        sl.Add('autozscale 0');
      s := 'chaos';
      for j := 0 to NumXForms+1 do begin
        s := s + format(' %g', [modWeights[j]]);
      end;
      sl.Add(s);

      sl.Add(format('opacity %g', [transOpacity]));
      sl.Add(format('var_color %g', [pluginColor]));

    end;
  FormatSettings.DecimalSeparator := OldDecimalSperator;
end;

procedure WriteDoubles(const handle: File; data: array of double);
var
  block: TBlock;
  i: integer;
begin
  for i := 0 to Length(data)-1 do begin
    DoubleToBlock(block, 0, 0); //pad to blocksize
    DoubleToBlock(block, 8, data[i]);
    BlockWrite(handle, block, 1);
  end;
end;
procedure WriteString(const handle: File; data: string);
var
  k, l, size, chunks: Integer;
  raw : THibRawString;
  block: TBlock;
begin
  size := Length(data);
  SetLength(raw, size);
  CopyMemory(@raw[0], @data[1], size);
  chunks := size div HIB_BLOCKSIZE;

  if size mod HIB_BLOCKSIZE > 0 then begin
    size := (1 + size div HIB_BLOCKSIZE) * HIB_BLOCKSIZE;
    chunks := chunks + 1;
  end;

  for k := 0 to chunks - 1 do begin
    for l := 0 to HIB_MAXOFFSET do
      if (k * HIB_BLOCKSIZE + l) < size then
        block[l] := raw[k * HIB_BLOCKSIZE + l]
      else block[l] := 0;
    BlockWrite(handle, block, 1);
  end;
end;

function CalcBinaryFlameSize(cp: TControlPoint): integer;
var
  (*str: string;
  i, nvariations, nvariables, nchaos: Integer;*)
  handle: File;
begin
  // I'm a bit ashamed but this hack has do to it for now...
  AssignFile(handle, AppData + 'CalcBinaryFlameSizeTemp.bin');
  ReWrite(handle, HIB_BLOCKSIZE);
  cp.SaveToBinary(handle);
  Result := FileSize(handle) * HIB_BLOCKSIZE;
  CloseFile(handle);
  DeleteFile(AppData + 'CalcBinaryFlameSizeTemp.bin');

  (*// CP data
  Result := 224;

  // Var list
  str := '';
  for i := 0 to NRVAR-1 do
    str := str + VarNames(i) + #13#10;
  for i:= 0 to GetNrVariableNames-1 do
    str := str + GetVariableNameAt(i)+ #13#10;
  str := trim(str);
  if Length(str) mod HIB_BLOCKSIZE > 0 then
    Result := Result + (1 + Length(str) div HIB_BLOCKSIZE) * HIB_BLOCKSIZE
  else
    Result := Result + Length(str);

  // XForm data
  nchaos := Min(cp.NumXForms+1, NXFORMS);
  nvariations := NRVAR;
  nvariables := GetNrVariableNames;
  if nvariations mod 2 > 0 then
    nvariations := nvariations + 1;
  if nvariables mod 2 > 0 then
    nvariables := nvariables + 1;
  if nchaos mod 2 > 0 then
    nchaos := nchaos + 1;
  Result := Result + (144 + (nvariations + nvariables + nchaos) * 8) *
    Min(cp.NumXForms+1, NXFORMS);  *)
end;

procedure TControlPoint.SaveToBinary(const handle: File);  
var
  i, j, nvariations, nvariables, nchaos: Integer;
  v: double;
  str: string;
  dbl: array of double;
  block: TBlock;
begin

  DoubleToBlock(block, 0, time);
  DoubleToBlock(block, 8, zoom);
  BlockWrite(handle, block, 1);

  DoubleToBlock(block, 0, FAngle);
  DoubleToBlock(block, 8, pixels_per_unit);
  BlockWrite(handle, block, 1);

  DoubleToBlock(block, 0, center[0]);
  DoubleToBlock(block, 8, center[1]);
  BlockWrite(handle, block, 1);

  DoubleToBlock(block, 0, cameraPitch);
  DoubleToBlock(block, 8, cameraYaw);
  BlockWrite(handle, block, 1);

  DoubleToBlock(block, 0, cameraRoll); // AV
  DoubleToBlock(block, 8, cameraZpos);
  BlockWrite(handle, block, 1);
 
  DoubleToBlock(block, 0, cameraPersp);
  DoubleToBlock(block, 8, cameraDOF);
  BlockWrite(handle, block, 1);

  DoubleToBlock(block, 0, contrast); // AV
  DoubleToBlock(block, 8, spatial_filter_radius);
  BlockWrite(handle, block, 1);

  DoubleToBlock(block, 0, sample_density);
  DoubleToBlock(block, 8, gamma_threshold);
  BlockWrite(handle, block, 1);

  //Int32ToBlock(block, 0, cmapindex);// TODO: replace
  Int32ToBlock(block, 4, spatial_oversample);
  Int32ToBlock(block, 8, Width);
  Int32ToBlock(block, 12, Height);
  BlockWrite(handle, block, 1);

  Int32ToBlock(block, 0, nbatches);
  //Int32ToBlock(block, 4, cmap_inter); // TODO: replace
  Int32ToBlock(block, 8, ifthen(finalxformenabled, 1, 0));
  Int32ToBlock(block, 12, soloXform);
  BlockWrite(handle, block, 1);

  //Int32ToBlock(block, 0, white_level); // TODO: replace
  Int32ToBlock(block, 4, background[0]);
  Int32ToBlock(block, 8, background[1]);
  Int32ToBlock(block, 12, background[2]);
  BlockWrite(handle, block, 1);

  DoubleToBlock(block, 0, Fbrightness * BRIGHT_ADJUST);
  DoubleToBlock(block, 8, gamma);
  BlockWrite(handle, block, 1);

  DoubleToBlock(block, 0, vibrancy);
  DoubleToBlock(block, 8, hue_rotation);
  BlockWrite(handle, block, 1);

  nchaos := Min(NumXForms+1, NXFORMS);
  nvariations := NRVAR;
  nvariables := GetNrVariableNames;

  Int32ToBlock(block, 0, nvariations);
  Int32ToBlock(block, 4, nvariables);
  Int32ToBlock(block, 8, nchaos);

// AV TODO: var_order   
  str := '';
  for i := 0 to nvariations-1 do
    str := str + VarNames(i) + #0;
  for i:= 0 to nvariables-1 do
    str := str + GetVariableNameAt(i)+ #0;
  str := trim(str);
  if Length(str) mod HIB_BLOCKSIZE > 0 then
    Int32ToBlock(block, 12, (1 + Length(str) div HIB_BLOCKSIZE) * HIB_BLOCKSIZE)
  else
    Int32ToBlock(block, 12, Length(str));
  BlockWrite(handle, block, 1);

  WriteString(handle, str);

  for i := 0 to nchaos - 1 do
    with xform[i] do begin
      DoubleToBlock(block, 0, density);
      DoubleToBlock(block, 8, color);
      BlockWrite(handle, block, 1);

      DoubleToBlock(block, 0, symmetry);
      DoubleToBlock(block, 8, transOpacity);
      BlockWrite(handle, block, 1);

      DoubleToBlock(block, 0, pluginColor);
      Int32ToBlock(block, 8, IfThen(postXswap, 1, 0));
      Int32ToBlock(block, 12, IfThen(autozscale, 1, 0));
      BlockWrite(handle, block, 1);

      DoubleToBlock(block, 0, c[0][0]);
      DoubleToBlock(block, 8, c[0][1]);
      BlockWrite(handle, block, 1);

      DoubleToBlock(block, 0, c[1][0]);
      DoubleToBlock(block, 8, c[1][1]);
      BlockWrite(handle, block, 1);

      DoubleToBlock(block, 0, c[2][0]);
      DoubleToBlock(block, 8, c[2][1]);
      BlockWrite(handle, block, 1);

      DoubleToBlock(block, 0, p[0][0]);
      DoubleToBlock(block, 8, p[0][1]);
      BlockWrite(handle, block, 1);

      DoubleToBlock(block, 0, p[1][0]);
      DoubleToBlock(block, 8, p[1][1]);
      BlockWrite(handle, block, 1);

      DoubleToBlock(block, 0, p[2][0]);
      DoubleToBlock(block, 8, p[2][1]);
      BlockWrite(handle, block, 1);

      SetLength(dbl, nvariations);
      for j := 0 to nvariations-1 do
        dbl[j] := GetVariation(j);
      WriteDoubles(handle, dbl);

      SetLength(dbl, nvariables);
      for j := 0 to nvariables-1 do
        GetVariable(GetVariableNameAt(j), dbl[j]);
      WriteDoubles(handle, dbl);

      SetLength(dbl, nchaos);
      for j := 0 to nchaos-1 do
        dbl[j] := xform[i].modWeights[j];
      WriteDoubles(handle, dbl);
    end;
end;

procedure TControlPoint.DirectCopy_AV(const cp1: TControlPoint); // AV
var
  i, j: smallint;
begin
	self.name := cp1.name;
	self.time := cp1.time;
  
	if cp1.cmapindex >= 0 then cmapindex := cp1.cmapindex;
	hue_rotation := cp1.hue_rotation;  // AV
	cmap := cp1.cmap;

	zoom := cp1.zoom;
	Fangle := cp1.FAngle;
	cameraPitch := cp1.cameraPitch;
	cameraYaw := cp1.cameraYaw;
	cameraRoll := cp1.cameraRoll;
	cameraPersp := cp1.cameraPersp;
	cameraZpos := cp1.cameraZpos;
	cameraDOF := cp1.cameraDOF;

	Width := cp1.Width;
	Height := cp1.Height;

	center[0]:= cp1.center[0];
	center[1]:= cp1.center[1];
	pixels_per_unit := cp1.pixels_per_unit;

	background[0] := cp1.background[0];
	background[1] := cp1.background[1];
	background[2] := cp1.background[2];

	for i := 0 to 3 do
		for j := 0 to 3 do begin
			curveWeights[i,j] := cp1.curveWeights[i,j];
			curvePoints[i,j].x := cp1.curvePoints[i,j].x;
			curvePoints[i,j].y := cp1.curvePoints[i,j].y;
		end;

	spatial_oversample := cp1.spatial_oversample;
	spatial_filter_radius := cp1.spatial_filter_radius;
	sample_density := cp1.sample_density;

	nbatches := cp1.nbatches;
	//white_level := cp1.white_level;  // <-- AV: it's const for now
	//jitters := cp1.jitters;  //<-- AV: currently unused
	//Transparency := cp1.Transparency;  // <-- do this?
	
	brightness := cp1.brightness;
	gamma := cp1.gamma;
	gamma_threshold := cp1.gamma_threshold;
	contrast := cp1.contrast; // AV
	vibrancy := cp1.vibrancy;

	finalXformEnabled := cp1.finalXformEnabled;
	soloXform := cp1.soloXform;

	estimator := cp1.estimator;
	estimator_min := cp1.estimator_min;
	estimator_curve := cp1.estimator_curve;
	enable_de := cp1.enable_de;

	comment := cp1.comment; // AV

	// used_plugins := cp1.used_plugins;  // <<--- fixed memory leak // AV
	used_plugins.Clear; // AV: that's how we must copy strings:
	used_plugins.AddStrings(cp1.used_plugins);
	
	for i := 0 to NXFORMS do // was: NXFORMS-1
		xform[i].assign(cp1.xform[i]);
end;

function TControlPoint.Clone: TControlPoint;
(*
var
  i, j: integer;
  sl: TStringList;
begin
  sl := TStringList.Create;
  SaveToStringlist(sl);
  Result := TControlPoint.Create;
  Result.ParseStringlist(sl);
  Result.Fangle := FAngle;
  Result.cmap := cmap;
  Result.name := name;
  //Result.nick := nick;
  //Result.url := url;
  Result.hue_rotation := hue_rotation;  // AV
  Result.Transparency := Transparency;
  Result.gamma_threshold := gamma_threshold;
  Result.estimator := estimator;
  Result.estimator_min := estimator_min;
  Result.estimator_curve := estimator_curve;
  Result.enable_de := enable_de;
  // Result.xdata := xdata;
  Result.comment := comment; // AV

  Result.Background[0] := background[0];
  Result.Background[1] := background[1];
  Result.Background[2] := background[2];

  for i := 0 to 3 do
    for j := 0 to 3 do begin
      Result.CurveWeights[i,j] := curveWeights[i,j];
      Result.curvePoints[i,j].x := curvePoints[i,j].x;
      Result.curvePoints[i,j].y := curvePoints[i,j].y;
    end;

  Result.used_plugins.Clear;
  Result.used_plugins.AddStrings(used_plugins);  // AV
 { for i := 0 to used_plugins.Count-1 do
    Result.used_plugins.Add(used_plugins[i]); }

  for i := 0 to NXFORMS - 1 do
    Result.xform[i].assign(xform[i]);

  sl.Free;
  *)
begin
  Result := TControlPoint.Create;
	Result.DirectCopy_AV(self);  // AV
  Result.Transparency := Transparency;
end;

procedure TControlPoint.Copy(cp1: TControlPoint; KeepSizes: boolean = false);
var
  //i, j: integer;
  //sl: TStringList;
  w, h: integer;
begin
  w := Width;
  h := Height;

  Clear;
  DirectCopy_AV(cp1);  // AV: made it faster than ever!

  if KeepSizes then  AdjustScale(w, h);
  (*
  sl := TStringList.Create;

  // --Z-- this is quite a weird and unoptimal way to copy things:
  cp1.SaveToStringlist(sl);  // <-- AV: so we must optimize it ;-)
  ParseStringlist(sl);

  Fangle := cp1.FAngle;
  center[0]:= cp1.center[0];
  center[1]:= cp1.center[1];
  pixels_per_unit := cp1.pixels_per_unit;
  cmap := cp1.cmap;
  name := cp1.name;
  //nick := cp1.nick;
  //url := cp1.url;
  hue_rotation := cp1.hue_rotation;  // AV
  gamma_threshold := cp1.gamma_threshold;
  estimator := cp1.estimator;
  estimator_min := cp1.estimator_min;
  estimator_curve := cp1.estimator_curve;
  enable_de := cp1.enable_de;
 // used_plugins := cp1.used_plugins;  // <<--- fixed memory leak // AV
  // xdata := cp1.xdata;
  comment := cp1.comment; // AV

  background[0] := cp1.background[0];
  background[1] := cp1.background[1];
  background[2] := cp1.background[2];

  for i := 0 to 3 do
    for j := 0 to 3 do begin
      CurveWeights[i,j] := cp1.curveWeights[i,j];
      curvePoints[i,j].x := cp1.curvePoints[i,j].x;
      curvePoints[i,j].y := cp1.curvePoints[i,j].y;
    end;

  if KeepSizes then
    AdjustScale(w, h);

  used_plugins.Clear;
  // AV: that's how we must copy strings:
  used_plugins.AddStrings(cp1.used_plugins);
 { for i := 0 to cp1.used_plugins.Count-1 do
    used_plugins.Add(cp1.used_plugins[i]); }

  // AV: to speed up the interpolation
  for i := 0 to nt {NXFORMS} do // was: NXFORMS-1
    xform[i].assign(cp1.xform[i]);
  finalXformEnabled := cp1.finalXformEnabled;

  sl.Free;
  *)
end;

procedure TControlPoint.ParseStringList(sl: TStringlist);
var
  s: string;
  i: integer;
begin
  finalXformEnabled := false;
  for i := 0 to sl.Count - 1 do begin
    s := s + sl[i] + ' ';
  end;
  ParseString(s);
end;

procedure TControlPoint.Clear;
var
  i, j: Integer;
begin
  symmetry := 0;
  cmapindex := -1;
  hue_rotation := 1; // AV
  zoom := 0;
  comment := '';  // AV
  for i := 0 to NXFORMS do xform[i].Clear;
  FinalXformEnabled := false;
  soloxform := -1;

  for i := 0 to 3 do
  begin
    curvePoints[i][0].x := 0.00; curvePoints[i][0].y := 0.00; curveWeights[i][0] := 1;
    curvePoints[i][1].x := 0.00; curvePoints[i][1].y := 0.00; curveWeights[i][1] := 1;
    curvePoints[i][2].x := 1.00; curvePoints[i][2].y := 1.00; curveWeights[i][2] := 1;
    curvePoints[i][3].x := 1.00; curvePoints[i][3].y := 1.00; curveWeights[i][3] := 1;
  end;

  if (used_plugins <> nil) then
      used_plugins.Clear;
   // else used_plugins := TStringlist.Create;  // AV: now it's not needed
end;

function TControlPoint.HasFinalXForm: boolean;
var
  i: integer;
begin
  with xform[NumXForms] do
  begin
    Result := (c[0,0]<>1) or (c[0,1]<>0) or (c[1,0]<>0) or (c[1,1]<>1) or (c[2,0]<>0) or (c[2,1]<>0) or
              (p[0,0]<>1) or (p[0,1]<>0) or (p[1,0]<>0) or (p[1,1]<>1) or (p[2,0]<>0) or (p[2,1]<>0) or
              (symmetry <> 1) or (GetVariation(0) <> 1);
    if Result = false then
      for i := 1 to NRVAR-1 do Result := Result or (GetVariation(i) <> 0);
  end;
end;

function add_symmetry_to_control_point(var cp: TControlPoint; sym: integer): integer;
var
  i, j, k: integer;
  a: double;
begin
  result := 0;
  if (0 = sym) then
    sym := random(13) - 6;

  if (1 = sym) or (0 = sym) then
    Exit(0);   // AV

  i := cp.NumXForms;  // AV

  if (i + abs(sym)>= NXFORMS) then  // AV: take into account additional xforms
    Exit(0);

  cp.symmetry := sym;

  if (sym < 0) then
  begin
    cp.xform[i].density := 1.0;
    cp.xform[i].symmetry := 1;
    cp.xform[i].SetVariation(0, 1.0);
    for j := 1 to NRVAR - 1 do
      cp.xform[i].SetVariation(j, 0.0);
    cp.xform[i].color := 1.0;  // AV: why if color_speed = 1?!
    cp.xform[i].c[0][0] := -1.0;
    cp.xform[i].c[0][1] := 0.0;
    cp.xform[i].c[1][0] := 0.0;
    cp.xform[i].c[1][1] := 1.0;
    cp.xform[i].c[2][0] := 0.0;
    cp.xform[i].c[2][1] := 0.0;

    inc(i);
    inc(result);
    sym := -sym;
  end;

  a := 2 * PI / sym;

  k := 1;
  while (k < sym) and (i < SymmetryNVars) do
  begin
    cp.xform[i].density := 1.0;
    cp.xform[i].SetVariation(0, 1);
    cp.xform[i].symmetry := 1;
    for j := 1 to NRVAR - 1 do
      cp.xform[i].SetVariation(j, 0);
    if sym < 3 then
      cp.xform[i].color := 0  // AV: why if color_speed = 1?!
    else
      cp.xform[i].color := (k - 1) / (sym - 2);

    while cp.xform[i].color > 1 do  // AV: why if color_speed = 1?!
        cp.xform[i].color := cp.xform[i].color - 1;

    cp.xform[i].c[0][0] := cos(k * a);
    cp.xform[i].c[0][1] := sin(k * a);
    cp.xform[i].c[1][0] := -cp.xform[i].c[0][1];
    cp.xform[i].c[1][1] := cp.xform[i].c[0][0];
    cp.xform[i].c[2][0] := 0.0;
    cp.xform[i].c[2][1] := 0.0;

    inc(i);
    inc(result);
    inc(k);
  end;
end;

procedure tile_control_point(var cp: TControlPoint; sym: integer); // AV
var
    k, i, j: smallint;
    t, cost, sint: double;
begin
  k := cp.NumXForms;
  if (k > NXFORMS - 5) then exit;

  case sym of
    1:
    for i := k to (k + 3) do
    begin
      cp.xform[i].density := 1.0;
      cp.xform[i].symmetry := 1;
      cp.xform[i].SetVariation(0, 1.0);
      for j := 1 to NRVAR - 1 do
        cp.xform[i].SetVariation(j, 0.0);
      cp.xform[i].c[0][0] := 1.0;
      cp.xform[i].c[0][1] := 0.0;
      cp.xform[i].c[1][0] := 0.0;
      cp.xform[i].c[1][1] := 1.0;
      t := (i - k) * pi * 0.5;
      SinCos(t, sint, cost);
      cp.xform[i].c[2][0] := RhombTR * cost;
      cp.xform[i].c[2][1] := RhombTR * sint;
    end;
    2:
    for i := k to (k + 3) do
    begin
      cp.xform[i].density := 1.0;
      cp.xform[i].symmetry := 1;
      cp.xform[i].SetVariation(0, 1.0);
      for j := 1 to NRVAR - 1 do
        cp.xform[i].SetVariation(j, 0.0);
      cp.xform[i].c[0][0] := 1.0;
      cp.xform[i].c[0][1] := 0.0;
      cp.xform[i].c[1][0] := 0.0;
      cp.xform[i].c[1][1] := 1.0;
      t := pi * ((i - k) * 0.5 + 0.25);
      SinCos(t, sint, cost);
      cp.xform[i].c[2][0] := SquareTR * cost;
      cp.xform[i].c[2][1] := SquareTR * sint;
    end;
    3:
    for i := k to (k + 5) do
    begin
      cp.xform[i].density := 1.0;
      cp.xform[i].symmetry := 1;
      cp.xform[i].SetVariation(0, 1.0);
      for j := 1 to NRVAR - 1 do
        cp.xform[i].SetVariation(j, 0.0);
      cp.xform[i].c[0][0] := 1.0;
      cp.xform[i].c[0][1] := 0.0;
      cp.xform[i].c[1][0] := 0.0;
      cp.xform[i].c[1][1] := 1.0;
      t := (i - k) * pi / 3;
      SinCos(t, sint, cost);
      cp.xform[i].c[2][0] := HexTR * cost;
      cp.xform[i].c[2][1] := HexTR * sint;
    end;
  end;
end;

 //***************************************************************

function det(a, b, c, d: double): double;
begin
  Result := (a * d - b * c);
end;

function solve3(x1, x2, x1h, y1, y2, y1h, z1, z2, z1h: double;
  var a, b, e: double): double;
var
  det1: double;
begin
  det1 := x1 * det(y2, 1.0, z2, 1.0) - x2 * det(y1, 1.0, z1, 1.0)
    + 1 * det(y1, y2, z1, z2);
  if (det1 = 0.0) then
  begin
    Result := det1;
    EXIT;
  end
  else
  begin
    a := (x1h * det(y2, 1.0, z2, 1.0) - x2 * det(y1h, 1.0, z1h, 1.0)
      + 1 * det(y1h, y2, z1h, z2)) / det1;
    b := (x1 * det(y1h, 1.0, z1h, 1.0) - x1h * det(y1, 1.0, z1, 1.0)
      + 1 * det(y1, y1h, z1, z1h)) / det1;
    e := (x1 * det(y2, y1h, z2, z1h) - x2 * det(y1, y1h, z1, z1h)
      + x1h * det(y1, y2, z1, z2)) / det1;
    a := Round6(a);
    b := Round6(b);
    e := Round6(e);
    Result := det1;
  end;
end;

function dist(x1, y1, x2, y2: double): double;
//var
//  d2: double;
begin
(*
  { From FDesign source
  { float pt_pt_distance(float x1, float y1, float x2, float y2) }
  d2 := (x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2);
  if (d2 = 0.0) then
  begin
    Result := 0.0;
    exit;
  end
  else
    Result := sqrt(d2);
*)

  // --Z-- This is just amazing... :-\
  // Someone needed an 'FDesign source' -  to compute distance between two points??!?

  Result := Hypot(x2-x1, y2-y1);
end;

function line_dist(x, y, x1, y1, x2, y2: double): double;
var
  a, b, e, c: double;
begin
  if ((x = x1) and (y = y1)) then
    a := 0.0
  else
    a := sqrt((x - x1) * (x - x1) + (y - y1) * (y - y1));
  if ((x = x2) and (y = y2)) then
    b := 0.0
  else
    b := sqrt((x - x2) * (x - x2) + (y - y2) * (y - y2));
  if ((x1 = x2) and (y1 = y2)) then
    e := 0.0
  else
    e := sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2));
  if ((a * a + e * e) < (b * b)) then
    Result := a
  else if ((b * b + e * e) < (a * a)) then
    Result := b
  else if (e <> 0.0) then
  begin
    c := (b * b - a * a - e * e) / (-2 * e);
    if ((a * a - c * c) < 0.0) then
      Result := 0.0
    else
      Result := sqrt(a * a - c * c);
  end
  else
    Result := a;
end;

//**************************************************************************//

procedure TControlPoint.ZoomtoRect(R: TSRect);
var
  scale, ppu: double;
  dx,dy: double;
begin
  scale := power(2, zoom);
  ppu := pixels_per_unit * scale;

  dx := ((r.Left + r.Right)/2 - Width/2) / ppu;
  dy := ((r.Top + r.Bottom)/2 - Height/2) / ppu;

  center[0] := center[0] + cos(FAngle) * dx - sin(FAngle) * dy;
  center[1] := center[1] + sin(FAngle) * dx + cos(FAngle) * dy;

  if PreserveQuality then
    zoom := Log2(scale * ( Width/(abs(r.Right - r.Left) + 1)))
  else
    pixels_per_unit := pixels_per_unit * Width / abs(r.Right - r.Left);
end;

procedure TControlPoint.ZoomOuttoRect(R: TSRect);
var
  ppu: double;
  dx, dy: double;
begin

  if PreserveQuality then
    zoom := Log2(power(2, zoom) / ( Width/(abs(r.Right - r.Left) + 1)))
  else
    pixels_per_unit := pixels_per_unit / Width * abs(r.Right - r.Left);
  ppu := pixels_per_unit * power(2, zoom);

  dx := ((r.Left + r.Right)/2 - Width/2) / ppu;
  dy := ((r.Top + r.Bottom)/2 - Height/2) / ppu;

  center[0] := center[0] - cos(FAngle) * dx + sin(FAngle) * dy;
  center[1] := center[1] - sin(FAngle) * dx - cos(FAngle) * dy;
end;

procedure TControlPoint.ZoomIn(Factor: double);
var
  scale: double;
begin
  scale := power(2, zoom);

  Scale := Scale / Factor;
  Zoom := Log2(Scale);
end;

procedure TControlPoint.MoveRect(R: TSRect);
var
  scale: double;
  ppux, ppuy: double;
  dx,dy: double;
begin
  scale := power(2, zoom);
  ppux := pixels_per_unit * scale;
  ppuy := pixels_per_unit * scale;

  dx := (r.Left - r.Right)/ppux;
  dy := (r.Top - r.Bottom)/ppuy;

  center[0] := center[0] + cos(FAngle) * dx - sin(FAngle) * dy;
  center[1] := center[1] + sin(FAngle) * dx + cos(FAngle) * dy ;
end;

procedure TControlPoint.Rotate(Angle: double);
begin
  FAngle := FAngle + Angle;
end;

{//////////////////////////////////////////////////////////////////////////////}

function TControlPoint.getppux: double;
begin
  result := pixels_per_unit * power(2, zoom)
end;

function TControlPoint.getppuy: double;
begin
  result := pixels_per_unit * power(2, zoom)
end;

{/////////////////////////////////////////////////////////////////////////////}
function TControlPoint.GetBrightness: double;
begin
  Result := Fbrightness;
end;

procedure TControlPoint.SetBrightness(br: double);
begin
  if br > 0 then begin
    if Fbrightness <> 0 then gamma_threshold := (gamma_threshold / Fbrightness) * br;
    Fbrightness := br;
  end;
end;

{//////////////////////////////////////////////////////////////////////////////}
function TControlPoint.GetRelativeGammaThreshold: double;
begin
  if Fbrightness <> 0 then
    Result := gamma_threshold / Fbrightness
  else
    Result := gamma_threshold;  
end;

procedure TControlPoint.SetRelativeGammaThreshold(gtr: double);
begin
  gamma_threshold := gtr * Fbrightness;
end;

{//////////////////////////////////////////////////////////////////////////////}
(*
var
  vdfilled: boolean = False;

procedure FillVarDisturb;
const
  startvar_distrib: array[0..26] of integer = (-1, -1, -1, -1, -1, -1, -1, 0, 0, 0, 0, 1, 1, 1, 2, 2, 2, 3, 3, 4, 4, 5, 5, 6, 6, 7, 7);
  startmixed_var_distrib: array[0..16] of integer = (0, 0, 0, 1, 1, 1, 2, 2, 2, 3, 3, 4, 4, 5, 6, 6, 7);
var
  i: integer;
begin
  if vdfilled then
    Exit;

  setlength(var_distrib, NRVAR + 19);
  setlength(mixed_var_distrib, NRVAR + 9);

  for i := 0 to High(startvar_distrib) do
    var_distrib[i] := startvar_distrib[i];

  for i := High(startvar_distrib) + 1 to high(var_distrib) do
    var_distrib[i] := 8 + i - High(startvar_distrib) - 1;

  for i := 0 to High(startmixed_var_distrib) do
    mixed_var_distrib[i] := startmixed_var_distrib[i];

  for i := High(startmixed_var_distrib) + 1 to high(mixed_var_distrib) do
    mixed_var_distrib[i] := 8 + i - High(startmixed_var_distrib) - 1;

  vdfilled := true;
end;
 *)
{////////////////////////////////////////////////////////////////////////////}
//
// --Z-- cp-specific functions moved here from MainForm
//

function TControlPoint.NumXForms: integer;
var
  i: smallint;
begin
//...
  Result := NXFORMS;
  for i := 0 to NXFORMS - 1 do
  begin
    if xform[i].density = 0 then
    begin
      Result := i;
      Break;
    end;
  end;
end;

function TControlPoint.TrianglesFromCP(var Triangles: TTriangles): integer;
{ Sets up the triangles from the IFS code }
var
  i, j: smallint;
  //temp_x, temp_y, xset, yset: double;
  //left, top, bottom, right: double;
begin
  Result := NumXForms;
{
  top := 0; bottom := 0; right := 0; left := 0;

  if ReferenceMode > 0 then
  begin
    for i := 0 to Result-1 do
    begin
      xset := 1.0;
      yset := 1.0;
      for j := 0 to 5 do
      with xform[i] do begin
        temp_x := xset * c[0][0] + yset * c[1][0] + c[2][0];
        temp_y := xset * c[0][1] + yset * c[1][1] + c[2][1];
        xset := temp_x;
        yset := temp_y;
      end;
      if (i = 0) then
      begin
        left := xset;
        right := xset;
        top := yset;
        bottom := yset;
      end
      else
      begin
        if (xset < left) then left := xset;
        if (xset > right) then right := xset;
        if (yset > top) then top := yset;
        if (yset < bottom) then bottom := yset;
      end;
    end;

    if ReferenceMode = 1 then
    begin
      Triangles[-1].x[0] := right-left;
      Triangles[-1].y[0] := 0;
      Triangles[-1].x[1] := 0;
      Triangles[-1].y[1] := 0;
      Triangles[-1].x[2] := 0;
      Triangles[-1].y[2] := -(top-bottom);
    end
    else begin
      Triangles[-1].x[0] := right;
      Triangles[-1].y[0] := -bottom;
      Triangles[-1].x[1] := left;
      Triangles[-1].y[1] := -bottom;
      Triangles[-1].x[2] := left;
      Triangles[-1].y[2] := -top;
    end;
  end
  else
}
  begin
    Triangles[-1].x[0] := 1; Triangles[-1].y[0] := 0; // "x"
    Triangles[-1].x[1] := 0; Triangles[-1].y[1] := 0; // "0"
    Triangles[-1].x[2] := 0; Triangles[-1].y[2] := -1; // "y"
  end;

  for j := 0 to Result do
  begin
    for i := 0 to 2 do
    with xform[j] do begin
      if postXswap then begin
        Triangles[j].x[i] := Triangles[-1].x[i] * p[0][0] + Triangles[-1].y[i] * p[1][0] + p[2][0];
        Triangles[j].y[i] := Triangles[-1].x[i] * p[0][1] + Triangles[-1].y[i] * p[1][1] + p[2][1];
      end
      else begin
        Triangles[j].x[i] := Triangles[-1].x[i] * c[0][0] + Triangles[-1].y[i] * c[1][0] + c[2][0];
        Triangles[j].y[i] := Triangles[-1].x[i] * c[0][1] + Triangles[-1].y[i] * c[1][1] + c[2][1];
      end;
    end;
  end;
  EnableFinalXform := FinalXformEnabled;

  // I don't like this... :-/
  for j := -1 to Result do // was: Result-1
    for i := 0 to 2 do
      Triangles[j].y[i] := -Triangles[j].y[i];
end;

//************ Weight distribution utils *************************//

(*  // AV: we have the save methods in RndFlame module
procedure TControlPoint.NormalizeWeights;
var
  i: integer;
  td: double;
begin
  td := 0.0;
  for i := 0 to NumXForms - 1 do
    td := td + xform[i].Density;
  if (td < 0.001) then
    EqualizeWeights
  else
    for i := 0 to NumXForms - 1 do
      xform[i].Density := xform[i].Density / td;
end;

procedure TControlPoint.ComputeWeights(Triangles: TTriangles; t: integer);
// Caclulate transform weight from triangle areas
var
  i: integer;
  total_area: double;
begin
  total_area := 0;
  for i := 0 to t - 1 do
  begin
    xform[i].Density := triangle_area(Triangles[i]);
    total_area := total_area + xform[i].Density;
  end;
  for i := 0 to t - 1 do
  begin
    xform[i].Density := xform[i].Density / total_area;
  end;
  //? cp1.NormalizeWeights;
end;
*)

 procedure TControlPoint.NormalizeProbabilities;
 // AV: useful for matematicians, useless for other people :)
var
  i: integer;
  sum, p: double;
begin
  sum := 0.0;
  for i := 0 to NumXForms - 1 do
    sum := sum + xform[i].Density;
  if sum = 1.0 then Exit;
  for i := 0 to NumXForms - 1 do
  begin
    p := xform[i].Density / sum;
	  if (p < 1E-6) then p := 1E-6;
	  xform[i].Density := p;
  end;
end;

procedure TControlPoint.CalculateWeights;
// AV: Calculate transform weight from its affine determinants
var
  i: smallint;
  deta: double;
begin
  for i := 0 to NumXForms - 1 do
  begin
	  with xform[i] do
      deta := abs(detC * detP);
    if (deta < 1E-6) then deta := 1E-6;

	  xform[i].density := deta;
  end;
end;

procedure TControlPoint.CalculateColorSpeed; // AV: experimental method
var
  i, t: smallint;
  maxw, sumw, w0: double;
  eqw: boolean;
begin
  t := NumXForms;
	if t = 1 then exit; // single xform
  eqw := True;
  w0 := xform[0].density;
  maxw := w0;   sumw := w0;
  for i := 1 to t - 1 do
  begin
    sumw := sumw + xform[i].density;
    if (xform[i].density > maxw) then maxw := xform[i].density;
    if eqw then if (xform[i].density <> w0) then eqw := False;
  end;

  if (not eqw) then
    for i := 0 to t - 1 do
    begin
      with xform[i] do
      begin
        if (GetVariation(0) = 1) then  // if affine - check contraction factor
          symmetry := min(abs(detC), 1)
        else if (density > 0.01) then
          symmetry := 0.85 * density / maxw
        else
          symmetry := -1;
      end;
    end
  else  // equal weights
    for i := 0 to t - 1 do
      with xform[i] do
        symmetry := min(abs(detC * detP), 0.85);
end;

//*********************************************************************//

procedure TControlPoint.GetFromTriangles(const Triangles: TTriangles; const t: integer);
var
  i: integer;
  v: double;
begin
  for i := 0 to t do
    if xform[i].postXswap then
    begin
      solve3(Triangles[-1].x[0], -Triangles[-1].y[0], Triangles[i].x[0],
             Triangles[-1].x[1], -Triangles[-1].y[1], Triangles[i].x[1],
             Triangles[-1].x[2], -Triangles[-1].y[2], Triangles[i].x[2],
             xform[i].p[0][0],    xform[i].p[1][0],   xform[i].p[2][0]);

      solve3(Triangles[-1].x[0], -Triangles[-1].y[0], -Triangles[i].y[0],
             Triangles[-1].x[1], -Triangles[-1].y[1], -Triangles[i].y[1],
             Triangles[-1].x[2], -Triangles[-1].y[2], -Triangles[i].y[2],
             xform[i].p[0][1],    xform[i].p[1][1],    xform[i].p[2][1]);
    end
    else begin
      solve3(Triangles[-1].x[0], -Triangles[-1].y[0], Triangles[i].x[0],
             Triangles[-1].x[1], -Triangles[-1].y[1], Triangles[i].x[1],
             Triangles[-1].x[2], -Triangles[-1].y[2], Triangles[i].x[2],
             xform[i].c[0][0],    xform[i].c[1][0],   xform[i].c[2][0]);

      solve3(Triangles[-1].x[0], -Triangles[-1].y[0], -Triangles[i].y[0],
             Triangles[-1].x[1], -Triangles[-1].y[1], -Triangles[i].y[1],
             Triangles[-1].x[2], -Triangles[-1].y[2], -Triangles[i].y[2],
             xform[i].c[0][1],    xform[i].c[1][1],    xform[i].c[2][1]);
      if xform[i].autoZscale then with xform[i] do begin
        v := c[0][0]*c[1][1] - c[0][1]*c[1][0];
        //n := GetVariationIndex('pre_zscale');
        if v = 1 then
          SetVariation(20, 0.0)  // pre_zscale not needed
        else
          SetVariation(20, sign(v) * sqrt(abs(v)));
      end;
   end;
  FinalXformEnabled := EnableFinalXform;
end;

procedure TControlPoint.GetTriangle(var Triangle: TTriangle; const n: integer);
var
  i, j: integer;
begin
  for i := 0 to 2 do
  with xform[n] do begin
    Triangle.x[i] :=  MainTriangles[-1].x[i] * c[0][0] - MainTriangles[-1].y[i] * c[1][0] + c[2][0];
    Triangle.y[i] := -MainTriangles[-1].x[i] * c[0][1] + MainTriangles[-1].y[i] * c[1][1] - c[2][1];
  end;
end;

procedure TControlPoint.GetPostTriangle(var Triangle: TTriangle; const n: integer);
var
  i, j: integer;
begin
  for i := 0 to 2 do
  with xform[n] do begin
    Triangle.x[i] :=  MainTriangles[-1].x[i] * p[0][0] - MainTriangles[-1].y[i] * p[1][0] + p[2][0];
    Triangle.y[i] := -MainTriangles[-1].x[i] * p[0][1] + MainTriangles[-1].y[i] * p[1][1] - p[2][1];
  end;
end;

////////////////////////////////////////////////////////////////////////////////
procedure TControlPoint.AdjustScale(w, h: integer);
begin
//  if width >= height then
  pixels_per_unit := pixels_per_unit * w/width;
//  else
//    pixels_per_unit := pixels_per_unit * h/height;
  width := w;
  height := h;
end;

end.

