ndp.pas

用pascal实现寻找最近点对问题

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}

{求两点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;

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

           寻找点集PointType中的最近点对

           u,v为最近点对在PointSet中的下标;
           函数返回值为P1,P2的距离。
                     written by starfish
                        [e-mail]

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

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  {初始化数组}
 begin
  X[i]:=i;
  Y[i]:=i;
 end;

for i:=0 to high(PointSet)-1 do  {对数组X和Y进行排序，使得X按照点的x坐标递增，Y按照点的y坐标递增}
 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;


{函数NDP找出X[a..b]中的最近点对u,v（u,v指示点在PointSet中的下标），返回最近点对的距离；
 数组Y将X中的点按照y坐标递增排序；X始终保持按照x坐标递增排序}

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个点，这两个点就是最近点对}
 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个点，就直接求出两两之间的距离，找到最近点对}
 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;
{对X[a..b]进行划分，划分为X[a..t]和X[t+1..b]，使两部分尽量平均}
t:=(a+b)div 2;
SetLength(Yl,0);
SetLength(Yr,0);
{将Y分割成Yl和Yr，使得Yl里的点属于X[a..t],Yr里的点属于X[t+1..b],并且按照y坐标递增排序}
for i:=0 to high(Y) do
 if PointSet[Y[i]].x<=PointSet[X[t]].x then  {说明点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]中的最近点对}
d1:=NDP(a,t,ul,vl,Yl);
d2:=NDP(t+1,b,ur,vr,Yr);
{使d1成为左右两个区间内最近点对的距离，ul,vl为最近点对}
if d2<d1 then
 begin
   d1:=d2;
   ul:=ur;
   vl:=vr;
 end;
{将Y中以分界线l=PointSet[X[t]].x为中线，宽度为2*d1的带状区域内的点存储在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;
{对于Yl中的每一个点Yl[i]，计算Yl中其它点到它的距离，找出比d1小的值；
 根据抽屉原理知最多只要计算Yl[i]之后的7个点就可以了}
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:='最近点对之间距离为：'+FloatToStr(nearest);
end;



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

end.