距离最近的点对

韩冰

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

例14-7 假设在一片金属上钻n 个大小一样的洞，如果洞太近，金属可能会断。若知道任意两个洞的最小距离，可估计金属断裂的概率。这种最小距离问题实际上也就是距离最近的点对问题。

, _/ g( |: [; x8 j+ J  F% D通过检查所有的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进行递归求解。
* ^" F6 ~# K: n3 o* W

+ G- L, F2 y- S: r7 e& fif (n较小) {用直接法寻找最近点对
+ j( y- L+ u: d7 R  c* q0 P
R e t u r n ; }

, `; y7 ]4 _' i4 o  x// n较大
3 Q# ]' f1 o: Y7 [! a
( s3 i5 w7 B1 z7 M# y- F4 R5 e将点集分成大致相等的两个部分A和B

) ~; D4 i  Y4 `2 w) E6 \8 J确定A和B中的最近点对

& \) j0 Z5 k1 N: V3 t确定一点在A中、另一点在B中的最近点对

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

; X" O; M4 Z5 K, Q' s) C图14-13 寻找最近的点对

. _. V7 G5 S+ ?
4 z2 h& P3 @5 D% u为了确定第三种情况下的最近点对，需要采用一种不同的方法。这种方法取决于点集是如何被划分成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。
9 S+ ~+ X' z- S' u3 w: l# B6 S: q
, B, k1 V0 v; U, f设是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）。

( g) p5 v  {& ]6 V: p例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) 是最近的点对。
  B# O/ H+ M6 E/ b
为了确定一个距离更小的第三类点，RA 中的每个点最多只需和RB 中的6个点比较，如图1 4 - 1 6所示。
  G7 P! f: b8 V$ X! n' }( A
1. 选择数据结构

! f/ {7 s0 b9 t5 P/ @: ]% I" h, G为了实现图1 4 - 1 3的分而治之算法，需要确定什么是“小问题”以及如何表示点。由于集合中少于两点时不存在最近点对，因此必须保证分解过程不会产生少于两点的点集。如果将少于四点的点集做为“小问题”，就可以避免产生少于两点的点集。
8 P, \( A1 e1 u$ U) O
# g+ ~4 G+ h. e# G  B每个点可有三个参数：标号， x 坐标，y 坐标。假设标号为整数，每个点可用P o i n t l类（见程序1 4 - 8）来表示。为了便于按x 坐标对各个点排序，可重载操作符＜=。归并排序程序如1 4 -3所示。
1 \) u* J! D% F, A: T* y+ z; A3 f
; D, r9 x  l2 f, n( }$ ?! P程序14-8 点类
' ]$ n. p4 x* l& |7 h
! H7 |- ?& M7 ~* `5 u# qclass Point1 {

1 ]7 R6 L- V9 W5 Y6 c: Ffriend float dist(const Point1&amp;, const Point1&amp;);

friend void close(Point1 *, Point2 *, Point2 *, int, int, Point1&amp;, Point1&amp;, float&amp;);

4 |( I3 s4 |( B2 z* ]: ]( }( tfriend bool closest(Point1 *, int, Point1&amp;, Point1&amp;,float&amp;);
4 S) V8 X: i' N4 X, C" \% c1 p
9 L/ {! l+ R7 E4 M9 Yfriend void main();

. n; L: e2 s3 ~3 W' h0 Kp u b l i c :
' n4 c8 i) ]' `. }* S) J
5 S* ~  i8 r$ r+ y' l, Bint operator&lt;=(Point1 a) const
# T; z, X: `4 j( T' X  M; a; f: d
. d+ ?& B& O9 B3 i+ Q+ T( k  q" `/ ]{return (x &lt;= a.x);}
' y3 k9 b: y& m6 F
p r i v a t e :

  i# \" l* I" U1 aint ID; // 点的编号

float x, y; // 点坐标

} ;
+ |, A& F+ ~/ M& [# ]$ R9 X
class Point2 {

7 \8 S; R) Z! ]7 r6 E8 E. x( {+ Mfriend float dist(const Point2&amp;, const Point2&amp;);
% v' V6 @. ]) T  A
friend void close(Point1 *, Point2 *, Point2 *, int, int, Point1&amp;, Point1&amp;, float&amp;);
3 |* T3 h. S8 N+ T
friend bool closest(Point1 *, int, Point1&amp;, Point1&amp;, float&amp;);

friend void main();
5 P. i# r/ [: {- p
p u b l i c :
  {% V$ j& S6 p
int operator&lt;=(Point2 a) const
# \' B1 ]1 l0 S$ }  Y
3 y* \( t3 {7 v- d- k2 ^' N{return (y &lt;= a.y);}

p r i v a t e :

int p; // 数组X中相同点的索引
! A0 O" J, g5 ]0 F: i3 d/ {
float x, y; // 点坐标
5 H: ^: n& S& B! d% I* x" r
} ;
/ R: d* C4 H" c7 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中的对应点。
' w& ]; r; n- d" h& l0 z9 `
$ V8 U3 R8 g# d- v, Q. m& Y确定了必要的数据结构之后，再来看看所要产生的代码。首先定义一个模板函数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类的友元。

5 Z3 c  q$ Q; C& B  w4 Z- I程序14-9 计算两点距离

template<CLASS T>
2 J- x" _: v% q7 H
inline float dist(const T&amp; u, const T&amp; v)
( i; m6 }+ g9 n& |
3 i* L" }% }% K- U7 ?{ / /计算点u 和v之间的距离

8 ?7 Z1 z4 |2 N  d7 @6 kfloat dx = u.x-v. x ;

float dy = u.y-v. y ;
" y6 N- L3 l# r
return sqrt(dx * dx + dy * dy);

}
' t, v8 p, w9 z( Y5 d6 L. k6 E
) W# n3 z: E8 ^2 K% a: R" g5 u如果点的数目少于两个，则函数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
, A2 Y5 e$ E3 C- M$ M& w# X0 Y
bool closest(Point1 X[], int n, Point1&amp; a, Point1&amp; b, float&amp; d)

) r6 w% o. U4 p, f{// 在n &gt;= 2 个点中寻找最近点对

// 如果少于2个点，则返回f a l s e

0 H% r  }$ j/ ~1 k8 ~/ k// 否则，在a 和b中返回距离最近的两个点

6 y$ I# w6 |' i* hif (n &lt; 2) return false;

// 按x坐标排序
* W: n. \1 g5 X2 b3 p2 R# p6 H
! `9 S* ?& ?3 r& `* y& bM e r g e S o r t ( X , n ) ;
0 G0 B6 U) o& @3 I' N. ~" _% J3 S/ u
; y0 h+ X" j' Y' n% a// 创建一个按y坐标排序的点数组

1 i& G" [3 W9 h' FPoint2 *Y = new Point2 [n];

/ X' _0 _$ M# w3 I, |% {$ P  u8 Rfor (int i = 0; i &lt; n; i++) {
: j3 x) J3 z4 j7 i, j- Q- {8 p
4 G1 P' g6 V" f7 x// 将点i 从X 复制到Y
8 I* r, h7 }9 Y  Q7 p/ E
% Y( h. U+ T# W/ mY.p = i;
; J" T9 q# ~- c/ m: e' O9 i
) \; w  c. x/ e1 ?6 L4 R  i- n/ X, {Y.x = X.x;

2 L% U9 X( o* L- G0 HY.y = X.y;
1 U  h* @' `* o
}

' b2 N9 |# h' _. _% k/ U8 g6 T: }M e r g e S o r t ( Y,n); // 按y坐标排序

. |8 V4 E! y7 V// 创建临时数组
' o/ I9 D: K0 {3 K5 H% @
$ A1 c6 C' z. F" bPoint2 *Z = new Point2 [n];
2 i, ^9 N- z" d
; G  ?( Z% Z4 B: i- z, s- A// 寻找最近点对
6 K; u% t) `' z2 }4 h8 w
9 g/ l% }- S5 \0 g$ vc l o s e ( X , Y, Z , 0 , n - 1 , a , b , d ) ;
5 l+ U- A; E4 A
// 删除数组并返回

, W) h1 Q1 R' L2 p& Y  {delete [] Y;
% @' v1 l* p) P$ P  }* }: T/ W
delete [] Z;
+ m: G, f% Z' y+ t8 ?( n
# ?' D/ _" C, X; Ereturn true;

* K+ L/ N8 v, p) K0 q# T$ H4 }}
- d% S* K. y! Y" n
程序1 4 - 11 计算最近点对
1 I* c. G/ ^: k- i# p' F& q
void close(Point1 X[], Point2 Y[], Point2 Z[], int l, int r, Point1&amp; a, Point1&amp; b, float&amp; d)
" w8 b) j6 V; g4 s2 E. G/ m
{//X[l:r] 按x坐标排序
  z- G: k' V% r' Z
+ n& C% n% a" y- a1 s. Z* T' X//Y[l:r] 按y坐标排序
' z) L) y' F5 L, Y- R( m! `
* Y2 |' Z7 y5 ^1 o5 r/ A2 ~if (r-l == 1) {// 两个点

3 q3 u  X  g5 P* Xa = X[l];

0 l3 b4 v: m1 R1 R5 Z# {b = X[r];
- k5 J5 e) C4 q& y; V
. |$ @5 x+ ?: Q3 b4 p! @/ Vd = dist(X[l], X[r]);
: _1 D! n& e- T
8 u( G$ @' k* }  j9 qr e t u r n ; }

! J0 t; f4 A* u0 o( H0 Dif (r-l == 2) {// 三个点

9 a, a7 H" f* F3 `// 计算所有点对之间的距离
( [" K1 x1 L  {# m1 F
; Q" k9 z5 x5 \float d1 = dist(X[l], X[l+1]);

float d2 = dist(X[l+1], X[r]);
% A' [% d7 Y1 w3 t, a
* R3 ^' P: \) W7 [% Ifloat d3 = dist(X[l], X[r]);

# [& K; S6 v( b& s" P6 ~// 寻找最近点对

/ u& `$ O1 D9 y9 O( ^" v% Iif (d1 &lt;= d2 &amp;&amp; d1 &lt;= d3) {

" [, C! Y. W7 [a = X[l];
7 }/ G% L6 b& U
5 ?6 W) q! F% E/ @b = X[l+1];
3 R5 l  c) g* e0 Q' k
d = d1;

r e t u r n ; }

if (d2 &lt;= d3) {a = X[l+1];
7 q# ]" t' {3 q. b2 n. u
, N$ V9 q( q/ W4 }' i+ ~b = X[r];
2 y0 t% g$ J0 x: O3 @
d = d2;}
; R0 l& K9 E1 J4 V3 r
: n5 E) U& }. Z2 \: Velse {a = X[l];
4 l" w) m6 J; f- u. B
b = X[r];

3 e: a1 i6 E) r) |$ g3 c% ^, @" Id = d3;}

+ J- H2 m- \) N0 @: m7 Cr e t u r n ; }
6 F: k9 M  J$ l; ]- k
  h" D1 H) O! V! O! f1 k/ /多于三个点，划分为两部分

8 b1 s4 \: w1 a9 _1 x1 Eint m = (l+r)/2; // X[l:m] 在A中，余下的在B中
( c5 D5 |4 Q- Q+ G' h5 e
5 y  M1 p) ^& p+ P- {// 在Z[l:m] 和Z [ m + 1 : r ]中创建按y排序的表

int f = l, // Z[l:m]的游标
4 S% m9 ?3 c* I+ K4 l: p
g = m+1; // Z[m+1:r]的游标
! X6 s9 H% Z* c: w' Y4 X9 C
+ D& ]0 f  Y6 P9 ~/ ufor (int i = l; i &lt;= r; i++)
( L4 g2 K! Q4 _# a# X' ?+ p, r6 i
if (Y.p &gt; m) Z[g++] = Y;
$ t8 x% z9 E' U1 A1 g8 \3 K/ h! N- A
8 b0 [2 E5 U- Oelse Z[f++] = Y;
$ u2 ^. o! Z8 T9 X0 x
// 对以上两个部分进行求解
9 A! u2 [. }9 F/ t6 h
& L4 M) F5 `9 j& w4 Sc l o s e ( X , Z , Y, l , m , a , b , d ) ;
) ^+ [7 b2 `" R( W7 X; j0 F1 v& Q" g
  @$ N# N6 M. N7 u* N8 [: Afloat dr;
- j4 a" e3 Q% U4 W8 r4 E, a/ w3 D
1 ^  j1 x; {5 k$ M3 F, {Point1 ar, br;

3 P2 _, d& ]% E% g3 Jc l o s e ( X , Z , Y, m + 1 , r, a r, b r, d r ) ;
* R# V! ]2 C8 `6 r$ u8 e. m% M
// (a,b) 是两者中较近的点对
" L' J5 v! @/ J  m8 Y
if (dr &lt; d) {a = ar;

b = br;
- ], }6 ?/ I6 ?* g
) a- p+ e. Z4 \3 X# `3 u8 Ud = dr;}

M e r g e ( Z , Y,l,m,r);// 重构Y

3 G$ s1 a. y* D7 U1 _/ /距离小于d的点放入Z
( F3 m; N9 D$ k* V
5 \: \" V) N' h. _/ nint k = l; // Z的游标

3 o6 ~  P, d( u- f- zfor (i = l; i &lt;= r; i++)
5 _& u0 a6 ~$ X0 d% `. T  ~  _) g
if (fabs(Y[m].x - Y.x) &lt; d) Z[k++] = Y;

7 c3 M1 m# f; L. e  L" w// 通过检查Z [ l : k - 1 ]中的所有点对，寻找较近的点对

& R. `. P  P8 e5 p' [for (i = l; i &lt; k; i++){
0 m: |+ x8 A6 D6 @5 ?. G5 N
  C+ f7 G' e% S4 l7 }for (int j = i+1; j &lt; k &amp;&amp; Z[j].y - Z.y &lt; d;
7 y1 v- |( K8 S
& r+ f+ [+ C* R3 i7 w3 Zj + + ) {
+ B7 B5 d+ @8 z$ @
- H+ ?' a7 e" u8 zfloat dp = dist(Z, Z[j]);
4 d! ^- u; s2 Y. z: ^
if (dp &lt; d) {// 较近的点对

5 z+ d9 }" W- P" j, H/ s3 F9 R3 ?! Od = dp;

5 f9 l  y& ~" D$ d, C/ ma = X[Z.p];

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

}

0 A" e. g* T, o}

}

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

  P  o1 k& j% I7 C3 E首先考察“小问题”，即少于四个点的点集。因为分割过程不会产生少于两点的数组，因此只需要处理两点和三点的情形。对于这两种情形，可以尝试所有的可能性。当点数超过三个时，通过计算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 L: G9 r' A' z7 B, U
为实现图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 S7 q: B2 Q+ i1 ~! I还是在RB中）与Z[j] 配对来实现，其中i＜j 且Z [ j ] . y - Z [ i ] . y＜ 。对每一个Z [ i ]，在2 × 区域内所检查的点如图1 4 - 1 7所示。由于在每个2 × 子区域内的点至少相距。因此每一个子区域中的点数不会超过四个，所以与Z [ i ]配对的点Z [ j ]最多有七个。
) H" B9 l8 K- h# C! R- A! F7 W
% @( H% A& v, b) N) j1 r2 F2. 复杂性分析
$ c5 p9 y" b7 a$ L7 V3 Z
令t (n) 代表处理n 个点时，函数close 所需要的时间。当n＜4时，t (n) 等于某个常数d。当n≥4时，需花费(n) 时间来完成以下工作：将点集划分为两个部分，两次递归调用后重构Y，淘汰距分割线很远的点，寻找更好的第三类点对。两次递归调用需分别耗时t (「n /2ù」和t (?n /2?).

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