距离最近的点对

韩冰

<>给定n 个点（xi，yi）（1≤i≤n），要求找出其中距离最近的两个点。
1 @* X# q+ n. C+ O$ o+ x  s3 f" o
例14-7 假设在一片金属上钻n 个大小一样的洞，如果洞太近，金属可能会断。若知道任意两个洞的最小距离，可估计金属断裂的概率。这种最小距离问题实际上也就是距离最近的点对问题。

5 Q8 b$ [8 h) M$ Z, T通过检查所有的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进行递归求解。


if (n较小) {用直接法寻找最近点对

, v2 C* Y: U- R3 V/ \  }R e t u r n ; }
3 l5 M3 D+ Z' u/ @
: i* o4 ~5 }8 L// n较大
! G1 J5 @  S% q  }
将点集分成大致相等的两个部分A和B

确定A和B中的最近点对
0 O7 Y2 X( m, G4 d0 a
确定一点在A中、另一点在B中的最近点对
2 N, X: i! {- E/ F( K: }
. A% p& v. ]1 Y, T# }2 ~从上面得到的三对点中，找出距离最小的一对点

图14-13 寻找最近的点对


为了确定第三种情况下的最近点对，需要采用一种不同的方法。这种方法取决于点集是如何被划分成A、B的。一个合理的划分方法是从xi（中间值）处划一条垂线，线左边的点属于A，线右边的点属于B。位于垂线上的点可在A和B之间分配，以便满足A、B的大小。
( a* v0 D2 c) i- |; h
, h2 ]: T5 b6 Y5 o5 _4 f+ M例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。
, P1 Z/ M; o: }0 L" ]# C0 t* o$ z: n
设是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）。

例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) 是最近的点对。
( y0 u% `1 k' n5 {3 T
; F2 q& Y1 H" z2 n0 d为了确定一个距离更小的第三类点，RA 中的每个点最多只需和RB 中的6个点比较，如图1 4 - 1 6所示。
2 w+ r  ]6 [! W5 U& M
" ^" }: F7 K& C% D2 ^& q! c1. 选择数据结构
  J$ Z/ C& a# f
为了实现图1 4 - 1 3的分而治之算法，需要确定什么是“小问题”以及如何表示点。由于集合中少于两点时不存在最近点对，因此必须保证分解过程不会产生少于两点的点集。如果将少于四点的点集做为“小问题”，就可以避免产生少于两点的点集。
+ d2 X- s9 z* d8 }5 d
& v& G7 t  Y9 Y+ w$ l; a0 u每个点可有三个参数：标号， x 坐标，y 坐标。假设标号为整数，每个点可用P o i n t l类（见程序1 4 - 8）来表示。为了便于按x 坐标对各个点排序，可重载操作符＜=。归并排序程序如1 4 -3所示。
% W+ e, a1 n2 D# v
1 R/ G* R8 {+ ?* w5 g程序14-8 点类

+ c) H2 z- C4 a) P! Q$ Yclass Point1 {

friend float dist(const Point1&amp;, const Point1&amp;);
$ Q1 j( \. t$ Y' `/ u% T, K
2 S2 P8 ]& p$ }7 Y  sfriend void close(Point1 *, Point2 *, Point2 *, int, int, Point1&amp;, Point1&amp;, float&amp;);
- D7 ~, P( L& Z4 \% m& Q
6 b/ ~9 a3 j: ?4 ^- bfriend bool closest(Point1 *, int, Point1&amp;, Point1&amp;,float&amp;);

0 D. H. O3 a& Ufriend void main();
/ P/ a5 D3 Y# w" Z
, k- e, K7 i9 ^0 Bp u b l i c :
+ ?/ |$ I/ w8 o- t
4 F% ]" b; @& q( T. k& Aint operator&lt;=(Point1 a) const

{return (x &lt;= a.x);}
, y4 q( ]& ], \9 ]' h* `
9 C  ~* @* _. T7 I7 Ep r i v a t e :
" X( Z2 m) M' }
' {; S: e$ F* H$ B2 H  [* {; s  sint ID; // 点的编号

- P* U, G% {$ H. m0 t3 h, `& y3 ~! efloat x, y; // 点坐标

. t% l5 H6 D, R/ G% x/ ~, }2 q} ;
- N8 e7 r  K- c; G: |' U
class Point2 {
+ S  x- J+ [" z
friend float dist(const Point2&amp;, const Point2&amp;);

8 i  h* H" M  b# Vfriend void close(Point1 *, Point2 *, Point2 *, int, int, Point1&amp;, Point1&amp;, float&amp;);
4 i+ M* C/ [# C: r  ~
  R; z( l$ G9 W* ~9 l7 J3 mfriend bool closest(Point1 *, int, Point1&amp;, Point1&amp;, float&amp;);
) i2 _9 k# X  U
2 V( ~1 r* Z) bfriend void main();
7 f7 e- t/ M  I( m% g2 h, R
+ X; ^# e8 R: B$ i( J+ Qp u b l i c :
. y5 K4 m6 U9 x+ Y
int operator&lt;=(Point2 a) const

# _' @/ y% r% B{return (y &lt;= a.y);}

p r i v a t e :
* Y' ^+ A) u1 M  o) F8 `
. Q; L, D" H2 t/ w# F( a. _# l' vint p; // 数组X中相同点的索引

5 ~0 ^, i8 N1 b9 nfloat x, y; // 点坐标

1 M7 H: a& ^! M1 |; Q} ;

5 d9 a3 K9 O5 {9 V1 b1 x所输入的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中的对应点。

( E9 S- B- [- V: Q7 a. c确定了必要的数据结构之后，再来看看所要产生的代码。首先定义一个模板函数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 计算两点距离
  b9 |! d1 a) B7 |  e9 h+ h7 S, l
template<CLASS T>
" j) h5 L; v/ R" J" I
inline float dist(const T&amp; u, const T&amp; v)

$ P8 _8 O$ c* _! V  B{ / /计算点u 和v之间的距离

float dx = u.x-v. x ;

float dy = u.y-v. y ;

2 M9 s8 i$ a& I' r! Q/ N4 i$ o+ Preturn sqrt(dx * dx + dy * dy);
+ U* t  Z* b6 `
}

+ I7 D, \3 L7 ^如果点的数目少于两个，则函数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 )，该函数实际求解最近点对。

程序14-10 预处理及调用c l o s e

bool closest(Point1 X[], int n, Point1&amp; a, Point1&amp; b, float&amp; d)
( H, d7 T& v2 p( ~) @
' k, Y! @- n$ T9 T+ F7 Y- T' k: s{// 在n &gt;= 2 个点中寻找最近点对
! z. s+ [+ V8 J0 Y7 a2 y9 }2 C
// 如果少于2个点，则返回f a l s e
+ x: k+ t+ j7 b
3 |  j! L$ h! j3 b// 否则，在a 和b中返回距离最近的两个点
, E+ {+ f3 B* Z. D& F: M; O5 m
if (n &lt; 2) return false;

// 按x坐标排序
8 ^; _7 s: k2 F7 w! W. S2 a
M e r g e S o r t ( X , n ) ;

// 创建一个按y坐标排序的点数组

$ w8 h3 K' h5 m8 ~Point2 *Y = new Point2 [n];

4 O2 P9 \4 m4 Q! Dfor (int i = 0; i &lt; n; i++) {
- a" i$ ?$ R! u! D8 h
! ~5 k* Z, z2 R7 ]: L7 A7 R3 {// 将点i 从X 复制到Y

Y.p = i;
5 c1 N$ x- G; A9 r7 {
2 @; a) M( r& g8 K& L/ |- `Y.x = X.x;
' ]9 t; R/ _9 G% Q
7 \! I8 M7 M1 L! @. |  W. oY.y = X.y;
3 @1 S' S; E& V2 h+ ]
* b( R( @6 M5 ^2 a}
: L2 O( T4 d0 q0 D$ I
M e r g e S o r t ( Y,n); // 按y坐标排序

// 创建临时数组

Point2 *Z = new Point2 [n];

// 寻找最近点对

; L0 w8 `6 y, ]& J: ^c l o s e ( X , Y, Z , 0 , n - 1 , a , b , d ) ;
  D( }& j7 w- U
// 删除数组并返回

: Y5 i' b  }! @$ Z$ odelete [] Y;

delete [] Z;
9 u0 H  X0 A: I( C( H
return true;

}
% p$ y- @+ U# M& H! J* q% C
程序1 4 - 11 计算最近点对

' N) @* @. G- S( T% C% n5 I7 ~void close(Point1 X[], Point2 Y[], Point2 Z[], int l, int r, Point1&amp; a, Point1&amp; b, float&amp; d)
+ Q+ H# x% G- R; ?; l. J
{//X[l:r] 按x坐标排序
: L7 o( X" _- O
//Y[l:r] 按y坐标排序

. R. o# M' w4 M7 G, L' |if (r-l == 1) {// 两个点
1 E! b; D) Z. N+ h. M. k5 W& Z
3 E, k0 A) X, U: r4 [; g6 Ga = X[l];

b = X[r];
- h; ?/ |; G- c- @! [0 y. o
d = dist(X[l], X[r]);
$ r  z0 \, ~9 e( E
r e t u r n ; }

% {- u% J* `- z1 ~4 x8 Hif (r-l == 2) {// 三个点
6 ~0 Y- s6 v( {( F- Q  n3 P
6 w8 N! j( @+ I( B' N( q9 }// 计算所有点对之间的距离

float d1 = dist(X[l], X[l+1]);
+ W0 ^& `& `5 {( R2 K: p, n4 [
2 o. q! K- g, gfloat d2 = dist(X[l+1], X[r]);

# M$ f( ~$ @& Yfloat d3 = dist(X[l], X[r]);
7 a0 v2 @1 ]$ O- X7 \  E
4 l: B% Z3 I; f: @" B: C// 寻找最近点对
) b5 f: e6 x) b/ D8 E
8 Q. E0 _* c- m/ hif (d1 &lt;= d2 &amp;&amp; d1 &lt;= d3) {

a = X[l];

7 B9 k' I8 R6 `' l1 I6 f6 A' b" ob = X[l+1];

d = d1;

r e t u r n ; }
$ r2 v9 E# E3 [7 Z  m: q3 }
  b; F! ?' E& @  R1 e1 Mif (d2 &lt;= d3) {a = X[l+1];

" _  |3 k% U" i5 X" J3 Yb = X[r];

d = d2;}

9 U3 Z: E% M7 A* H3 }% S% oelse {a = X[l];
9 h0 v% P* C( g/ f/ ^* r2 S$ Y/ d- i
3 T4 Z3 ?7 _3 bb = X[r];
5 Q0 F# G8 L9 t. x# \; `
# E  Z0 z' @- Rd = d3;}

8 r4 I5 R* i. L4 ^r e t u r n ; }

: o9 n) [; K+ c  C/ /多于三个点，划分为两部分

int m = (l+r)/2; // X[l:m] 在A中，余下的在B中

0 _: M& V- i. h3 g1 d3 \' \$ ?// 在Z[l:m] 和Z [ m + 1 : r ]中创建按y排序的表
* I0 E$ ~6 S. d( f& }( E. E
int f = l, // Z[l:m]的游标
* b3 K# y7 _$ ~$ B/ g! ~; ^
g = m+1; // Z[m+1:r]的游标

for (int i = l; i &lt;= r; i++)

if (Y.p &gt; m) Z[g++] = Y;

else Z[f++] = Y;
8 |  |5 B$ B/ V; {3 l
9 y1 |* o: {2 Y( P/ r8 q// 对以上两个部分进行求解
) w# N4 F& l! {' \/ M& j% z1 j
- B" Q; E8 O3 d! ac l o s e ( X , Z , Y, l , m , a , b , d ) ;

float dr;
# Z4 O0 `' _+ Y" y% m" J; x9 O* W
3 g, B+ H; {! bPoint1 ar, br;

! U0 ^$ x2 k, a7 |. \6 @c l o s e ( X , Z , Y, m + 1 , r, a r, b r, d r ) ;

// (a,b) 是两者中较近的点对
$ L; Q8 _& g! I- G% r
6 u; s1 y8 [4 b- d3 ]if (dr &lt; d) {a = ar;
# ~4 D7 B8 U' P4 L/ ~
& @% f$ m5 v, m7 z3 P/ Wb = br;
. C& E% d9 v3 F5 e8 ~" M: y& @
2 _8 c6 n2 m+ ?* ]9 l, C, Hd = dr;}
, j/ {2 M. r$ E6 u% |  b
: S7 a6 u: ?# d' E- PM e r g e ( Z , Y,l,m,r);// 重构Y
- z% v  G5 T3 P1 _
4 Q0 P3 w1 m$ h4 z$ m/ /距离小于d的点放入Z
& k4 O. O3 ~; _7 S) m
int k = l; // Z的游标
6 q+ q/ L7 R/ P& v1 d8 Y9 e
; Q9 B- a: p* o+ efor (i = l; i &lt;= r; i++)
0 Y/ O4 j1 s$ s
if (fabs(Y[m].x - Y.x) &lt; d) Z[k++] = Y;
: p8 |4 |% M2 f8 h: R- k
// 通过检查Z [ l : k - 1 ]中的所有点对，寻找较近的点对

for (i = l; i &lt; k; i++){
: |+ r. E* p5 V; I7 D/ o0 A) p
for (int j = i+1; j &lt; k &amp;&amp; Z[j].y - Z.y &lt; d;
3 ]/ T3 g4 g' X9 o2 x6 J# ]5 E
, O: h. B. y9 a' \$ Bj + + ) {

float dp = dist(Z, Z[j]);

8 ~- q0 _8 ?8 Q9 ]9 V3 }if (dp &lt; d) {// 较近的点对

d = dp;
" H2 F3 |# d( Z
9 Q) y/ U" X  b4 c+ Y7 Xa = X[Z.p];
, X: s% ~& v' J  s+ K) {. _
3 w; w0 ]. s4 Ub = X[Z[j].p];}

}
+ [/ ]  n! e) y* d; j, j9 r
}
4 k; H& S' [) L- e3 ?
8 m8 @% J: d! k0 {" P# t. t}

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

首先考察“小问题”，即少于四个点的点集。因为分割过程不会产生少于两点的数组，因此只需要处理两点和三点的情形。对于这两种情形，可以尝试所有的可能性。当点数超过三个时，通过计算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 ]。
  R7 X  M: X% M' Y$ K0 k6 B9 x
为实现图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

还是在RB中）与Z[j] 配对来实现，其中i＜j 且Z [ j ] . y - Z [ i ] . y＜ 。对每一个Z [ i ]，在2 × 区域内所检查的点如图1 4 - 1 7所示。由于在每个2 × 子区域内的点至少相距。因此每一个子区域中的点数不会超过四个，所以与Z [ i ]配对的点Z [ j ]最多有七个。

2. 复杂性分析
3 ~( W. w% m5 [6 Z/ k  W
令t (n) 代表处理n 个点时，函数close 所需要的时间。当n＜4时，t (n) 等于某个常数d。当n≥4时，需花费(n) 时间来完成以下工作：将点集划分为两个部分，两次递归调用后重构Y，淘汰距分割线很远的点，寻找更好的第三类点对。两次递归调用需分别耗时t (「n /2ù」和t (?n /2?).
) {7 e; a% \: h( g/ |8 A) Y0 E& L
& m6 k8 f2 @4 A( }% _/ z这个递归式与归并排序的递归式完全一样，其结果为t (n) = (nl o gn)。另外，函数c l o s e s t还需耗时(nl o gn)来完成如下额外工作：对X进行排序，创建Y和Z，对Y进行排序。因此分而治之最近点对求解算法的时间复杂性为(nl o gn)。</P>