距离最近的点对

韩冰

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

例14-7 假设在一片金属上钻n 个大小一样的洞，如果洞太近，金属可能会断。若知道任意两个洞的最小距离，可估计金属断裂的概率。这种最小距离问题实际上也就是距离最近的点对问题。
& H$ a+ P, n; {
通过检查所有的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进行递归求解。
/ r0 S* c" e" Z0 H8 Q+ X, r# N

% S/ z" z% X7 x* `. Wif (n较小) {用直接法寻找最近点对

( J& H! \! N% _( G  A! mR e t u r n ; }

) M& U: z2 @. t; B+ O7 T3 B) L7 b// n较大

; T" k4 z# ^; S, q9 _! m+ v将点集分成大致相等的两个部分A和B

) q9 P. s3 v3 k' \确定A和B中的最近点对

4 i' s( V, D* M& z0 F2 C. g确定一点在A中、另一点在B中的最近点对

从上面得到的三对点中，找出距离最小的一对点

图14-13 寻找最近的点对
1 ?/ |) A/ [4 Y. N8 ^1 \, C

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

' I6 v$ i' p( j例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。
4 Z7 V3 |. \5 ~- X$ t% [! R
% I3 a- G' \- K  T- M7 F  z, t1 s0 k设是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）。
6 a1 s8 _+ _2 C# x1 B, p- i, K! ]
8 G+ D2 A+ W+ ~4 Y; z% R例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) 是最近的点对。
  a6 H2 p1 y- V) [# U8 K
8 o' @- u8 ?6 ~4 s: w* N为了确定一个距离更小的第三类点，RA 中的每个点最多只需和RB 中的6个点比较，如图1 4 - 1 6所示。
9 M7 ?! S# c5 T0 N
1. 选择数据结构
7 `8 ^! ]2 U1 q6 m: ]  G7 A1 N
& g: H9 t  a0 z; u- z4 u6 j为了实现图1 4 - 1 3的分而治之算法，需要确定什么是“小问题”以及如何表示点。由于集合中少于两点时不存在最近点对，因此必须保证分解过程不会产生少于两点的点集。如果将少于四点的点集做为“小问题”，就可以避免产生少于两点的点集。
) ]7 H( Q; G9 {2 p  |
3 \/ S( T3 G! G$ [每个点可有三个参数：标号， x 坐标，y 坐标。假设标号为整数，每个点可用P o i n t l类（见程序1 4 - 8）来表示。为了便于按x 坐标对各个点排序，可重载操作符＜=。归并排序程序如1 4 -3所示。

2 S) ^; H0 I( F: `, V: F3 j( u& `程序14-8 点类

class Point1 {
: j+ q" m8 _2 W# Y1 m. D
) c5 ]4 O. ?* j7 m' Y0 T8 D; p- Lfriend float dist(const Point1&amp;, const Point1&amp;);

9 Y. }7 U1 f2 A$ L, P/ ofriend void close(Point1 *, Point2 *, Point2 *, int, int, Point1&amp;, Point1&amp;, float&amp;);

0 d9 u* h) ?* O% F8 ]friend bool closest(Point1 *, int, Point1&amp;, Point1&amp;,float&amp;);

% E4 Y6 @% W) f2 [( {1 L1 xfriend void main();

5 z( x. a. P+ @. s) tp u b l i c :

# a- z* j/ z3 U+ {8 Y( Y2 @: t- ~. Z  qint operator&lt;=(Point1 a) const
9 ^* m, q9 k0 ~; x
- }# l9 L) M+ {, R4 e% e{return (x &lt;= a.x);}
' e2 g) L& @& I3 L' v" z- ?
8 @# v; V3 n7 W2 j" Dp r i v a t e :
0 G* J' |% D6 ~$ @3 a3 E# U
int ID; // 点的编号
' U& U/ F) D; O( S
float x, y; // 点坐标
$ ]$ Q% B% O. N0 M1 B
} ;

class Point2 {

& S2 e  f' X' c! {  Z: Xfriend float dist(const Point2&amp;, const Point2&amp;);

friend void close(Point1 *, Point2 *, Point2 *, int, int, Point1&amp;, Point1&amp;, float&amp;);
2 x* B" @& y: }: i/ @9 |- P- D5 P7 \
& ]7 e, c4 x7 f. }" q) ?friend bool closest(Point1 *, int, Point1&amp;, Point1&amp;, float&amp;);

" N; p8 ]' b* }9 ~8 O- s. ~friend void main();

p u b l i c :
+ w; Q/ Q0 m# B, ], o. E3 v9 G+ M
int operator&lt;=(Point2 a) const
# K' m7 R5 E4 b6 b; Q$ R" j
{return (y &lt;= a.y);}

" g. ], R' }* q9 `p r i v a t e :

int p; // 数组X中相同点的索引

* n$ p# d  A$ q2 ?5 mfloat x, y; // 点坐标
$ B3 n- ~4 Y" j% t; J2 k
} ;

所输入的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中的对应点。
4 ?. s2 s( Q& d: W
确定了必要的数据结构之后，再来看看所要产生的代码。首先定义一个模板函数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类的友元。
+ Z- L1 k/ G# }' O, ~" J2 Z. f9 r
程序14-9 计算两点距离

template<CLASS T>

inline float dist(const T&amp; u, const T&amp; v)
; J) X7 k- Q# ^4 ?6 q7 T6 a( u% Q+ z
5 S" N7 d% ?) c3 E& h$ O$ `$ U  ?+ |{ / /计算点u 和v之间的距离

float dx = u.x-v. x ;

float dy = u.y-v. y ;

return sqrt(dx * dx + dy * dy);
' B( @' A* h( N) i& O6 M$ T
}
0 I  {- b2 A/ n& u7 q- 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 )，该函数实际求解最近点对。

程序14-10 预处理及调用c l o s e
3 `5 c( z( a0 G2 j5 F! U
bool closest(Point1 X[], int n, Point1&amp; a, Point1&amp; b, float&amp; d)
' n* `) V3 I. a: N
& M4 n4 z. {( n* R+ F0 x+ y{// 在n &gt;= 2 个点中寻找最近点对
0 F4 b8 a2 a" T/ @8 y2 m
0 r3 ~) L) A( R1 P" q' ?/ {// 如果少于2个点，则返回f a l s e

" X4 ?% U% Q  o0 P7 s// 否则，在a 和b中返回距离最近的两个点

if (n &lt; 2) return false;
) E- Q) p* s/ ?& v
3 }- G; p' [$ ?% R$ c// 按x坐标排序

' ], e. K9 S0 F' j" S; e8 cM e r g e S o r t ( X , n ) ;

( g7 J8 d7 A. |  q9 a0 q& Y7 k// 创建一个按y坐标排序的点数组
' g$ K' o$ G" O  Y) T; d
# K2 V( F- ~, `3 M, uPoint2 *Y = new Point2 [n];
- j" H; x0 w. V% F
5 h6 o- m( X% M( a6 jfor (int i = 0; i &lt; n; i++) {

// 将点i 从X 复制到Y

Y.p = i;
- M& q" L% U* y% t1 L0 Y( ]) B9 o2 n, g
Y.x = X.x;

Y.y = X.y;
- c! k5 J7 ^4 p( _1 U
}

M e r g e S o r t ( Y,n); // 按y坐标排序

; G! V0 c2 Y5 m+ \0 p. u1 x// 创建临时数组

Point2 *Z = new Point2 [n];
/ y/ J5 y) w. c) @
2 q) X6 ]% A4 `// 寻找最近点对
- H! z2 M" B4 _+ r& I% x- k
c l o s e ( X , Y, Z , 0 , n - 1 , a , b , d ) ;

9 B! P1 k5 L; w( k, y3 W$ E// 删除数组并返回

delete [] Y;

) j& M8 _3 X2 b& z+ odelete [] Z;

4 @* q- @$ l2 M+ N+ J0 a8 n; hreturn true;

2 w. c& S' U5 _; Z7 s; z9 c" Z/ b3 H( x}
' f; D8 ~7 B" x* G2 X
程序1 4 - 11 计算最近点对
, O  N( e7 x' b- f: Q
void close(Point1 X[], Point2 Y[], Point2 Z[], int l, int r, Point1&amp; a, Point1&amp; b, float&amp; d)

$ |7 [. y% A& m( _, L3 J5 }; r{//X[l:r] 按x坐标排序

8 K' S% R9 _/ R9 k( o//Y[l:r] 按y坐标排序
. m, e3 x' R# ^( N* V5 n* ^
! F% A: t; B" I4 X: J+ g8 b  Rif (r-l == 1) {// 两个点
: ^7 }6 J) v; }
. a) ^9 ?( n: L" K  O; G( z9 Pa = X[l];
! f: h+ }1 y( Y/ b
: B4 {" Y) Q7 M6 P" rb = X[r];
) @" {4 f' N4 u( y! a( A+ U3 `/ E
d = dist(X[l], X[r]);

r e t u r n ; }

if (r-l == 2) {// 三个点

8 `; T  y7 U5 w  Q  U) ?) @1 \// 计算所有点对之间的距离

1 f1 W8 ~% }0 w% x# A9 Bfloat d1 = dist(X[l], X[l+1]);

, h, E6 r) v, Q: O/ F4 `( Tfloat d2 = dist(X[l+1], X[r]);

% D/ U/ O9 N" g: I& C! B% l# l( Bfloat d3 = dist(X[l], X[r]);

// 寻找最近点对

if (d1 &lt;= d2 &amp;&amp; d1 &lt;= d3) {
" A$ ]' @4 Q- ~2 ^, e
% s( t$ m, B( V- Fa = X[l];

b = X[l+1];

d = d1;
- q' Y: W" }3 ]; ?. ~8 q1 G8 d7 @+ \
, P+ c9 |7 _" o9 P& ]. x( |r e t u r n ; }

: k1 v9 ^' C1 e' n& ]( C2 aif (d2 &lt;= d3) {a = X[l+1];

8 A8 q8 a7 L$ mb = X[r];
- n! l8 b3 \7 Y! _* o1 F
. O; C+ h5 }9 M3 I3 @! S- P' Zd = d2;}

else {a = X[l];

5 w0 m* r0 [4 m: l, nb = X[r];
# W3 q9 I  }2 e! p# J- Q  C
2 R1 I  H  w4 r; h2 yd = d3;}
0 K' t: f$ r; V& }
r e t u r n ; }
/ ^1 Q9 l& U6 D. o' w
0 L& d7 n9 s5 h$ s' p5 q4 ?/ /多于三个点，划分为两部分
) q* P1 a& }8 ?" i" I
: h3 u0 f: E; _int m = (l+r)/2; // X[l:m] 在A中，余下的在B中
) S; }3 T/ ?" X) f) b9 t4 R# A' @
( e, p0 m* p' Z) c& L// 在Z[l:m] 和Z [ m + 1 : r ]中创建按y排序的表
; ?' [4 D" G0 n7 o* w( R1 M, z
2 u$ w/ H0 U) f& _int f = l, // Z[l:m]的游标
5 r, [$ E* U( F. Y' J
g = m+1; // Z[m+1:r]的游标
! L* `' w/ Q& ^' ^- Y
for (int i = l; i &lt;= r; i++)

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

else Z[f++] = Y;
  T( j6 t/ X; R& H% Q$ A
// 对以上两个部分进行求解
9 }5 t+ I$ ~2 L& s* l- |
5 ^0 c. o# [4 _# D, u4 x  k% jc l o s e ( X , Z , Y, l , m , a , b , d ) ;
" |2 X7 @% j) R4 z, X/ U
; m9 B- o9 ~% ?float dr;

, l. P% u* [- r0 G2 X; J3 p$ ~# rPoint1 ar, br;

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

// (a,b) 是两者中较近的点对

if (dr &lt; d) {a = ar;
( \2 C& x' v* n* Z% i
b = br;

. B) ~0 i* P  Z( Cd = dr;}
; F% ^7 z$ r& m0 R9 D$ I4 W9 P
7 z1 X. o- }* TM e r g e ( Z , Y,l,m,r);// 重构Y

) K" D+ Q" Y- n; h! {1 o! ~; ?/ /距离小于d的点放入Z

int k = l; // Z的游标
0 Q( I, T: b' y
for (i = l; i &lt;= r; i++)
0 j4 e* O2 l7 j9 W1 S& Y/ J$ U" M
if (fabs(Y[m].x - Y.x) &lt; d) Z[k++] = Y;
3 F" r# z0 p' A5 t% r
) ^/ t9 Z6 @+ q/ i. P% t// 通过检查Z [ l : k - 1 ]中的所有点对，寻找较近的点对
1 q# L( n: p8 o4 ?& o
2 E9 _7 O) \9 O' gfor (i = l; i &lt; k; i++){

+ |/ W4 o+ u, |& v) ~$ E# ^; K! n* y; ]for (int j = i+1; j &lt; k &amp;&amp; Z[j].y - Z.y &lt; d;
) j4 p# R; ?0 q3 F: R6 P
j + + ) {

float dp = dist(Z, Z[j]);
2 V8 k6 i- L6 y$ T; i5 w: U* j
if (dp &lt; d) {// 较近的点对
- @( a2 e! @% e. h. v7 t" P
d = dp;

9 `' w: ~  b, s$ e) E  Na = X[Z.p];

: h4 Y! Y# Q" g! f8 E( Z  fb = X[Z[j].p];}
" J& P" F0 c! l1 V0 O# g! E
}

}

6 o9 o' ^& ?9 w* J}

% }* y9 v7 X1 V* f, H8 x" u2 e函数c l o s e（见程序1 4 - 11）用来确定X[1:r] 中的最近点对。假定这些点按x 坐标排序。在Y [ 1 : r ]中对这些点按y 坐标排序。Z[ 1 : r ]用来存放中间结果。找到最近点对以后，将在a, b中返回最近点对，在d 中返回距离，数组Y被恢复为输入状态。函数并未修改数组X。
2 D3 o! w. I4 S9 r& ]% R
首先考察“小问题”，即少于四个点的点集。因为分割过程不会产生少于两点的数组，因此只需要处理两点和三点的情形。对于这两种情形，可以尝试所有的可能性。当点数超过三个时，通过计算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 ]。

0 m  c7 J! S, `& c为实现图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
" n8 F) K% ~1 j7 ]) _
还是在RB中）与Z[j] 配对来实现，其中i＜j 且Z [ j ] . y - Z [ i ] . y＜ 。对每一个Z [ i ]，在2 × 区域内所检查的点如图1 4 - 1 7所示。由于在每个2 × 子区域内的点至少相距。因此每一个子区域中的点数不会超过四个，所以与Z [ i ]配对的点Z [ j ]最多有七个。
6 O7 G0 z. S0 |- l* h
2. 复杂性分析

, v5 L( }0 f8 O. r' ^令t (n) 代表处理n 个点时，函数close 所需要的时间。当n＜4时，t (n) 等于某个常数d。当n≥4时，需花费(n) 时间来完成以下工作：将点集划分为两个部分，两次递归调用后重构Y，淘汰距分割线很远的点，寻找更好的第三类点对。两次递归调用需分别耗时t (「n /2ù」和t (?n /2?).
* z8 }& u' T: \8 e
  R# @' |1 m) E* L! T+ X, z+ L这个递归式与归并排序的递归式完全一样，其结果为t (n) = (nl o gn)。另外，函数c l o s e s t还需耗时(nl o gn)来完成如下额外工作：对X进行排序，创建Y和Z，对Y进行排序。因此分而治之最近点对求解算法的时间复杂性为(nl o gn)。</P>