Quick demonstration how space partition methods could help to get significant gain. Code uses simple 64x64 grid with cells 32x32. Each cell contains indexes of points belonging to this cell. At first brute-force method counts (doubled) number of point pairs with distance<= 4. The second code piece checks only cells neighbouring to the cell containing examined point. For 64K points gain is about 200 times. Complexity is O(N^2) vs O(N^3/2) (for proper cell size and uniform point distribution). More advanced structures like kd-tree might provide better performance for complex cases too (worse distribution etc).
var
Wdt, Hgt, NPts: Integer;
i, j, k, r, m: Integer;
t: DWord;
sqd, maxsqd, cnt, gx, gy: Integer;
P: TArray<TPoint>;
Grid: array[0..63, 0..63] of TArray<Integer>;
begin
Wdt := 2048;
Hgt := 2048;
NPts := 64 * 1024;
RandSeed := 42;
SetLength(P, NPts);
for i := 0 to NPts - 1 do begin
P[i] := Point(Random(Wdt), Random(Hgt));
Grid[P[i].Y shr 5, P[i].X shr 5] := Grid[P[i].Y shr 5, P[i].X shr 5] + [i];
end;
maxsqd := 16;
cnt := 0;
t := GetTickCount;
for j := 0 to NPts - 1 do
for i := 0 to NPts - 1 do begin
sqd := Sqr(P[i].X - P[j].X) + Sqr(P[i].Y - P[j].Y);
if sqd <= maxsqd then
Inc(cnt);
end;
Memo1.Lines.Add(Format('%d %d', [cnt, GetTickCount - t]));
cnt := 0;
t := GetTickCount;
for i := 0 to NPts - 1 do begin
gx := P[i].X shr 5;
gy := P[i].Y shr 5;
for j := Max(0, gy - 1) to Min(gy + 1, 63) do
for k := Max(0, gx - 1) to Min(gx + 1, 63) do
for m := 0 to High(Grid[j, k]) do begin
r := Grid[j, k, m];
sqd := Sqr(P[i].X - P[r].X) + Sqr(P[i].Y - P[r].Y);
if sqd <= maxsqd then
Inc(cnt);
end;
end;
Memo1.Lines.Add(Format('%d %d', [cnt, GetTickCount - t]));
count time (milliseconds)
115406 11466
115406 62