ASContours.pas

unit ASContours;

{ **************************************************************************** }
{ Rejbrand AlgoSim numerics library: implicit plots and contours               }
{ Copyright © 2017–2025 Andreas Rejbrand                                       }
{ https://english.rejbrand.se/                                                 }
{ **************************************************************************** }

{$DEFINE REALCHECK}
{$WARN SYMBOL_PLATFORM OFF}
{$WARN DUPLICATE_CTOR_DTOR OFF}

interface

uses
  SysUtils, Types, Math, Generics.Defaults, Generics.Collections, GenHelpers,
  DoublePoint, ASNum;

type
  TSampleMeta = record
    xmin,
    xmax,
    ymin,
    ymax: TASR;
    nx,
    ny: Cardinal;
    function xspan: TASR; inline;
    function yspan: TASR; inline;
    function xdelta: TASR; inline;
    function ydelta: TASR; inline;
    function avgspan: TASR; inline;
    function maxspan: TASR; inline;
    function minspan: TASR; inline;
    function Aspect: TASR; inline;
    function Rect: TRectD; inline;
    function ScreenToLogical(const P: TPoint): TASR2;
    function LogicalToScreen(const P: TASR2): TPoint;
  end;

  TSamples = TArray<TASR>;

  TContourMeta = record
    Automatic: Boolean;
    Count: Cardinal;
    Levels: TArray<TASR>;
    Interpolate: Boolean;
    Simplify: Boolean;
  end;

  TContour = TArray<TASR2>;

  TConicRec = record
    a, b, c, d, e, f: TASR;
  end;

const
  ZeroContour: TContourMeta =
    (
      Automatic: False;
      Count: 0;
      Levels: [0.0];
      Interpolate: False;
      Simplify: False;
    );
  ZeroContourI: TContourMeta =
    (
      Automatic: False;
      Count: 0;
      Levels: [0.0];
      Interpolate: True;
      Simplify: False;
    );
  ZeroContourS: TContourMeta =
    (
      Automatic: False;
      Count: 0;
      Levels: [0.0];
      Interpolate: False;
      Simplify: True;
    );
  ZeroContourIS: TContourMeta =
    (
      Automatic: False;
      Count: 0;
      Levels: [0.0];
      Interpolate: True;
      Simplify: True;
    );

  AutoContours: TContourMeta =
    (
      Automatic: True;
      Count: 12;
      Levels: nil;
      Interpolate: False;
      Simplify: False;
    );
  AutoContoursI: TContourMeta =
    (
      Automatic: True;
      Count: 12;
      Levels: nil;
      Interpolate: True;
      Simplify: False;
    );
  AutoContoursS: TContourMeta =
    (
      Automatic: True;
      Count: 12;
      Levels: nil;
      Interpolate: False;
      Simplify: True;
    );
  AutoContoursIS: TContourMeta =
    (
      Automatic: True;
      Count: 12;
      Levels: nil;
      Interpolate: True;
      Simplify: True;
    );


function Contours(const ASampleMeta: TSampleMeta; const ASamples: TArray<TASR>;
  const AContourMeta: TContourMeta): TContour;
procedure SimplifyCurves(const ASampleMeta: TSampleMeta; var AContour: TContour);
procedure LineDetect(const ASampleMeta: TSampleMeta; var AContour: TContour);
procedure CircleDetect(const ASampleMeta: TSampleMeta; const AContour: TContour;
  ACircles: TDictionary<Integer, TCircleRec>);
procedure ConicDetect(const ASampleMeta: TSampleMeta; const AContour: TContour;
  AConics: TDictionary<Integer, TConicRec>);
function ContourRadius(const AContour: TContour): TASR;
function ContourBoundingBox(const AContour: TContour): TRectD;

implementation

type
  TRealPoint = TASR2;

function RP(const X, Y: TASR): TASR2; inline;
begin
  Result := TRealPoint.Create(X, Y);
end;

function ipl1(const A, B: TASR): TASR; inline;
begin
  Result := 2 * (0 - A) / (B - A);
end;

function ipl0(const A, B: TASR): TASR; inline;
begin
  Result := 1;
end;

function Contours(const ASampleMeta: TSampleMeta; const ASamples: TArray<TASR>;
  const AContourMeta: TContourMeta): TContour;

var
  ipl: function(const A, B: TASR): TASR;  
  
type
  TEdge = 0..3;
  TCase = 0..15;

const

  //  Edge indices:
  //
  //        0
  //      ┌───┐
  //    3 │   │ 1
  //      └───┘
  //        2


  StartEdge: array[Boolean] of array[TCase] of TEdge =
    (        {0, 1, 2, 3, 4, 5, 6, 7, 8, 9, A, B, C, D, E, F}
      {False}(0, 3, 2, 3, 0, 3, 0, 3, 3, 0, 0, 0, 3, 2, 3, 0),
      {True }(0, 2, 1, 1, 1, 0, 2, 0, 0, 2, 1, 1, 1, 1, 2, 0)
    );

  OtherEdge: array[TCase] of array[TEdge] of TEdge =
    (     {0, 1, 2, 3}
      { 0}(0, 0, 0, 0),
      { 1}(0, 0, 3, 2),
      { 2}(0, 2, 1, 0),
      { 3}(0, 3, 0, 1),
      { 4}(1, 0, 0, 0),
      { 5}(3, 2, 1, 0),
      { 6}(2, 0, 0, 0),
      { 7}(3, 0, 0, 0),
      { 8}(3, 0, 0, 0),
      { 9}(2, 0, 0, 0),
      {10}(1, 0, 3, 2),
      {11}(1, 0, 0, 0),
      {12}(0, 3, 0, 1),
      {13}(0, 2, 1, 0),
      {14}(0, 0, 3, 2),
      {15}(0, 0, 0, 0)
    );

  OtherEdgeNextCell: array[TCase] of array[TEdge] of TEdge =
    (     {0, 1, 2, 3}
      { 0}(2, 2, 2, 2),
      { 1}(2, 2, 1, 0),
      { 2}(2, 0, 3, 2),
      { 3}(2, 1, 2, 3),
      { 4}(3, 2, 2, 2),
      { 5}(1, 0, 3, 2),
      { 6}(0, 2, 2, 2),
      { 7}(1, 2, 2, 2),
      { 8}(1, 2, 2, 2),
      { 9}(0, 2, 2, 2),
      {10}(3, 2, 1, 0),
      {11}(3, 2, 2, 2),
      {12}(2, 1, 2, 3),
      {13}(2, 0, 3, 2),
      {14}(2, 2, 1, 0),
      {15}(2, 2, 2, 2)
    );

var
  BoxCountX, BoxCountY, VertexCountX, VertexCountY: Integer;
  LogicalBoxSizeX, LogicalBoxSizeY, LogicalHalfBoxSizeX, LogicalHalfBoxSizeY,
    DX, DY: TASR;
  S: TArray<Boolean>;
  R: TArray<TASR>;
  B: TArray<Byte>;

  function NextBox(var AOrigin: Integer; ACode: Byte; var AStartEdge: TEdge): Boolean;
  begin

    ACode := ACode and %1111;

    var x := AOrigin mod BoxCountX;
    var y := AOrigin div BoxCountX;
    var LStartEdge := AStartEdge;

    AStartEdge := OtherEdge[ACode, LStartEdge];     // in case of Result = False

    case ACode of
      1, 14:
        begin
          if LStartEdge = 3 then
            Inc(y)
          else
            Dec(x);
        end;
      2, 13:
        begin
          if LStartEdge = 2 then
            Inc(x)
          else
            Inc(y);
        end;
      3, 12:
        begin
          if LStartEdge = 3 then
            Inc(x)
          else
            Dec(x);
        end;
      4, 11:
        begin
          if LStartEdge = 0 then
            Inc(x)
          else
            Dec(y);
        end;
      5:
        begin
          case LStartEdge of
            0:
              Dec(x);
            1:
              Inc(y);
            2:
              Inc(x);
            3:
              Dec(y);
          end;
        end;
      10:
        begin
          case LStartEdge of
            0:
              Inc(x);
            1:
              Dec(y);
            2:
              Dec(x);
            3:
              Inc(y);
          end;
        end;
      6, 9:
        begin
          if LStartEdge = 0 then
            Inc(y)
          else
            Dec(y);
        end;
      7, 8:
        begin
          if LStartEdge = 3 then
            Dec(y)
          else
            Dec(x);
        end;
    else
      Exit(False);
    end;

    if (x < 0) or (y < 0) or (x > BoxCountX - 1) or (y > BoxCountY - 1) then
      Exit(False);

    Result := True;
    AOrigin := y * BoxCountX + x;
    AStartEdge := OtherEdgeNextCell[ACode, LStartEdge];

  end;

  function GetCoords(i, j: Integer; ACode: Byte; AEdge: TEdge): TRealPoint;
  begin
    var P := ASampleMeta.ScreenToLogical(Point(i, j)); // top-left vertex of box (i, j)
    case ACode and %1111 of
      1, 14 {╲ }:
        begin
          if AEdge = 3 then
            Result := P + RP(0, DY * ipl(R[j * VertexCountX + i], R[Succ(j) * VertexCountX + i]))
          else
            Result := P + RP(DX * ipl(R[Succ(j) * VertexCountX + i], R[Succ(j) * VertexCountX + Succ(i)]), 2*DY);
        end;
      2, 13 { ╱}:
        begin
          if AEdge = 2 then
            Result := P + RP(DX * ipl(R[Succ(j) * VertexCountX + i], R[Succ(j) * VertexCountX + Succ(i)]), 2*DY)
          else
            Result := P + RP(2*DX, DY * ipl(R[j * VertexCountX + Succ(i)], R[Succ(j) * VertexCountX + Succ(i)]));
        end;
      3, 12 {─}:
        begin
          if AEdge = 3 then
            Result := P + RP(0, DY * ipl(R[j * VertexCountX + i], R[Succ(j) * VertexCountX + i]))
          else
            Result := P + RP(2*DX, DY * ipl(R[j * VertexCountX + Succ(i)], R[Succ(j) * VertexCountX + Succ(i)]));
        end;
      4, 11: { ╲}
        begin
          if AEdge = 0 then
            Result := P + RP(DX * ipl(R[j * VertexCountX + i], R[j * VertexCountX + Succ(i)]), 0)
          else
            Result := P + RP(2*DX, DY * ipl(R[j * VertexCountX + Succ(i)], R[Succ(j) * VertexCountX + Succ(i)]));
        end;
      5 {╱╱}:
        begin
          case AEdge of
            0:
              Result := P + RP(DX * ipl(R[j * VertexCountX + i], R[j * VertexCountX + Succ(i)]), 0);
            1:
              Result := P + RP(2*DX, DY * ipl(R[j * VertexCountX + Succ(i)], R[Succ(j) * VertexCountX + Succ(i)]));
            2:
              Result := P + RP(DX * ipl(R[Succ(j) * VertexCountX + i], R[Succ(j) * VertexCountX + Succ(i)]), 2*DY);
            3:
              Result := P + RP(0, DY * ipl(R[j * VertexCountX + i], R[Succ(j) * VertexCountX + i]));
          else
            Result := P + RP(DX, DY);
          end;
        end;
      10 {╲╲}:
        begin
          case AEdge of
            0:
              Result := P + RP(DX * ipl(R[j * VertexCountX + i], R[j * VertexCountX + Succ(i)]), 0);
            1:
              Result := P + RP(2*DX, DY * ipl(R[j * VertexCountX + Succ(i)], R[Succ(j) * VertexCountX + Succ(i)]));
            2:
              Result := P + RP(DX * ipl(R[Succ(j) * VertexCountX + i], R[Succ(j) * VertexCountX + Succ(i)]), 2*DY);
            3:
              Result := P + RP(0, DY * ipl(R[j * VertexCountX + i], R[Succ(j) * VertexCountX + i]))
          else
            Result := P + RP(DX, DY);
          end;
        end;
      6, 9 {│}:
        begin
          if AEdge = 0 then
            Result := P + RP(DX * ipl(R[j * VertexCountX + i], R[j * VertexCountX + Succ(i)]), 0)
          else
            Result := P + RP(DX * ipl(R[Succ(j) * VertexCountX + i], R[Succ(j) * VertexCountX + Succ(i)]), 2*DY);
          end;
      7, 8 {╱ }:
        begin
          if AEdge = 3 then
            Result := P + RP(0, DY * ipl(R[j * VertexCountX + i], R[Succ(j) * VertexCountX + i]))
          else
            Result := P + RP(DX * ipl(R[j * VertexCountX + i], R[j * VertexCountX + Succ(i)]), 0);
        end;
    else
      Result := P + RP(DX, DY);
    end;
  end;

type
  PCurveBuf = ^TCurveBuf;
  TCurveBuf = record
    Data: TArray<TRealPoint>;
    Idx: Integer;
  end;
  TCurveBufs = array[Boolean] of TCurveBuf;

begin

  const MaxCurveCount = 10_000;
  const RegionsOnly = True;

  if AContourMeta.Interpolate then
    ipl := @ipl1
  else
    ipl := @ipl0;  

  var LCurveCount := 0;
  var k := 0;
  try

    var tc1, tc2, freq: Int64;

    VertexCountX := ASampleMeta.nx;
    VertexCountY := ASampleMeta.ny;
    BoxCountX := VertexCountX - 1;
    BoxCountY := VertexCountY - 1;
    LogicalBoxSizeX := ASampleMeta.xspan / BoxCountX;
    LogicalBoxSizeY := ASampleMeta.yspan / BoxCountY;
    LogicalHalfBoxSizeX := LogicalBoxSizeX / 2;
    LogicalHalfBoxSizeY := LogicalBoxSizeY / 2;
    DX :=  LogicalHalfBoxSizeX;
    DY := -LogicalHalfBoxSizeY;

    if VertexCountX * VertexCountY <> Length(ASamples) then
      raise Exception.Create('Wrong number of samples supplied to contour generator.');

    Result := nil;

    if BoxCountX < 2 then
      Exit;

    if BoxCountY < 2 then
      Exit;

    SetLength(Result, 2 * BoxCountX * BoxCountY);

    var LCurveBufs := Default(TCurveBufs);
    SetLength(LCurveBufs[False].Data, 2 * BoxCountX * BoxCountY + 64);
    SetLength(LCurveBufs[True ].Data, 2 * BoxCountX * BoxCountY + 64);

    var LLevels := TArray<TASR>(nil);

    var SampledValues := ASamples;

    if AContourMeta.Automatic then
    begin
      const LevelCount = AContourMeta.Count;
      if LevelCount = 0 then
        Exit(nil);
      var LMin := SampledValues[0];
      var LMax := SampledValues[0];
      for var x in SampledValues do
      begin
        if x > LMax then
          LMax := x;
        if x < LMin then
          LMin := x;
      end;
      if LMin = LMax then
        Exit(nil);
      const LLevelDistance = (LMax - LMin) / (LevelCount + 1);
      SetLength(LLevels, LevelCount);
      for var i := 1 to LevelCount do
        LLevels[i - 1] := LMin + i * LLevelDistance;
    end
    else
      LLevels := AContourMeta.Levels;

    for var LLevel in LLevels do
    begin

      S := nil;
      R := nil;
      B := nil;

      SetLength(S, VertexCountX * VertexCountY);
      SetLength(R, VertexCountX * VertexCountY);

      begin
        var idx := 0;
        for var j := 0 to VertexCountY - 1 do
          for var i := 0 to VertexCountX - 1 do
          begin
            R[idx] := SampledValues[idx] - LLevel;
            if R[idx] > 0 then
              S[idx] := True;
            Inc(idx);
          end;
      end;

      SetLength(B, BoxCountX * BoxCountY);

      for var j := 0 to BoxCountY - 1 do
        for var i := 0 to BoxCountX - 1 do
        begin
          var LCode :=
            Ord(S[     j  * VertexCountX +      i ])    shl 3
              or
            Ord(S[     j  * VertexCountX + Succ(i)])    shl 2
              or
            Ord(S[Succ(j) * VertexCountX + Succ(i)])    shl 1
              or
            Ord(S[Succ(j) * VertexCountX +      i ])    shl 0;
//          if LCode in [5, 10] then
//          begin
//            var P := ASampleMeta.ScreenToLogical(Point(i, j));
//            P.X := P.X + LogicalBoxSizeX / 2;
//            P.Y := P.Y - LogicalBoxSizeY / 2;
//            const FcnValBoxCtr = F(P.X, P.Y) - LLevel;                         // Request additional sample
//            if (LCode = 5) and (FcnValBoxCtr <= 0) then
//              LCode := 10
//            else if (LCode = 10) and (FcnValBoxCtr <= 0) then
//              LCode := 5;
//          end;
          B[j * BoxCountX + i] := LCode;
        end;

      for var j := 0 to BoxCountY - 1 do
        for var i := 0 to BoxCountX - 1 do
        begin
          var BoxIdx := j * BoxCountX + i;
          var Value := B[BoxIdx];
          if Value and 128 <> 0 then
            Continue;
          Value := Value and %1111;
          if Value in [0, 15] then
            Continue;
          // New curve
          Inc(LCurveCount);
          if LCurveCount > MaxCurveCount then
            raise Exception.Create('Maximum curve count exceeded.');
          if k = Length(Result) then
            SetLength(Result, 2 * Succ(Length(Result)));
          Result[k].X := NaN;  // header
          Result[k].Y := 0;    // if closed: +curve length; else -curve length
          const k0 = k;        // idx of header
          Inc(k);
          for var LReverse := False to True do
          begin
            const LCurveBuf: PCurveBuf = @LCurveBufs[LReverse];
            LCurveBuf.Idx := 0;
            var x, y: Integer;
            var idx := BoxIdx;
            var Edge := StartEdge[LReverse, Value and %1111];
            if LReverse and not NextBox(idx, B[idx], Edge) or (B[idx] and 128 <> 0) or (B[idx] in [0, 15]) then
              Break;
            repeat
              x := idx mod BoxCountX;
              y := idx div BoxCountX;
              B[idx] := B[idx] or 128;
              LCurveBuf.Data[LCurveBuf.Idx] := GetCoords(x, y, B[idx], Edge);
              Inc(LCurveBuf.Idx);
              if LCurveBuf.Idx = Length(LCurveBuf.Data) then
                raise Exception.Create('Buffer overflow.');
              if not NextBox(idx, B[idx], Edge) then
                Break;
              if B[idx] and 128 <> 0 then
                Break;
              if B[idx] in [0, 15] then
                Break;
            until False;
            x := idx mod BoxCountX;
            y := idx div BoxCountX;
            LCurveBuf.Data[LCurveBuf.Idx] := GetCoords(x, y, B[idx], Edge);
            Inc(LCurveBuf.Idx);
            if LCurveBuf.Idx = Length(LCurveBuf.Data) then
              raise Exception.Create('Buffer overflow.');
          end;
          var LTotal := LCurveBufs[False].Idx + LCurveBufs[True].Idx;
          if k + LTotal > Length(Result) then
            SetLength(Result, 2 * Succ(Length(Result)));
          if LCurveBufs[True].Idx >= 1 then
            for var q := LCurveBufs[True].Idx - 1 downto 0 do
            begin
              Result[k] := LCurveBufs[True].Data[q];
              Inc(k);
            end;
          if LCurveBufs[False].Idx >= 1 then
          begin
            Move(LCurveBufs[False].Data[0], Result[k], LCurveBufs[False].Idx * SizeOf(TRealPoint));
            Inc(k, LCurveBufs[False].Idx);
          end;
          if ((k - k0) >= 3) and (Result[k0 + 1] = Result[k - 1]) then
            Result[k0].Y := LTotal
          else
            Result[k0].Y := -LTotal
        end;

    end;

  finally
    SetLength(Result, k);          
    if AContourMeta.Simplify then
    begin
      LineDetect(ASampleMeta, Result);
      SimplifyCurves(ASampleMeta, Result);
    end;
  end;

end;

procedure SimplifyCurves(const ASampleMeta: TSampleMeta; var AContour: TContour);

  function SameDirection(const u, v: TRealPoint; const Denom: TASR): Boolean;
  begin
    if u.IsZero or v.IsZero then
      Exit(True);
    const alpha = (u / u.Norm) * (v / v.Norm);
    if SameValue(alpha, 1.0) then
      Exit(True)
    else if SameValue(alpha, -1.0) then
      Exit(False);
    Result := ArcCos(alpha) < Pi / Denom;
  end;

  function SameDirectionSmall(const u, v: TRealPoint): Boolean; inline;
  begin
    Result := SameDirection(u, v, 16);
  end;

  function SameDirectionLarge(const u, v: TRealPoint): Boolean; inline;
  begin
    Result := SameDirection(u, v, 1024);
  end;

  function SameDirectionMed(const u, v: TRealPoint): Boolean; inline;
  begin
    Result := SameDirection(u, v, 64);
  end;

  var LSegmentAccumLength: TASR;

  function AcceptableLengthReduction(const ANewStep, ASuggestedSegmentDisplacement: TRealPoint): Boolean;
  begin
    LSegmentAccumLength := LSegmentAccumLength + ANewStep.Norm;
    if IsZero(LSegmentAccumLength) then
      Exit(False);
    Result := ASuggestedSegmentDisplacement.Norm / LSegmentAccumLength > 0.99999;
  end;

  function SubsegmentCheck(AFirstIndex, ALastIndex: Integer): Boolean;
  begin
    const Span = ALastIndex - AFirstIndex;
    if Span < 2 then
      Exit(True);
    const TotalDisplacement = AContour[ALastIndex] - AContour[AFirstIndex];
    for var n := 2 to 5 do
    begin
      const IdxStepWidth = Span div n;
      const FirstStep = AContour[AFirstIndex + IdxStepWidth] - AContour[AFirstIndex];
      for var i := 2 to n do
      begin
        const ThisStep = AContour[AFirstIndex + i * IdxStepWidth] - AContour[AFirstIndex + Pred(i) * IdxStepWidth];
        if not SameDirectionLarge(ThisStep, FirstStep) then
          Exit(False);
      end;
      if not SameDirectionMed(TotalDisplacement, FirstStep) then
        Exit(False);
    end;
    Result := True;
  end;

begin

  if AContour = nil then
    Exit;

  var LBuf := TArray<TRealPoint>(nil);
  SetLength(LBuf, Length(AContour));
  var LActualLength := 0;
  var LNewCurveStart := -1;
  try
    var LCurveStart := 0;
    while (LCurveStart <= High(AContour)) and IsNan(AContour[LCurveStart].X) do
    begin
      // new curve, first finalise old one:
      if LNewCurveStart <> -1 then
        LBuf[LNewCurveStart].Y := Sign(LBuf[LNewCurveStart].Y) * (LActualLength - LNewCurveStart - 1);
      // now the new curve:
      LNewCurveStart := LActualLength; // start idx of new curve in LBuf
      LBuf[LActualLength] := AContour[LCurveStart]; // copy header
      Inc(LActualLength);
      const LCurveLength = Abs(Round(AContour[LCurveStart].Y));
//      const LClosed = AContour[LCurveStart].Y > 0;
      if LCurveStart + LCurveLength > High(AContour) then
        raise Exception.Create('Incomplete curve.');
      if LCurveLength >= 2 then
      begin
        // add first point of new curve:
        LBuf[LActualLength] := AContour[Succ(LCurveStart)];
        Inc(LActualLength);
        var i := Succ(Succ(LCurveStart)); // second point of new curve, gives us first displacement vector
        while i <= LCurveStart + LCurveLength do
        begin
          var LSegmentFirstIndex := Pred(i);
          var LSegmentFirstPoint := AContour[Pred(i)];
          var LSegmentFirstVector := AContour[i] - AContour[Pred(i)];
          LSegmentAccumLength := LSegmentFirstVector.Norm;
          while
            (i < LCurveStart + LCurveLength)
              and
            SameDirectionSmall(LSegmentFirstVector, AContour[Succ(i)] - AContour[i])
              and
            SameDirectionLarge(LSegmentFirstVector, AContour[Succ(i)] - LSegmentFirstPoint)
              and
            AcceptableLengthReduction(AContour[Succ(i)] - AContour[i], AContour[Succ(i)] - LSegmentFirstPoint)
              and
            SubsegmentCheck(LSegmentFirstIndex, Succ(i))
          do
            Inc(i);
          LBuf[LActualLength] := AContour[i]; // line segment end point (=next segment start point)
          Inc(LActualLength);
          Inc(i);
        end;
      end;
      // while loop epilogue:
      Inc(LCurveStart, LCurveLength + 1);
    end;
    if LNewCurveStart <> -1 then
      LBuf[LNewCurveStart].Y := Sign(LBuf[LNewCurveStart].Y) * (LActualLength - LNewCurveStart - 1);
  finally
    SetLength(LBuf, LActualLength);
    AContour := LBuf;
  end;

end;

procedure LineDetect(const ASampleMeta: TSampleMeta; var AContour: TContour);
type
  TLineRec = record
    P, Q: TRealPoint;
  end;
begin

  if AContour = nil then
    Exit;

  var LLines := TDictionary<Integer, TLineRec>.Create;
  try

    var LCurveStart := 0;
    while (LCurveStart <= High(AContour)) and IsNan(AContour[LCurveStart].X) do
    begin
      const LCurveLength = Abs(Round(AContour[LCurveStart].Y));
      const LClosed = AContour[LCurveStart].Y > 0;
      if LCurveStart + LCurveLength > High(AContour) then
        raise Exception.Create('Incomplete curve.');
      if (LCurveLength >= 2) and not LClosed then
      begin
        const LDisplacement = AContour[LCurveStart + LCurveLength] - AContour[LCurveStart + 1];
        if not LDisplacement.IsZero then
        begin
          const u = LDisplacement / LDisplacement.Norm;
          var LMaxDistance: TASR := 0.0;
          for var i := Succ(0) to Pred(LCurveLength - 1) do
          begin
            const LVector = AContour[LCurveStart + 1 + i] - AContour[LCurveStart + 1];
            if not LVector.IsZero then
            begin
              const LProjection = (LVector * u) * u;
              const LDistance = (LVector - LProjection).Norm;
              if LDistance > LMaxDistance then
                LMaxDistance := LDistance;
            end;
          end;
          if LMaxDistance < ASampleMeta.avgspan / 1000 then
          begin
            var L := Default(TLineRec);
            L.P := AContour[LCurveStart + 1];
            L.Q := AContour[LCurveStart + LCurveLength];
            LLines.Add(LCurveStart, L);
          end;
        end;
      end;
      // while loop epilogue:
      Inc(LCurveStart, LCurveLength + 1);
    end;

    if LLines.Count > 0 then
    begin

      var LBuf := TArray<TRealPoint>(nil);
      SetLength(LBuf, Length(AContour));
      LCurveStart := 0;
      var LActualLength := 0;
      try
        while (LCurveStart <= High(AContour)) and IsNan(AContour[LCurveStart].X) do
        begin
          const LCurveLength = Abs(Round(AContour[LCurveStart].Y));
          var L: TLineRec;
          if LLines.TryGetValue(LCurveStart, L) then
          begin
            if LActualLength + 2 > High(LBuf) then
              raise Exception.Create('Buffer overflow.');
            LBuf[LActualLength].X := NaN;
            LBuf[LActualLength].Y := -2;
            Inc(LActualLength);
            LBuf[LActualLength] := L.P;
            Inc(LActualLength);
            LBuf[LActualLength] := L.Q;
            Inc(LActualLength);
          end
          else
          begin
            if LActualLength + LCurveLength > High(LBuf) then
              raise Exception.Create('Buffer overflow.');
            Move(AContour[LCurveStart], LBuf[LActualLength], (1 + LCurveLength) * SizeOf(TRealPoint));
            Inc(LActualLength, 1 + LCurveLength);
          end;
          // while loop epilogue:
          Inc(LCurveStart, LCurveLength + 1);
        end;
      finally
        SetLength(LBuf, LActualLength);
        AContour := LBuf;
      end;

    end;

  finally
    LLines.Free;
  end;

end;

procedure CircleDetect(const ASampleMeta: TSampleMeta; const AContour: TContour;
  ACircles: TDictionary<Integer, TCircleRec>);
begin

  if AContour = nil then
    Exit;

  if ACircles = nil then
    Exit;

  var LCurveStart := 0;
  while (LCurveStart <= High(AContour)) and IsNan(AContour[LCurveStart].X) do
  begin
    const LCurveLength = Abs(Round(AContour[LCurveStart].Y));
    const LClosed = AContour[LCurveStart].Y > 0;
    if LCurveStart + LCurveLength > High(AContour) then
      raise Exception.Create('Incomplete curve.');
    const N = 16;
    if (LCurveLength >= 128) and LClosed then
    begin
      var LCandidate := True;
      var LOrigins: TArray<TRealPoint>;
      SetLength(LOrigins, N);
      for var i := 0 to N - 1 do
      begin
        const LOffset = LCurveLength div (3*Succ(N));
        const LStride = LCurveLength div 3;
        const A = AContour[LCurveStart + 1 + i*LOffset + 0 * LStride];
        const B = AContour[LCurveStart + 1 + i*LOffset + 1 * LStride];
        const C = AContour[LCurveStart + 1 + i*LOffset + 2 * LStride];
        const z1 = ASC(A.X, A.Y);
        const z2 = ASC(B.X, B.Y);
        const z3 = ASC(C.X, C.Y);
        if SameValue2(z1, z2) then
        begin
          LCandidate := False;
          Break;
        end;
        const w = (z3 - z1) / (z2 - z1);
        if IsZero(w.Im) then
        begin
          LCandidate := False;
          Break;
        end;
        const CplCtr = (z2 - z1) * (w - w.ModSqr) / (2*ImaginaryUnit*w.Im) + z1;
        LOrigins[i].X := CplCtr.Re;
        LOrigins[i].Y := CplCtr.Im;
        if (i > 0) and ((LOrigins[i] - LOrigins[0]).Norm > ASampleMeta.maxspan / 1000) then
        begin
          LCandidate := False;
          Break;
        end;
      end;
      if LCandidate then
      begin
        var LOrigin := Default(TRealPoint);
        for var P in LOrigins do
          LOrigin := LOrigin + (1/N) * P;
        var LRadii: TArray<TASR>;
        SetLength(LRadii, LCurveLength);
        for var i := 0 to LCurveLength - 1 do
          LRadii[i] := (AContour[LCurveStart + 1 + i] - LOrigin).Norm;
        const r = Math.Mean(LRadii);
        if not IsZero(r) then
        begin
          for var i := 0 to LCurveLength  - 1 do
            if Abs((AContour[LCurveStart + 1 + i] - LOrigin).Norm - r) >= ASampleMeta.maxspan / 1000 then
            begin
              LCandidate := False;
              Break;
            end;
        end;
        if LCandidate then
        begin
          var C := Default(TCircleRec);
          C.P := LOrigin;
          C.r := r;
          ACircles.Add(LCurveStart, C);
        end;
      end;
    end;
    // while loop epilogue:
    Inc(LCurveStart, LCurveLength + 1);
  end;

end;

procedure ConicDetect(const ASampleMeta: TSampleMeta; const AContour: TContour;
  AConics: TDictionary<Integer, TConicRec>);
begin

  if AContour = nil then
    Exit;

  if AConics = nil then
    Exit;

  var LCurveStart := 0;
  while (LCurveStart <= High(AContour)) and IsNan(AContour[LCurveStart].X) do
  begin

    const LCurveLength = Abs(Round(AContour[LCurveStart].Y));
    const LClosed = AContour[LCurveStart].Y > 0;
    if LCurveStart + LCurveLength > High(AContour) then
      raise Exception.Create('Incomplete curve.');

    if (LCurveLength >= 32) and LClosed then
    begin

      for var i := LCurveStart + 2 to LCurveStart + LCurveLength do
        if
          (AContour[i] - AContour[Pred(i)]).Norm > ASampleMeta.minspan / 100
        then
          Exit;

      const N = 16;

      var LCandidate := True;
      var LControls: TArray<TRealPoint>;
      SetLength(LControls, 5);
      var LConics: TArray<TRealVector>;
      SetLength(LConics, N);

      for var i := 0 to N - 1 do
      begin
        const LOffset = LCurveLength div (Length(LControls)*Succ(N));
        const LStride = LCurveLength div Length(LControls);
        for var j := 0 to High(LControls) do
          LControls[j] := AContour[LCurveStart + 1 + i * LOffset + j * LStride];

        var A := TRealMatrix.Create(TMatrixSize.Create(5, 5),
          function(i, j: Integer): TASR
          begin
            case j of
              1:
                Result := LControls[i - 1].X * LControls[i - 1].Y;
              2:
                Result := Sqr(LControls[i - 1].Y);
              3:
                Result := LControls[i - 1].X;
              4:
                Result := LControls[i - 1].Y;
              5:
                Result := 1;
            else
              Result := 0;
            end;
          end);

        var Y := -TRealVector.Create([
          Sqr(LControls[0].X),
          Sqr(LControls[1].X),
          Sqr(LControls[2].X),
          Sqr(LControls[3].X),
          Sqr(LControls[4].X)
        ]);

        if not TrySysSolve(A, Y, LConics[i]) then
        begin
          LCandidate := False;
          Break;
        end;

        LConics[i].Insert(0, 1);

        if LConics[i][1]*LConics[i][1] - 4*LConics[i][0]*LConics[i][2] >= 0 then
        begin
          LCandidate := False;
          Break;
        end;

        if (i > 0) and ((LConics[i] - LConics[0]).MaxNorm > LConics[0].MaxNorm / 100) then
        begin
          LCandidate := False;
          Break;
        end;

      end;

      if LCandidate then
      begin
        var LMeanConic := ZeroVector(6);
        for var LConic in LConics do
          LMeanConic := LMeanConic + (1/N) * LConic;
        const LConicMaxNorm = LMeanConic.MaxNorm;
        var LMaxErr := 0.0;
        for var i := 0 to LCurveLength  - 1 do
        begin
          const P = AContour[LCurveStart + 1 + i];
          const LGrad =
            TRealVector.Create(
              [
                2 * LMeanConic[0] * P.X + LMeanConic[1] * P.Y + LMeanConic[3],
                2 * LMeanConic[2] * P.Y + LMeanConic[1] * P.X + LMeanConic[4]
              ]
            );
          const MaxXDev = LGrad[0] * (ASampleMeta.xspan / 10000);
          const MaxYDev = LGrad[1] * (ASampleMeta.yspan / 10000);
          const MaxAbsDev = Abs(MaxXDev) + Abs(MaxYDev);
          var LThisErr := Abs(TRealVector.Create(P.ConicFactor) * LMeanConic);
          if LThisErr > LMaxErr then
            LMaxErr := LThisErr;
          if LThisErr > MaxAbsDev then
          begin
            LCandidate := False;
            Break;
          end;
        end;

        if LCandidate then
        begin
          var C := Default(TConicRec);
          C.a := LMeanConic.Data[0];
          C.b := LMeanConic.Data[1];
          C.c := LMeanConic.Data[2];
          C.d := LMeanConic.Data[3];
          C.e := LMeanConic.Data[4];
          C.f := LMeanConic.Data[5];
          AConics.Add(LCurveStart, C);
        end;

      end;

    end;

    // while loop epilogue:
    Inc(LCurveStart, LCurveLength + 1);

  end;

end;

function ContourRadius(const AContour: TContour): TASR;
begin
  Result := 0.0;
  for var P in AContour do
    if not P.IsNAN then
    begin
      const r = P.Norm;
      if r > Result then
        Result := r;
    end;
end;

function ContourBoundingBox(const AContour: TContour): TRectD;
begin

  var LHasNonNAN := False;
  Result := Default(TRectD);
  for var i := 0 to High(AContour) do
  begin
    if AContour[i].IsNAN then
      Continue;
    if not LHasNonNAN or (AContour[i].X < Result.Left) then
      Result.Left := AContour[i].X;
    if not LHasNonNAN or (AContour[i].X > Result.Right) then
      Result.Right := AContour[i].X;
    if not LHasNonNAN or (AContour[i].Y < Result.Bottom) then
      Result.Bottom := AContour[i].Y;
    if not LHasNonNAN or (AContour[i].Y > Result.Top) then
      Result.Top := AContour[i].Y;
    LHasNonNAN := True;
  end;

end;

{ TSampleMeta }

function TSampleMeta.Aspect: TASR;
begin
  Result := xspan / yspan;
end;

function TSampleMeta.avgspan: TASR;
begin
  Result := (xspan + yspan) / 2;
end;

function TSampleMeta.LogicalToScreen(const P: TASR2): TPoint;
begin
  Result.X := Round( (P.X - xmin) / xdelta);
  Result.Y := Round(-(P.Y - ymax) / ydelta);
end;

function TSampleMeta.maxspan: TASR;
begin
  Result := Max(xspan, yspan);
end;

function TSampleMeta.minspan: TASR;
begin
  Result := Min(xspan, yspan);
end;

function TSampleMeta.Rect: TRectD;
begin
  Result := TRectD.Create(xmin, ymax, xmax, ymin);
end;

function TSampleMeta.ScreenToLogical(const P: TPoint): TASR2;
begin
  Result.X := xmin + xdelta * P.X;
  Result.Y := ymax - ydelta * P.Y;
end;

function TSampleMeta.xdelta: TASR;
begin
  Result := xspan / (nx - 1);
end;

function TSampleMeta.xspan: TASR;
begin
  Result := xmax - xmin;
end;

function TSampleMeta.ydelta: TASR;
begin
  Result := yspan / (ny - 1);
end;

function TSampleMeta.yspan: TASR;
begin
  Result := ymax - ymin;
end;

end.