{==============================================================================]
@action: Sorts TPA using Heatmap methods for analyzing/detecting "hot" & "cold" points in TPA.
         Hot points = a lot of neighbour points | Cold points = less/none neighbour points.
@note: Supports custom order with TSortOrder!
@contributors: Janilabo, slacky
[==============================================================================}
procedure TPASortByHeatmapEx(var TPA: TPointArray; order: TSortOrder); callconv
  procedure ShellSort(var arr: TPointArray; heatmap: T2DExtendedArray; area: TBox; order: TSortOrder);
  var
    x, a, b, l: Integer;
    t: TPoint;
  begin
    l := Length(arr);
    if (l > 1) then
    begin
      x := 0;
      while (x < (l div 3)) do
        x := ((x * 3) + 1);
      case order of
        so_HighToLow:
        if (l > 2) then
        begin
          while (x >= 1) do
          begin
            for a := x to (l - 1) do
            begin
              b := a;
              while ((b >= x) and (heatmap[(arr[b].X - area.X1)][(arr[b].Y - area.Y1)] > heatmap[(arr[(b - x)].X - area.X1)][(arr[(b - x)].Y - area.Y1)])) do
              begin
                t := arr[b];
                arr[b] := arr[(b - x)];
                arr[(b - x)] := t;
                b := (b - x);
              end;
            end;
            x := (x div 3);
          end;
        end else
          if (heatmap[(arr[0].X - area.X1)][(arr[0].Y - area.Y1)] < heatmap[(arr[1].X - area.X1)][(arr[1].Y - area.Y1)]) then
          begin
            t := arr[0];
            arr[0] := arr[1];
            arr[1] := t;
          end;
        so_LowToHigh:
        if (l > 2) then
        begin
          while (x >= 1) do
          begin
            for a := x to (l - 1) do
            begin
              b := a;
              while ((b >= x) and (heatmap[(arr[b].X - area.X1)][(arr[b].Y - area.Y1)] < heatmap[(arr[(b - x)].X - area.X1)][(arr[(b - x)].Y - area.Y1)])) do
              begin
                t := arr[b];
                arr[b] := arr[(b - x)];
                arr[(b - x)] := t;
                b := (b - x);
              end;
            end;
            x := (x div 3);
          end;
        end else
          if (heatmap[(arr[0].X - area.X1)][(arr[0].Y - area.Y1)] > heatmap[(arr[1].X - area.X1)][(arr[1].Y - area.Y1)]) then
          begin
            t := arr[0];
            arr[0] := arr[1];
            arr[1] := t;
          end;
      end;
    end;
  end;
var
  x, y, l, i, s: Integer;
  c: T2DIntegerArray;
  p: T2DExtendedArray;
  t: TPointArray;
  b: TBox;
  d: Extended;
begin
  l := Length(TPA);
  if (l > 1) then
  begin
    b := TPABounds(TPA);
    SetLength(p, ((b.X2 - b.X1) + 1));
    SetLength(c, ((b.X2 - b.X1) + 1));
    for x := 0 to (b.X2 - b.X1) do
    begin
      SetLength(p[x], ((b.Y2 - b.Y1) + 1));
      SetLength(c[x], ((b.Y2 - b.Y1) + 1));
      for y := 0 to (b.Y2 - b.Y1) do
      begin
        p[x][y] := 0.0;
        c[x][y] := 0;
      end;
    end;
    SetLength(t, l);
    s := 0;
    for i := 0 to (l - 1) do
    begin
      if (c[(TPA[i].X - b.X1)][(TPA[i].Y - b.Y1)] = 0) then
      begin
        t[s] := TPA[i];
        s := (s + 1);
        p[(TPA[i].X - b.X1)][(TPA[i].Y - b.Y1)] := 1;
      end;
      Inc(c[(TPA[i].X - b.X1)][(TPA[i].Y - b.Y1)]);
    end;
    SetLength(t, s);
    for i := 0 to (s - 1) do
      for x := ((t[i].X - b.X1) - 5) to ((t[i].X - b.X1) + 5) do
        for y := ((t[i].Y - b.Y1) - 5) to ((t[i].Y - b.Y1) + 5) do
          if ((x >= 0) and (y >= 0) and (x <= (b.X2 - b.X1)) and (y <= (b.Y2 - b.Y1))) then
            if (c[x][y] > 0) then
              if ((x <> (t[i].X - b.X1)) or (y <> (t[i].Y - b.Y1))) then
              begin
                d := Sqrt(Sqr((t[i].X - b.X1) - x) + Sqr((t[i].Y - b.Y1) - y));
                if (d <= 5.5) then
                  p[(t[i].X - b.X1)][(t[i].Y - b.Y1)] := (p[(t[i].X - b.X1)][(t[i].Y - b.Y1)] + (0.5 / d));
              end else
                p[(t[i].X - b.X1)][(t[i].Y - b.Y1)] := (p[(t[i].X - b.X1)][(t[i].Y - b.Y1)] + 1);
    ShellSort(t, p, b, order);
    l := 0;
    for x := 0 to (s - 1) do
    begin
      for y := 0 to (c[(t[x].X - b.X1)][(t[x].Y - b.Y1)] - 1) do
        TPA[(l + y)] := t[x];
      l := (l + c[(t[x].X - b.X1)][(t[x].Y - b.Y1)]);
    end;
  end;
end;