距离最近的点对

韩冰

<>给定n 个点（xi，yi）（1≤i≤n），要求找出其中距离最近的两个点。

( M5 b, L$ Y- e- |例14-7 假设在一片金属上钻n 个大小一样的洞，如果洞太近，金属可能会断。若知道任意两个洞的最小距离，可估计金属断裂的概率。这种最小距离问题实际上也就是距离最近的点对问题。
$ k% A# |( u) G! }4 @$ E( u9 d
, ~! {4 N: Z: p/ J+ M& [1 p! g通过检查所有的n(n- 1 ) / 2对点，并计算每一对点的距离，可以找出距离最近的一对点。这种方法所需要的时间为(n2 )。我们称这种方法为直接方法。图1 4 - 1 3中给出了分而治之求解算法的伪代码。该算法对于小的问题采用直接方法求解，而对于大的问题则首先把它划分为两个较小的问题，其中一个问题（称为A）的大小为「n /2ù，另一个问题（称为B）的大小为「n /2ù。初始时，最近的点对可能属于如下三种情形之一： 1) 两点都在A中（即最近的点对落在A中）；2) 两点都在B中；3) 一点在A，一点在B。假定根据这三种情况来确定最近点对，则最近点对是所有三种情况中距离最小的一对点。在第一种情况下可对A进行递归求解，而在第二种情况下可对B进行递归求解。

( E( y4 }) c3 z3 u. J9 C
if (n较小) {用直接法寻找最近点对

R e t u r n ; }

// n较大
0 a  h5 v$ W5 a$ @! a
# Z0 [1 ~; h. X( L2 s  ^5 e5 n+ a" O将点集分成大致相等的两个部分A和B
  Z( R$ F" C) D, q/ X7 G$ k2 A) A
确定A和B中的最近点对

确定一点在A中、另一点在B中的最近点对

从上面得到的三对点中，找出距离最小的一对点
' V. F/ ?- j3 T5 e* b
( E" Q2 x# u; T7 N6 s图14-13 寻找最近的点对
  Z6 R6 e% s2 H- S9 Z3 Z, q$ Q
9 R, U- a0 Z- O' ^4 S/ |
为了确定第三种情况下的最近点对，需要采用一种不同的方法。这种方法取决于点集是如何被划分成A、B的。一个合理的划分方法是从xi（中间值）处划一条垂线，线左边的点属于A，线右边的点属于B。位于垂线上的点可在A和B之间分配，以便满足A、B的大小。

0 K+ L: v* z! P( r例2-8 考察图14-14a 中从a到n的1 4个点。这些点标绘在图14-14b 中。中点xi = 1，垂线x = 1如图14-14b 中的虚线所示。虚线左边的点(如b, c, h, n, i)属于A，右边的点(如a, e, f, j, k, l) 属于B。d, g, m 落在垂线上，可将其中两个加入A, 另一个加入B，以便A、B中包含相同的点数。假设d ,m加入A，g加入B。

' ^% J1 d" e# r: G  S% k# y设是i 的最近点对和B的最近点对中距离较小的一对点。若第三种情况下的最近点对比小。则每一个点距垂线的距离必小于，这样，就可以淘汰那些距垂线距离≥ 的点。图1 4 - 1 5中的虚线是分割线。阴影部分以分割线为中线，宽为2 。边界线及其以外的点均被淘汰掉，只有阴影中的点被保留下来，以便确定是否存在第三类点对（对应于第三种情况）其距离小于。用RA、RB 分别表示A和B中剩下的点。如果存在点对(p,q)，p?A, q?B且p, q 的距离小于，则p?RA，q?RB。可以通过每次检查RA 中一个点来寻找这样的点对。假设考察RA 中的p 点，p的y 坐标为p.y，那么只需检查RB 中满足p.y- ＜q.y＜p.y+ 的q 点，看是否存在与p 间距小于的点。在图14-16a 中给出了包含这种q 点的RB 的范围。因此，只需将RB 中位于×2 阴影内的点逐个与p 配对，以判断p 是否是距离小于的第三类点。这个×2 区域被称为是p 的比较区（comparing region）。
& h9 u' M0 s+ H% C2 |: F
! p( \" }  h  j, ^) P& O; a. \例2-9 考察例2 - 8中的1 4个点。A中的最近点对为(b,h)，其距离约为0 . 3 1 6。B中最近点对为(f, j)，其距离为0 . 3，因此= 0 . 3。当考察是否存在第三类点时，除d, g, i, l, m 以外的点均被淘汰，因为它们距分割线x= 1的距离≥ 。RA ={d, i, m}，RB= {g, l}，由于d 和m 的比较区中没有点，只需考察i即可。i 的比较区中仅含点l。计算i 和l的距离，发现它小于，因此(i, l) 是最近的点对。
9 N3 e, e7 T5 i. l
) b6 z5 ~  A% H8 d为了确定一个距离更小的第三类点，RA 中的每个点最多只需和RB 中的6个点比较，如图1 4 - 1 6所示。

( \% X0 ]. d0 V. O0 {- g/ B1. 选择数据结构
! o5 D- g3 p  ]# `& H
7 |& B7 V( H: [9 }! d为了实现图1 4 - 1 3的分而治之算法，需要确定什么是“小问题”以及如何表示点。由于集合中少于两点时不存在最近点对，因此必须保证分解过程不会产生少于两点的点集。如果将少于四点的点集做为“小问题”，就可以避免产生少于两点的点集。

5 z. I% s' O$ i# W: e# V. k! k3 I每个点可有三个参数：标号， x 坐标，y 坐标。假设标号为整数，每个点可用P o i n t l类（见程序1 4 - 8）来表示。为了便于按x 坐标对各个点排序，可重载操作符＜=。归并排序程序如1 4 -3所示。

程序14-8 点类
" g/ b" e. N+ }5 Q4 u: q/ O
  F* N& Y9 p6 i: oclass Point1 {

friend float dist(const Point1&amp;, const Point1&amp;);

" F! \. h, |6 I3 nfriend void close(Point1 *, Point2 *, Point2 *, int, int, Point1&amp;, Point1&amp;, float&amp;);

friend bool closest(Point1 *, int, Point1&amp;, Point1&amp;,float&amp;);

friend void main();
7 o& q# [1 s, R* K% \
1 u, D5 M1 R# g* Z5 X9 u4 w" cp u b l i c :

int operator&lt;=(Point1 a) const
5 u) u: i  `& B4 X
{return (x &lt;= a.x);}
: n" V6 F0 E: z( N% _, ~3 K
* l9 O! Y/ o, j8 c5 E0 gp r i v a t e :

int ID; // 点的编号
8 i- S! g% \% U2 U8 J
float x, y; // 点坐标
6 i" }% P$ e& [" A) B
( B% z2 h5 L- O# m, Z" ]3 h} ;

class Point2 {

6 c/ e1 K6 t& x% I9 E$ i6 Cfriend float dist(const Point2&amp;, const Point2&amp;);
4 `7 H. w/ w; @) O/ v' M! N7 ~
friend void close(Point1 *, Point2 *, Point2 *, int, int, Point1&amp;, Point1&amp;, float&amp;);
# Q) M: t# ?; }) m/ K/ j
& b1 f5 L) f/ _+ m2 {) ?friend bool closest(Point1 *, int, Point1&amp;, Point1&amp;, float&amp;);

! B$ C3 Z: p4 K6 p! ^3 F- |. zfriend void main();
$ k0 G1 E8 r/ E& a6 B' R" c+ z
1 L+ P, C4 M! Op u b l i c :
( w. ~) N) T5 V# {7 o* O% G
int operator&lt;=(Point2 a) const

{return (y &lt;= a.y);}
* s( r! M/ D& o4 G6 \' @
: F+ o5 w* }, S/ {6 Ip r i v a t e :

int p; // 数组X中相同点的索引
  |( s4 t8 X5 Y$ l. e
float x, y; // 点坐标
" Z+ z7 ], p# _% @) T
; V5 t* A+ u! `6 ]* i2 ~} ;
0 H$ b3 j7 B) A; y
所输入的n 个点可以用数组X来表示。假设X中的点已按照x 坐标排序，在分割过程中如果当前考察的点是X [l :r]，那么首先计算m= (l+r) / 2，X[ l:m]中的点属于A，剩下的点属于B。计算出A和B中的最近点对之后，还需要计算RA 和RB，然后确定是否存在更近的点对，其中一点属于RA，另一点属于RB。如果点已按y 坐标排序，那么可以用一种很简单的方式来测试图1 4 - 1 6。按y 坐标排序的点保存在另一个使用类P o i n t 2 (见程序14-8) 的数组中。注意到在P o i n t 2类中，为了便于y 坐标排序，已重载了操作符＜＝。成员p 用于指向X中的对应点。

8 H" a+ E' l) c+ N3 w' j确定了必要的数据结构之后，再来看看所要产生的代码。首先定义一个模板函数d i s t (见程序1 4 - 9 )来计算点a, b 之间的距离。T可能是P o i n t 1或P o i n t 2，因此d i s t必须是P o i n t 1和P o i n t 2类的友元。

程序14-9 计算两点距离

template<CLASS T>

inline float dist(const T&amp; u, const T&amp; v)
8 w/ f% ^% U* K% ]- d/ z  @
4 i9 [- k3 D, S- _{ / /计算点u 和v之间的距离
( d% k: f5 ?3 A% v9 h- ]
& J: w( E& a% }) V7 k6 {8 ?/ Bfloat dx = u.x-v. x ;
5 M/ X9 X: C7 ~& b* E
4 k5 \3 B! ^" K3 rfloat dy = u.y-v. y ;
, G# m( }2 \8 \2 f( M# P& V
( q4 @  z7 A% ^  F1 O4 w- Nreturn sqrt(dx * dx + dy * dy);
5 K& |/ S3 Y0 F# X1 \8 t
}
  f8 n6 v9 ~, ^9 d
7 u4 N1 I( b% i, `如果点的数目少于两个，则函数c l o s e s t (见程序1 4 - 1 0 )返回f a l s e，如果成功时函数返回t r u e。当函数成功时，在参数a 和b 中返回距离最近的两个点，在参数d 中返回距离。代码首先验证至少存在两点，然后使用M e rg e S o r t函数(见程序14-3) 按x 坐标对X中的点排序。接下来把这些点复制到数组Y中并按y 坐标进行排序。排序完成时，对任一个i，有Y [i ] . y≤Y [i+ 1 ] . y，并且Y [i ] .p给出了点i 在X中的位置。上述准备工作做完以后，调用函数close (见程序1 4 - 11 )，该函数实际求解最近点对。
( U: y0 D7 r# D2 e$ a2 W
$ P2 n3 X! x7 |" A0 a$ c& X  V程序14-10 预处理及调用c l o s e

bool closest(Point1 X[], int n, Point1&amp; a, Point1&amp; b, float&amp; d)

{// 在n &gt;= 2 个点中寻找最近点对
/ A2 l7 z" v- H9 ^7 J9 ^+ s
// 如果少于2个点，则返回f a l s e

: n1 ?8 f2 D, Q8 H# M* ^  [// 否则，在a 和b中返回距离最近的两个点
  [7 x- J, s8 j9 x6 g/ \! ?
* r* I" p% ^4 w- g2 }8 y3 @; xif (n &lt; 2) return false;

// 按x坐标排序

2 Y; e9 o" \- J1 L$ nM e r g e S o r t ( X , n ) ;
; O3 V/ `6 m! m) G) u: g
: k( {! t/ Y* O' D  ?2 ?# r3 V// 创建一个按y坐标排序的点数组

; u; G( S  t3 V. yPoint2 *Y = new Point2 [n];
; J7 j8 u1 z# V3 R0 T9 R4 ~
& j8 _  G. p" S. i! _- A* y5 ]for (int i = 0; i &lt; n; i++) {
+ j4 d4 o3 k/ T+ G( L" g, z
// 将点i 从X 复制到Y
6 o) i4 ~& Z  ]- l/ _, I6 `
& o! X  V2 X, {6 a# n9 YY.p = i;

( j( K0 b$ H1 b" BY.x = X.x;

Y.y = X.y;

}
6 \$ ^# b8 G+ ?* n) G. x* D) Y& S2 g. e
M e r g e S o r t ( Y,n); // 按y坐标排序

: K' l) D7 f  r// 创建临时数组

/ {! Z* j) u7 v1 Z0 [Point2 *Z = new Point2 [n];

2 i$ n5 P- t8 T, V8 F* s9 p// 寻找最近点对

c l o s e ( X , Y, Z , 0 , n - 1 , a , b , d ) ;
# Q( S4 Y. h) P. X  k
// 删除数组并返回
- p& m) {4 w/ _+ n1 R# N8 h
6 V; D4 J1 [- I/ C) Gdelete [] Y;
5 G; W" _# G4 X! |- p) e$ B- s
delete [] Z;
6 j% Q9 x# r6 }+ h
4 g# d6 ~$ o$ c& u" r% B3 U# Vreturn true;

}

6 v% c" Z9 v+ l5 j) R+ x程序1 4 - 11 计算最近点对

# `" U( \+ H+ t. mvoid close(Point1 X[], Point2 Y[], Point2 Z[], int l, int r, Point1&amp; a, Point1&amp; b, float&amp; d)
/ n: U! S8 o6 }  w, ~
& K1 }+ b/ v' P) l4 Y% k{//X[l:r] 按x坐标排序
- {- m; b/ }. H& J+ H3 c9 y
+ [0 ^3 _/ |) F# P7 z1 x6 f//Y[l:r] 按y坐标排序

if (r-l == 1) {// 两个点
! G& d/ U6 d- F/ q6 X9 {2 v" @
a = X[l];

- G2 t! E& z/ F! l+ {$ k7 G% s! Fb = X[r];

8 k2 U/ |+ c1 A3 Jd = dist(X[l], X[r]);

/ S) y  ~2 h+ r8 Pr e t u r n ; }

2 L  |; i9 t  ^. K1 _/ vif (r-l == 2) {// 三个点

// 计算所有点对之间的距离

$ ~5 |' k4 p. }/ d, E! [' Ifloat d1 = dist(X[l], X[l+1]);

& _& N0 m# h' O, Q) U5 Lfloat d2 = dist(X[l+1], X[r]);

! q4 J' q5 z. l6 ~- o$ Efloat d3 = dist(X[l], X[r]);

, |9 _/ u, \) [0 i- n// 寻找最近点对

3 |6 p) v, n# O. \7 _5 R4 Sif (d1 &lt;= d2 &amp;&amp; d1 &lt;= d3) {
& W. o# O/ S/ e- @# a- H1 @4 g
( k8 Z( m2 n/ ^$ ~( E* na = X[l];

b = X[l+1];
7 w. V- ~6 C" H
d = d1;

r e t u r n ; }

if (d2 &lt;= d3) {a = X[l+1];

& u0 I+ ?  H6 \1 `b = X[r];
) |1 k4 C6 x* v( U. X
4 q+ n: j2 x) j3 D. k9 F$ fd = d2;}

+ j4 f  _# o+ g1 C/ Q9 K2 Ielse {a = X[l];
1 ?5 o1 ^6 t8 k6 g
' R5 D5 c) I" I; Yb = X[r];

# A& G* H$ P1 W1 o4 h! q5 ^d = d3;}
/ e; o) V- F) @8 `) j- L
r e t u r n ; }
! `$ P1 O1 ^( n+ z$ p
' e3 N5 O* [. J" A% V/ _+ r/ /多于三个点，划分为两部分
; _* c0 S# I5 I, e/ ]3 k0 N
+ B: Q2 S: q' v( C( Yint m = (l+r)/2; // X[l:m] 在A中，余下的在B中

// 在Z[l:m] 和Z [ m + 1 : r ]中创建按y排序的表

- Y/ j2 B, l+ Q2 a% L( Kint f = l, // Z[l:m]的游标
! n3 R/ S4 b9 ?5 G  o
g = m+1; // Z[m+1:r]的游标
% G' h' t1 P; N/ `
for (int i = l; i &lt;= r; i++)
- ^6 u9 c8 S* c
if (Y.p &gt; m) Z[g++] = Y;
2 q6 ^( w: B$ J* i5 K6 {
else Z[f++] = Y;

5 ?! c$ L* B" }' t/ s9 s/ O// 对以上两个部分进行求解

c l o s e ( X , Z , Y, l , m , a , b , d ) ;

4 [% i9 X' B7 Kfloat dr;

) _% p# k1 y" `) g  e9 sPoint1 ar, br;
+ x8 o4 f9 y; S
7 n* @- S2 x+ H8 X( E0 ?# Sc l o s e ( X , Z , Y, m + 1 , r, a r, b r, d r ) ;

; @6 H* Q( |% c1 X- h// (a,b) 是两者中较近的点对

if (dr &lt; d) {a = ar;

b = br;
" u1 `% l5 \8 K  U$ C$ q
d = dr;}

M e r g e ( Z , Y,l,m,r);// 重构Y
: Q, Y  ?' R( H& n; _8 ?" t. s
4 d; c  l1 w" W) m8 p/ /距离小于d的点放入Z

) \- y( C& p( ~$ y! e' hint k = l; // Z的游标

6 I+ |2 x& }9 w" x) e, z9 Yfor (i = l; i &lt;= r; i++)

if (fabs(Y[m].x - Y.x) &lt; d) Z[k++] = Y;
; r/ o1 m0 h) d7 S
5 S; D* n9 r0 v( Q( c// 通过检查Z [ l : k - 1 ]中的所有点对，寻找较近的点对
: N. ~9 h0 t# A
for (i = l; i &lt; k; i++){

for (int j = i+1; j &lt; k &amp;&amp; Z[j].y - Z.y &lt; d;

j + + ) {
8 {# I4 \3 q  J  l4 p
% M  ?3 ]9 ~  T) o  B& Tfloat dp = dist(Z, Z[j]);
' ^' }+ ^* D9 n- [8 u9 L
if (dp &lt; d) {// 较近的点对

+ l  q: `- p7 y1 Gd = dp;
/ Z9 F$ {0 S/ c: C: |; c
% I) i0 e! w4 x8 l+ f# \a = X[Z.p];

b = X[Z[j].p];}

}

}

( H2 o' w: w5 p1 R+ T7 e3 d5 |}

# d2 j( ]! X9 @- Y% S函数c l o s e（见程序1 4 - 11）用来确定X[1:r] 中的最近点对。假定这些点按x 坐标排序。在Y [ 1 : r ]中对这些点按y 坐标排序。Z[ 1 : r ]用来存放中间结果。找到最近点对以后，将在a, b中返回最近点对，在d 中返回距离，数组Y被恢复为输入状态。函数并未修改数组X。

, z$ j& V7 d, }0 N- U. q首先考察“小问题”，即少于四个点的点集。因为分割过程不会产生少于两点的数组，因此只需要处理两点和三点的情形。对于这两种情形，可以尝试所有的可能性。当点数超过三个时，通过计算m = ( 1 + r ) / 2把点集分为两组A和B，X [ 1 : m ]属于A，X [ m + 1 : r ]属于B。通过从左至右扫描Y中的点以及确定哪些点属于A，哪些点属于B，可以创建分别与A组和B组对应的，按y 坐标排序的Z [ 1 : m ]和Z [ m + 1 : r ]。此时Y和Z的角色互相交换，依次执行两个递归调用来获取A和B中的最近点对。在两次递归调用返回后，必须保证Z不发生改变，但对Y则无此要求。不过，仅Y [ l : r ]可能会发生改变。通过合并操作（见程序1 4 - 5）可以以Z [ 1 : r ]重构Y [ 1 : r ]。
% }4 r. t. O" K1 a* ^6 L
+ v9 r. n3 w) `2 C' I为实现图1 4 - 1 6的策略，首先扫描Y [ 1 : r ]，并收集距分割线小于的点，将这些点存放在Z [ 1 : k - 1 ]中。可按如下两种方式来把RA中点p 与p 的比较区内的所有点进行配对：1) 与RB 中y 坐标≥p.y 的点配对；2) 与y 坐标≤p.y 的点配对。这可以通过将每个点Z [ i ]（1≤i &lt; k，不管该点是在RA
0 u) o6 H& K0 {9 m  v& c0 m
9 s2 o, ~4 j0 c4 E还是在RB中）与Z[j] 配对来实现，其中i＜j 且Z [ j ] . y - Z [ i ] . y＜ 。对每一个Z [ i ]，在2 × 区域内所检查的点如图1 4 - 1 7所示。由于在每个2 × 子区域内的点至少相距。因此每一个子区域中的点数不会超过四个，所以与Z [ i ]配对的点Z [ j ]最多有七个。
6 C" m- o8 B; |4 Y
* h8 G( H, {/ q5 x2. 复杂性分析

* L3 X" s. E: _, Q5 q: R令t (n) 代表处理n 个点时，函数close 所需要的时间。当n＜4时，t (n) 等于某个常数d。当n≥4时，需花费(n) 时间来完成以下工作：将点集划分为两个部分，两次递归调用后重构Y，淘汰距分割线很远的点，寻找更好的第三类点对。两次递归调用需分别耗时t (「n /2ù」和t (?n /2?).
* Q$ m3 @- \9 Z% R7 \( ]6 h5 n$ t5 m
这个递归式与归并排序的递归式完全一样，其结果为t (n) = (nl o gn)。另外，函数c l o s e s t还需耗时(nl o gn)来完成如下额外工作：对X进行排序，创建Y和Z，对Y进行排序。因此分而治之最近点对求解算法的时间复杂性为(nl o gn)。</P>