unit NDP;

interface

uses
  Windows, Messages, SysUtils, Classes, Graphics, Controls, Forms, Dialogs,
  StdCtrls,math;

type
  TForm1 = class(TForm)
    Button1: TButton;
    procedure Button1Click(Sender: TObject);
  private
    { Private declarations }
    procedure Draw; 
  public
    { Public declarations }
  end;

const
PointSize=3;
e=10;
PointRight=400;
PointBottom=400;
PointTop=20;
PointLeft=20;
MaxPointsCount=100;
var
  Form1: TForm1;

implementation

{$R *.DFM}

{求兩點(diǎn)P1,P2之間的距離}                      
function Distance(P1,P2:TPoint):real;
begin
result:=sqrt((P1.x-P2.x)*(P1.x-P2.x)+(P1.y-P2.y)*(P1.y-P2.y));
end;

{==================================================

           尋找點(diǎn)集PointType中的最近點(diǎn)對(duì)

           u,v為最近點(diǎn)對(duì)在PointSet中的下標(biāo);
           函數(shù)返回值為P1,P2的距離。
                     written by starfish
                        starfish.h@china.com

===================================================}

function Nearest_Dotted_Pairs(var PointSet:array of TPoint;var u,v:integer):real;
var
X,Y:array of integer;

procedure Swap(var a,b:integer);  {交換a,b}
var
tmp:integer;
begin
tmp:=a;
a:=b;
b:=tmp;
end;

procedure Init_NDP;   {初始化}
var
i,j,Xk,Yk:integer;
begin
setlength(X,Length(PointSet));
setlength(Y,Length(PointSet));
for i:=0 to high(PointSet) do  {初始化數(shù)組}
 begin
  X[i]:=i;
  Y[i]:=i;
 end;

for i:=0 to high(PointSet)-1 do  {對(duì)數(shù)組X和Y進(jìn)行排序，使得X按照點(diǎn)的x坐標(biāo)遞增，Y按照點(diǎn)的y坐標(biāo)遞增}
 begin
   Xk:=i;
   Yk:=i;
   for j:=i+1 to high(PointSet) do
    begin
    if PointSet[X[j]].x<PointSet[X[Xk]].x then Xk:=j;
    if PointSet[Y[j]].y<PointSet[Y[Yk]].y then Yk:=j;
    end;
   Swap(X[i],X[Xk]);
   Swap(Y[i],Y[Yk]);
 end;
end;


{函數(shù)NDP找出X[a..b]中的最近點(diǎn)對(duì)u,v（u,v指示點(diǎn)在PointSet中的下標(biāo)），返回最近點(diǎn)對(duì)的距離；
 數(shù)組Y將X中的點(diǎn)按照y坐標(biāo)遞增排序；X始終保持按照x坐標(biāo)遞增排序}

function NDP(a,b:integer;var u,v:integer;var Y:array of integer):real;
var
Yl,Yr:array of integer;
d1,d2,d3:real;
i,j,t,ul,vl,ur,vr:integer;
begin
if b-a=1 then {如果X[a..b]中只有2個(gè)點(diǎn)，這兩個(gè)點(diǎn)就是最近點(diǎn)對(duì)}
 begin
  result:=distance(PointSet[X[a]],PointSet[X[b]]);
  u:=X[a];
  v:=X[b];
  exit;
 end;
if b-a=2 then  {如果X[a..b]中只有3個(gè)點(diǎn)，就直接求出兩兩之間的距離，找到最近點(diǎn)對(duì)}
 begin
  d1:=distance(PointSet[X[a]],PointSet[X[a+1]]);
  d2:=distance(PointSet[X[a]],PointSet[X[a+2]]);
  d3:=distance(PointSet[X[a+1]],PointSet[X[a+2]]);
  if (d1<=d2)and(d1<=d3) then   {d1最小}
    begin
     u:=X[a];
     v:=X[a+1];
     result:=d1;
    end
   else
    if (d2<=d1)and(d2<=d3) then  {d2最小}
    begin
     u:=X[a];
     v:=X[a+2];
     result:=d2;
    end
    else      {d3最小}
    begin
     u:=X[a+1];
     v:=X[a+2];
     result:=d3;
    end;
   exit;
 end;
{對(duì)X[a..b]進(jìn)行劃分，劃分為X[a..t]和X[t+1..b]，使兩部分盡量平均}
t:=(a+b)div 2;
SetLength(Yl,0);
SetLength(Yr,0);
{將Y分割成Yl和Yr，使得Yl里的點(diǎn)屬于X[a..t],Yr里的點(diǎn)屬于X[t+1..b],并且按照y坐標(biāo)遞增排序}
for i:=0 to high(Y) do
 if PointSet[Y[i]].x<=PointSet[X[t]].x then  {說(shuō)明點(diǎn)Y[i]屬于X[a..t]}
      begin
       SetLength(Yl,Length(Yl)+1);
       Yl[High(Yl)]:=Y[i];
      end
 else begin
       SetLength(Yr,Length(Yr)+1);
       Yr[High(Yr)]:=Y[i];
      end;
{遞歸求出X[a..t]和X[t+1..b]中的最近點(diǎn)對(duì)}
d1:=NDP(a,t,ul,vl,Yl);
d2:=NDP(t+1,b,ur,vr,Yr);
{使d1成為左右兩個(gè)區(qū)間內(nèi)最近點(diǎn)對(duì)的距離，ul,vl為最近點(diǎn)對(duì)}
if d2<d1 then
 begin
   d1:=d2;
   ul:=ur;
   vl:=vr;
 end;
{將Y中以分界線l=PointSet[X[t]].x為中線，寬度為2*d1的帶狀區(qū)域內(nèi)的點(diǎn)存儲(chǔ)在Yl中}
setlength(Yl,0);
for i:=0 to high(Y) do
 if abs(PointSet[Y[i]].x-PointSet[X[t]].x)<=2*d1 then
  begin
   SetLength(Yl,Length(Yl)+1);
   Yl[High(Yl)]:=Y[i];
  end;
{對(duì)于Yl中的每一個(gè)點(diǎn)Yl[i]，計(jì)算Yl中其它點(diǎn)到它的距離，找出比d1小的值；
 根據(jù)抽屜原理知最多只要計(jì)算Yl[i]之后的7個(gè)點(diǎn)就可以了}
for i:=0 to high(Yl)-1 do
 for j:=i+1 to min(i+7,High(Yl)) do
  begin
   d2:=distance(PointSet[Yl[i]],PointSet[Yl[j]]);
   if d2<d1 then
    begin
     d1:=d2;
     ul:=Yl[i];
     vl:=Yl[j];
    end;
  end;
result:=d1;
u:=ul;
v:=vl;
end;


begin
Init_NDP; {初始化}
result:=NDP(0,high(X),u,v,Y);
end;

{==========================================================}

procedure GernateTestData(var PointSet:array of TPoint;count:integer);

var
i,j:integer;
begin

for i:=0 to count-1 do
  repeat
    PointSet[i].x:=trunc(random(Pointright-PointLeft))+PointLeft;
    PointSet[i].y:=trunc(random(Pointbottom-PointTop))+PointTop;
    j:=0;
    while (j<=i-1)and( (abs(PointSet[j].x-PointSet[i].x)>e)or(abs(PointSet[j].y-PointSet[i].y)>e))
     do inc(j);
  until j=i;
end;

procedure TForm1.Draw;

var
PointSet:array of TPoint;
count,len,i:integer;
u,v:integer;
nearest:real;
begin
canvas.Rectangle(PointLeft-5,PointTop-5,Pointright+5,Pointbottom+5);
randomize;
repeat
count:=trunc(random(MaxPointsCount));
until count>3;
setlength(PointSet,count);
GernateTestData(PointSet,count);
for i:=0 to count-1 do canvas.Ellipse(PointSet[i].x-PointSize,PointSet[i].y-PointSize,PointSet[i].x+PointSize,PointSet[i].y+PointSize);

nearest:=Nearest_Dotted_Pairs(PointSet,u,v);
canvas.Pen.Color:=clred;
canvas.Pen.width:=3;
canvas.moveTo(PointSet[u].x,PointSet[u].y);
canvas.lineto(PointSet[v].x,PointSet[v].y);
canvas.Pen.width:=1;
canvas.Pen.Color:=clBlack;
caption:='最近點(diǎn)對(duì)之間距離為：'+FloatToStr(nearest);
end;



procedure TForm1.Button1Click(Sender: TObject);
begin
draw;
end;

end.