距离最近的点对

韩冰

<>给定n 个点（xi，yi）（1≤i≤n），要求找出其中距离最近的两个点。
; N1 X( z; G" Q7 G  _4 C
例14-7 假设在一片金属上钻n 个大小一样的洞，如果洞太近，金属可能会断。若知道任意两个洞的最小距离，可估计金属断裂的概率。这种最小距离问题实际上也就是距离最近的点对问题。
- t. `& u. Y0 A# u: ^( U4 F
通过检查所有的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+ ^1 B  m8 D. u8 y
& Q1 @- T+ x. K0 ~9 b  Yif (n较小) {用直接法寻找最近点对
* ^" c, L( i3 O9 X7 ]  _
R e t u r n ; }

0 V$ ^5 _* L: q! P8 I& i- ^// n较大
7 P5 R3 `. Q& |; ^& h) l) @, j
将点集分成大致相等的两个部分A和B

2 w/ R! |4 v; v) M& @. |- l确定A和B中的最近点对
3 a4 j5 d7 v9 Z) ~! g% J
1 R' u; L) Q: e% u( H确定一点在A中、另一点在B中的最近点对

+ p9 Z3 ~3 o4 ^' Z& F5 v& ]从上面得到的三对点中，找出距离最小的一对点
! Z+ p. T/ J; N) I2 E6 D3 m
" M+ r' W" P2 Z* {, V, J& l5 X图14-13 寻找最近的点对


: x5 D) |* C* F3 M$ ?1 p为了确定第三种情况下的最近点对，需要采用一种不同的方法。这种方法取决于点集是如何被划分成A、B的。一个合理的划分方法是从xi（中间值）处划一条垂线，线左边的点属于A，线右边的点属于B。位于垂线上的点可在A和B之间分配，以便满足A、B的大小。

例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。
( F& j( D  r& _- n
; O$ r! L3 T, @. W0 g0 S. T设是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）。
! K( V+ o6 _* b: Q- w  G' H. [
3 {! y$ o) p9 D4 ]2 W例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) 是最近的点对。

为了确定一个距离更小的第三类点，RA 中的每个点最多只需和RB 中的6个点比较，如图1 4 - 1 6所示。

: H" L/ H  H- k1. 选择数据结构
3 A9 Q! @( `' d6 W* {
3 ^8 M9 ?3 ?/ u为了实现图1 4 - 1 3的分而治之算法，需要确定什么是“小问题”以及如何表示点。由于集合中少于两点时不存在最近点对，因此必须保证分解过程不会产生少于两点的点集。如果将少于四点的点集做为“小问题”，就可以避免产生少于两点的点集。

/ M! _9 \% t2 q每个点可有三个参数：标号， x 坐标，y 坐标。假设标号为整数，每个点可用P o i n t l类（见程序1 4 - 8）来表示。为了便于按x 坐标对各个点排序，可重载操作符＜=。归并排序程序如1 4 -3所示。

8 j; d" U; j, y$ c程序14-8 点类

class Point1 {
! {0 g" ?; x, P7 W5 ?2 z5 [. d- J# A
. m; `/ k- e& {2 N% |& efriend float dist(const Point1&amp;, const Point1&amp;);
/ a. S6 T! u6 |
friend void close(Point1 *, Point2 *, Point2 *, int, int, Point1&amp;, Point1&amp;, float&amp;);

" m% O3 q1 Z, F/ r8 U+ s2 @friend bool closest(Point1 *, int, Point1&amp;, Point1&amp;,float&amp;);
! K& i; J1 N  I
friend void main();

  v5 q6 A: s3 A/ q: j* jp u b l i c :
: H6 U- b! A+ v" z7 A
int operator&lt;=(Point1 a) const

{return (x &lt;= a.x);}
# G1 i% _# V+ J7 r8 e( n
p r i v a t e :
- S, j% z) k; ^5 J5 F% u
6 p9 |& i9 C  C) r. Wint ID; // 点的编号
; U) z6 f' f+ d4 R; ~$ i* S
float x, y; // 点坐标
, o. `0 P. n, k. w  l8 c
* Y" Z# V  K, b' o} ;
! S3 G2 J6 u: E/ q5 }
class Point2 {
& x- H5 A/ n1 M' ]( D9 J) E  {
' q3 [0 S* r2 L3 I" c5 gfriend float dist(const Point2&amp;, const Point2&amp;);
/ W7 @- m% b' Y) _: ^8 i- _9 ^1 j1 b
% e; f% H' L" X7 g' |& Ofriend void close(Point1 *, Point2 *, Point2 *, int, int, Point1&amp;, Point1&amp;, float&amp;);

friend bool closest(Point1 *, int, Point1&amp;, Point1&amp;, float&amp;);
" }$ r! `! N' k7 M# C) L4 X& `
  A$ O  H. A; L) B' \friend void main();

p u b l i c :

int operator&lt;=(Point2 a) const
2 J/ T0 T& V! S( Z' x& D9 m
{return (y &lt;= a.y);}

p r i v a t e :
) w/ s3 Q0 y, {" z: A
int p; // 数组X中相同点的索引
$ q) l1 z) Z2 ^4 V8 h
5 t# ^# a0 ]4 e9 Cfloat x, y; // 点坐标

6 I1 Y5 O3 J! l7 f; n, t} ;

! d* \% ]. m- \/ w* I8 d所输入的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中的对应点。

确定了必要的数据结构之后，再来看看所要产生的代码。首先定义一个模板函数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类的友元。
1 ]4 J1 h, ^3 _( {  E
. Q, q; }& U5 e0 h! z3 f3 I  f6 F程序14-9 计算两点距离
2 q5 M9 R5 u* J0 s' i# a8 b0 d' @
7 A7 K" m. Y: ztemplate<CLASS T>

inline float dist(const T&amp; u, const T&amp; v)
: e) Z! L5 `  c0 n+ }  R4 B
. g) p' P' }* N' A$ D: B  [{ / /计算点u 和v之间的距离
7 n  o' O) s6 K# P, n3 a
float dx = u.x-v. x ;
* m5 ]7 A2 k3 \9 [' A. r6 R+ z
float dy = u.y-v. y ;

return sqrt(dx * dx + dy * dy);

}

3 _- K$ P- h; M1 b/ |& g/ L, 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 )，该函数实际求解最近点对。
8 h' K8 d; P1 D& C4 d: g
% o2 p* a5 l: R/ W% }程序14-10 预处理及调用c l o s e
8 k* E2 D1 m( a7 n% X, @
bool closest(Point1 X[], int n, Point1&amp; a, Point1&amp; b, float&amp; d)
* ]1 e0 H8 _3 t- E& ~+ x& U
6 O; S. @& Y4 n# |3 P, A0 g( j' q{// 在n &gt;= 2 个点中寻找最近点对
/ s" M; u- \! v9 Y
/ Y$ h6 y; j  z// 如果少于2个点，则返回f a l s e

" Y% `7 i' Q' R// 否则，在a 和b中返回距离最近的两个点

if (n &lt; 2) return false;

) d! N4 M% E+ ~' z/ D2 s7 m1 b* x// 按x坐标排序

" h( |9 C3 n; t( R3 {4 qM e r g e S o r t ( X , n ) ;
; c; R0 P4 `2 R8 N0 }
/ G0 A% S0 X( `// 创建一个按y坐标排序的点数组
8 @, `7 [+ t0 U+ x8 O" B
Point2 *Y = new Point2 [n];

for (int i = 0; i &lt; n; i++) {
; ^5 A& h" `. U/ l( Z7 P
) E; X9 a, f6 e, L; D  D6 g& P* W// 将点i 从X 复制到Y

Y.p = i;
" |8 T( i5 Y. k2 |& @7 k  A
Y.x = X.x;

Y.y = X.y;
) V0 w# z9 b$ }, w( f" g
}
. Y3 U5 d4 w" v8 r' X) ?! S
M e r g e S o r t ( Y,n); // 按y坐标排序
5 G; a0 W: Q3 t# S6 k! k! x; O* T* i
// 创建临时数组

0 X5 S- Z0 H' l, hPoint2 *Z = new Point2 [n];
/ o" q$ K3 w! i1 g
// 寻找最近点对
3 a" j" e2 h# K7 P4 p5 O
c l o s e ( X , Y, Z , 0 , n - 1 , a , b , d ) ;
) y7 f( m: X" i' V% \; k4 \
// 删除数组并返回
$ w! g) `2 b) {. `" I
delete [] Y;
" D6 ], d0 j& I7 }; ]4 f
delete [] Z;

return true;

" m* ?3 S7 g, z4 ?# s1 t+ o* h- ]3 o}

4 ~" [; v) N2 t' b% Q# k程序1 4 - 11 计算最近点对
. i! c2 v/ N' K2 n% ^+ C: C
& p9 w% d5 l4 e% ~void close(Point1 X[], Point2 Y[], Point2 Z[], int l, int r, Point1&amp; a, Point1&amp; b, float&amp; d)

{//X[l:r] 按x坐标排序
- T; ^5 S$ |0 `6 |& y  @
) X9 |- i$ ]! [( {  ^//Y[l:r] 按y坐标排序
% x8 ^3 _& ]4 _
$ I0 w5 o9 v6 i, x/ X8 tif (r-l == 1) {// 两个点
' t: A  k' n- s# J) F+ {- S3 `
+ r% c' }- N- e5 T# M$ P3 za = X[l];
( z2 Y- o) z8 q
b = X[r];
2 b* @9 N# M( ]- {! k. y3 T
d = dist(X[l], X[r]);

9 W1 G, Q6 b  f* c7 `! Tr e t u r n ; }

- M: b. E. |4 K& M/ B/ rif (r-l == 2) {// 三个点

7 B3 b3 W5 e) f' }9 J. t* V( V// 计算所有点对之间的距离
$ C% y, q  v8 t4 T: e
1 s+ C$ ?% S, j+ w) Zfloat d1 = dist(X[l], X[l+1]);
9 c6 s, a: q( b! {2 l! D0 G$ T
float d2 = dist(X[l+1], X[r]);
( k% f) d. `# W& k* v$ i
float d3 = dist(X[l], X[r]);
( U2 I4 a& y/ Q# }" L
: f& ?9 Y9 T  n8 q( O. o* j// 寻找最近点对

if (d1 &lt;= d2 &amp;&amp; d1 &lt;= d3) {
$ A6 Q  p" O% ^7 r
; k6 y) L+ g, z: na = X[l];
4 \1 W" G- f- C# W! ?+ C+ N2 A1 Y
b = X[l+1];
( p, Z1 Y8 ^( o
1 i( q1 Q# o) H" ]4 d% i8 O, fd = d1;

r e t u r n ; }
) J3 W& r. L/ Y9 I* u5 }
7 j, U( P/ z- n1 ~% E$ F' h7 Rif (d2 &lt;= d3) {a = X[l+1];

b = X[r];

7 X2 K  Q  j, w+ E* w' D  sd = d2;}
% p% _' v: B, `1 C8 o/ ~$ Z9 m
else {a = X[l];
  y& g. N$ b8 g. h7 v0 U8 b# y9 d
b = X[r];

d = d3;}
, w' B$ f3 K' y! v" I8 r
r e t u r n ; }

/ /多于三个点，划分为两部分
7 g; |0 B6 ]; C$ D5 r
! k3 |4 s( g- l7 q# j2 q1 n; Hint m = (l+r)/2; // X[l:m] 在A中，余下的在B中
) t* f) u+ t. g( Z
  D: j1 Z8 g5 ?4 R: X. G2 i/ \// 在Z[l:m] 和Z [ m + 1 : r ]中创建按y排序的表
. ]! }+ E3 @( x  h+ X3 E2 E
int f = l, // Z[l:m]的游标

g = m+1; // Z[m+1:r]的游标

& S6 F, i8 K$ S. Ifor (int i = l; i &lt;= r; i++)
3 S2 O6 k3 [8 i9 [! T& h
if (Y.p &gt; m) Z[g++] = Y;

* u- Z* D9 c) V( q9 C, j/ r# felse Z[f++] = Y;

, D) g7 e1 d; b4 U' t6 P/ o// 对以上两个部分进行求解
( d# m* E8 x& H
c l o s e ( X , Z , Y, l , m , a , b , d ) ;

+ d/ T! ^! F- Hfloat dr;
4 d/ z( [! S7 G0 s, a, N$ @% `# u
$ q6 Q- Z% _. i. X7 s0 Y8 ~& |Point1 ar, br;

4 ~. P4 @5 ]( C* \c l o s e ( X , Z , Y, m + 1 , r, a r, b r, d r ) ;

1 M7 K! v4 c' R' S2 u% a6 T1 g  \// (a,b) 是两者中较近的点对

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

b = br;

d = dr;}

+ A% C0 J) _3 }9 L  [M e r g e ( Z , Y,l,m,r);// 重构Y
: Q5 i! N- P$ |" K# K3 Y8 v# \
/ /距离小于d的点放入Z
& n" S1 D* h7 c5 l; H
int k = l; // Z的游标
5 i' N! R: ]( @7 L4 g$ K& h
! q6 t1 B+ P' y7 ?for (i = l; i &lt;= r; i++)
9 Q6 Q9 b4 w" p8 w( X
if (fabs(Y[m].x - Y.x) &lt; d) Z[k++] = Y;

// 通过检查Z [ l : k - 1 ]中的所有点对，寻找较近的点对
5 E5 J* p1 h! X3 R
4 t# t& c0 h; ]8 M$ Qfor (i = l; i &lt; k; i++){
  B8 E& j% M+ r
; N6 m; s' E5 @3 n  C5 m" ~- G, Bfor (int j = i+1; j &lt; k &amp;&amp; Z[j].y - Z.y &lt; d;
  g; S/ K0 v% i' h9 U
j + + ) {
% _3 [1 x' u- v/ m& g" v
float dp = dist(Z, Z[j]);
. o2 V! E* u& ?. B
  r2 h/ N0 Z& b4 Pif (dp &lt; d) {// 较近的点对

d = dp;

a = X[Z.p];
4 n  w4 A# B7 a* \4 ~/ o1 w
b = X[Z[j].p];}
/ V  R$ R* L" B- h$ J
! g! W( M4 V/ p: Y0 @3 s}

}

$ K" |# O- S( m! T5 J* A+ 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。
1 B3 u4 F: Z& X) T) D8 P' G
首先考察“小问题”，即少于四个点的点集。因为分割过程不会产生少于两点的数组，因此只需要处理两点和三点的情形。对于这两种情形，可以尝试所有的可能性。当点数超过三个时，通过计算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 ]。
7 r  d% t/ T5 q- d. v) E
为实现图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
) W5 n* A' [8 ^
2 E& B( L% d' C. P+ E8 j还是在RB中）与Z[j] 配对来实现，其中i＜j 且Z [ j ] . y - Z [ i ] . y＜ 。对每一个Z [ i ]，在2 × 区域内所检查的点如图1 4 - 1 7所示。由于在每个2 × 子区域内的点至少相距。因此每一个子区域中的点数不会超过四个，所以与Z [ i ]配对的点Z [ j ]最多有七个。
( F* P( w7 D/ k& F8 U; j) |% o& M
. |5 O% j/ m( \" C+ Y+ V) t) e2. 复杂性分析

3 u4 T- q0 z! c( e- y) r+ [令t (n) 代表处理n 个点时，函数close 所需要的时间。当n＜4时，t (n) 等于某个常数d。当n≥4时，需花费(n) 时间来完成以下工作：将点集划分为两个部分，两次递归调用后重构Y，淘汰距分割线很远的点，寻找更好的第三类点对。两次递归调用需分别耗时t (「n /2ù」和t (?n /2?).

. J. ?% J8 K# v3 _+ T4 ~这个递归式与归并排序的递归式完全一样，其结果为t (n) = (nl o gn)。另外，函数c l o s e s t还需耗时(nl o gn)来完成如下额外工作：对X进行排序，创建Y和Z，对Y进行排序。因此分而治之最近点对求解算法的时间复杂性为(nl o gn)。</P>