Function Version

{==============================================================================]
 Splits TPA with minDist, maxDist.
[==============================================================================}
function TPADistributeEx(TPA: TPointArray; minDist, maxDist: Extended): T2DPointArray;
var
  h, i, l, c, s, x, y, o, r, d, m: Integer;
  p: T2DIntegerArray;
  q: TPointArray;
  a, b: TBox;
  e: Extended;
  z: TPoint;
begin
  SetLength(Result, 0);
  h := High(TPA);
  if ((minDist <= maxDist) and (h > -1)) then
    if (h > 0) then // We make sure array is bigger than 1 item...
    begin
      b.X1 := TPA[0].X; // START: Capturing bounds of TPA
      b.Y1 := TPA[0].Y;
      b.X2 := TPA[0].X;
      b.Y2 := TPA[0].Y;
      r := 0;
      for i := 1 to h do
      begin
        if (TPA[i].X < b.X1) then
          b.X1 := TPA[i].X
        else
          if (TPA[i].X > b.X2) then
            b.X2 := TPA[i].X;
        if (TPA[i].Y < b.Y1) then
          b.Y1 := TPA[i].Y
        else
          if (TPA[i].Y > b.Y2) then
            b.Y2 := TPA[i].Y;
      end; // END: Bounds captured!
      SetLength(p, ((b.X2 - b.X1) + 1)); // We create 2D Integer table, for storing point counts (if array contains duplicates). This part is the 1D.
      for i := 0 to (b.X2 - b.X1) do // Creating 2D side with loop.
      begin
        SetLength(p[i], ((b.Y2 - b.Y1) + 1));
        for c := 0 to (b.Y2 - b.Y1) do
          p[i][c] := 0; // Setting 0 as default value for all array items.
      end;
      if (maxDist < 0.0) then // Making sure maxDist isn't lower than 0
        maxDist := 0.0; // ..when it isn't, we set it to 0.
      d := Ceil(maxDist); // Storing distance by maxDist to d variable, which will be used for scanning the point groups to result.
      m := Max(((b.X2 - b.X1) + 1), ((b.Y2 - b.Y1) + 1)); // We get the maximum radius that we should be using.
      if (d > m) then // We make sure d isn't greater than maximum radius could be in the first place.. To reduce checking a bit, when user hits with long distances that aren't even needed.
        d := m;
      for i := 0 to h do // We loop through whole array increasing the Integer table values with points, this we will get the point counts.
        Inc(p[(TPA[i].X - b.X1)][(TPA[i].Y - b.Y1)]);
      for i := 0 to h do
        if (p[(TPA[i].X - b.X1)][(TPA[i].Y - b.Y1)] > 0) then // Time to start capturing the groups
        begin
          c := Length(Result);
          SetLength(Result, (c + 1)); // Creating first/next group of points for result.
          SetLength(Result[c], p[(TPA[i].X - b.X1)][(TPA[i].Y - b.Y1)]); // We set the length for all points at this coordinate, by point count.
          for o := 0 to (p[(TPA[i].X - b.X1)][(TPA[i].Y - b.Y1)] - 1) do // Loop for planting all the points from current TPA index to this new group.
            Result[c][o] := TPA[i];
          SetLength(q, 1); // Creating new queue
          q[0] := TPA[i]; // ..we set the current TPA item as the object for it...
          r := (r + p[(TPA[i].X - b.X1)][(TPA[i].Y - b.Y1)]); // R variable is for keeping track of points that we have "met" so far, this part increases it by the total points
          if (r > h) then // We can simply exit, there's no more to go through, because there was so many duplicates in the array, probably.
            Exit;
          p[(TPA[i].X - b.X1)][(TPA[i].Y - b.Y1)] := 0; // Setting the point count to 0, that way we don't hit the same place twice, which would fuck up things pretty bad!
          s := 1; // Simply needed for the while..do loop to work.
          while (s > 0) do // While there's still queue items to check through...
          begin
            s := High(q); // Highest queue
            z := q[s]; // Storing the point to Z variable for checking.
            a.X1 := (z.X - d); // Creating the area from point Z (by QUEUE) with radius, for scanning the points.
            a.Y1 := (z.Y - d);
            a.X2 := (z.X + d);
            a.Y2 := (z.Y + d);
            SetLength(q, s); // We decrease queue by 1, pretty much.
            if (a.X1 < b.X1) then // START: Lets constrain TBox a (scan area) to the boundaries of b, which is the variable for our TPA bounds.
              a.X1 := b.X1
            else
              if (a.X1 > b.X2) then
                a.X1 := b.X2;
            if (a.Y1 < b.Y1) then
              a.Y1 := b.Y1
            else
              if (a.Y1 > b.Y2) then
                a.Y1 := b.Y2;
            if (a.X2 < b.X1) then
              a.X2 := b.X1
            else
              if (a.X2 > b.X2) then
                a.X2 := b.X2;
            if (a.Y2 < b.Y1) then
              a.Y2 := b.Y1
            else
              if (a.Y2 > b.Y2) then
                a.Y2 := b.Y2; // END: DONE! This way, there's no looping outside the bounds, speeds up things.
            for x := a.X1 to a.X2 do // Lets loop through area for new points to group.
              for y := a.Y1 to a.Y2 do
                if (p[(x - b.X1)][(y - b.Y1)] > 0) then // When we hit Integer table slot which contains more than 0 points, we possibly hitted a match... Time to capture the point(s)!
                begin
                  e := Sqrt(Sqr(z.X - x) + Sqr(z.Y - y)); // Calculating the distance between Z (center point) and this new point from table.
                  if ((e >= minDist) and (e <= maxDist)) then // ..and if it's within minDist and maxDist, we take it in, else we will leave it, because its not valid.
                  begin
                    l := Length(Result[c]);
                    SetLength(Result[c], (l + p[(x - b.X1)][(y - b.Y1)])); // Once again, capturing the point(s) to group, pretty much the same actions as with the first group item.
                    for o := 0 to (p[(x - b.X1)][(y - b.Y1)] - 1) do
                    begin
                      Result[c][(l + o)].X := x;
                      Result[c][(l + o)].Y := y;
                    end;
                    r := (r + p[(x - b.X1)][(y - b.Y1)]);
                    if (r > h) then
                      Exit;
                    p[(x - b.X1)][(y - b.Y1)] := 0;
                    SetLength(q, (s + 1));
                    q[s] := Result[c][l]; // Added the point to checking/scanning queu.
                    Inc(s);
                  end;
                end;
          end;
        end;
    end else
    begin // When we have only 1 point, there's no need to do any looping, so lets just create the array simply this easy way...
      SetLength(Result, 1);
      SetLength(Result[0], 1);
      Result[0][0] := TPA[0];
    end;
end;

function TPADistribute(TPA: TPointArray; dist: Extended): T2DPointArray;
begin
  Result := TPADistributeEx(TPA, 0, dist);
end;