procedure TPAClassify(TPA: TPointArray; dist: Integer; var output: T2DPointArray); callconv
var
  h, i, l, c, s, x, y, o: Integer;
  p: T2DIntegerArray;
  q: TPointArray;
  a, b: TBox;
begin
  SetLength(output, 0);
  h := High(TPA);
  if (h > -1) then
  begin
    b := TPABounds(TPA);
    CreateATIA(0, ((b.X2 - b.X1) + 1), ((b.Y2 - b.Y1) + 1), p);
    EnsureRange(dist, 0, Max(((b.X2 - b.X1) + 1), ((b.Y2 - b.Y1) + 1)));
    dist := (dist + 1);
    for i := 0 to h do
      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
      begin
        c := Length(output);
        SetLength(output, (c + 1));
        SetLength(output[c], p[(TPA[i].X - b.X1)][(TPA[i].Y - b.Y1)]);
        for o := 0 to (p[(TPA[i].X - b.X1)][(TPA[i].Y - b.Y1)] - 1) do
          output[c][o] := TPA[i];
        SetLength(q, 1);
        q[0] := TPA[i];
        p[(TPA[i].X - b.X1)][(TPA[i].Y - b.Y1)] := 0;
        s := 1;
        while (s > 0) do
        begin
          s := High(q);
          SetLength(q, s);
          l := High(output[c]);
          a.X1 := (output[c][l].X - dist);
          a.Y1 := (output[c][l].Y - dist);
          a.X2 := (output[c][l].X + dist);
          a.Y2 := (output[c][l].Y + dist);
          BoxConstraint(a, b);
          for x := a.X1 to a.X2 do
            for y := a.Y1 to a.Y2 do
              if (p[(x - b.X1)][(y - b.Y1)] > 0) then
              begin
                l := Length(output[c]);
                SetLength(output[c], (l + p[(x - b.X1)][(y - b.Y1)]));
                for o := 0 to (p[(x - b.X1)][(y - b.Y1)] - 1) do
                begin
                  output[c][(l + o)].X := x;
                  output[c][(l + o)].Y := y;
                end;
                p[(x - b.X1)][(y - b.Y1)] := 0;
                SetLength(q, (s + 1));
                q[s] := output[c][l];
                Inc(s);
              end;
        end;
      end;
  end;
end;