距离最近的点对

韩冰

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

& I7 \7 R) p+ W. x2 U例14-7 假设在一片金属上钻n 个大小一样的洞，如果洞太近，金属可能会断。若知道任意两个洞的最小距离，可估计金属断裂的概率。这种最小距离问题实际上也就是距离最近的点对问题。
; z( M, p; f: F/ K9 E
通过检查所有的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进行递归求解。
7 V- p! T) F) l/ F, ^% s$ n( q

6 h9 Z/ J3 r1 J$ l( l; u" p7 @0 Oif (n较小) {用直接法寻找最近点对

R e t u r n ; }
* A5 b3 ]; a  u& S2 t; m1 N
# @; v. G! L9 S9 l// n较大
) n- q5 C8 z- X& W
& v9 d  R- t4 Y. m0 @( b' M将点集分成大致相等的两个部分A和B
5 p! T2 m& L' Y1 R
: x1 D" t. Q1 C确定A和B中的最近点对
# _6 i  F6 y% N
确定一点在A中、另一点在B中的最近点对

8 I% M% l  L$ O) v, q) n; q从上面得到的三对点中，找出距离最小的一对点
- {" S( }6 c- y5 y5 ~2 M
) f9 ^& U' |/ R图14-13 寻找最近的点对


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

' G- X- x/ v$ Y5 t* k例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。

设是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) 是最近的点对。

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

1. 选择数据结构
% Y! U7 r) T5 h5 T% M0 C
" [) q% T# @. V. w" m8 F6 J为了实现图1 4 - 1 3的分而治之算法，需要确定什么是“小问题”以及如何表示点。由于集合中少于两点时不存在最近点对，因此必须保证分解过程不会产生少于两点的点集。如果将少于四点的点集做为“小问题”，就可以避免产生少于两点的点集。
- i2 O$ ~# ]) N
每个点可有三个参数：标号， x 坐标，y 坐标。假设标号为整数，每个点可用P o i n t l类（见程序1 4 - 8）来表示。为了便于按x 坐标对各个点排序，可重载操作符＜=。归并排序程序如1 4 -3所示。

/ E4 N8 z  g' \& }& D  ~程序14-8 点类
- U0 F! W! b3 U0 S
class Point1 {

  }) N- U  H( X' Q5 H' ffriend float dist(const Point1&amp;, const Point1&amp;);
: G4 w. D6 u' o4 E- y2 ~
# V4 }( w. d7 K: ?friend void close(Point1 *, Point2 *, Point2 *, int, int, Point1&amp;, Point1&amp;, float&amp;);
( t  k% n( m) s0 U4 `, |+ k
/ V! o% {; f6 U8 h" l/ Z# x9 rfriend bool closest(Point1 *, int, Point1&amp;, Point1&amp;,float&amp;);
1 b- x9 h5 U7 N! _/ c+ }+ X
friend void main();

p u b l i c :

1 D$ m0 u/ i) ?. Yint operator&lt;=(Point1 a) const

{return (x &lt;= a.x);}
+ v; g: ?* b& Z( A! J5 \7 O
p r i v a t e :
3 S3 n' z6 U5 \" j- ^
int ID; // 点的编号
. [1 h: C% i9 @$ L5 n$ P% l
3 m4 p) M/ C! J# \float x, y; // 点坐标
2 S6 w. y; f2 ?% E) P5 _" h
% K3 W, j+ T# e8 R* r} ;

2 d$ d, d! M# M7 ?  d. N5 Aclass Point2 {
% ~7 V+ \: M2 G7 u* H
) Y$ V8 l& a- R; N9 {( ?9 efriend float dist(const Point2&amp;, const Point2&amp;);
; ], M# A! K4 }9 o4 }
friend void close(Point1 *, Point2 *, Point2 *, int, int, Point1&amp;, Point1&amp;, float&amp;);

% T2 T% Y6 U$ V) U6 `friend bool closest(Point1 *, int, Point1&amp;, Point1&amp;, float&amp;);

friend void main();

" }3 ]" D4 j7 R7 up u b l i c :

9 B+ G* ?4 ?* M. b8 C* Tint operator&lt;=(Point2 a) const
5 P, Z; O8 }2 D3 w
. {9 C+ O; h: U+ x$ x; O% n{return (y &lt;= a.y);}
4 W& U* ?  U8 l1 j1 k
4 ?& @) E5 _7 [, Gp r i v a t e :

int p; // 数组X中相同点的索引
" r9 T8 F' n/ g/ ~, @) |# `
  c* F' p- E; m- q5 p. U+ efloat x, y; // 点坐标

0 n1 q5 I* z6 {5 b1 @} ;

. `9 \0 ]& g5 W2 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中的对应点。

, l7 A, M$ X% K: d* 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类的友元。
! X8 L+ R9 U! K, L) p  t$ w
$ ?) l! s- I- z! c/ X0 H4 _  S程序14-9 计算两点距离

4 b3 p  X$ m( T# m# y+ t+ ?template<CLASS T>
  Z0 M4 G& y" l! w# c/ h
inline float dist(const T&amp; u, const T&amp; v)

- m* z3 M* X% \3 G1 y- B+ O{ / /计算点u 和v之间的距离
4 B9 e3 h/ @9 C0 [1 V- C+ y
- j; O, k# x5 D, y! V7 G6 Y6 |float dx = u.x-v. x ;
" x5 y9 P  J+ h/ R: @! }  Q. ]
float dy = u.y-v. y ;

# {/ R- y) E1 X+ r3 q( |5 ^return sqrt(dx * dx + dy * dy);

}
9 b/ s" b0 w# r; _: d; l; L0 ~
+ m( a1 s; j  y  f6 Q1 {/ b% S3 a如果点的数目少于两个，则函数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 )，该函数实际求解最近点对。
3 R( [; K: i9 a! i! K/ O
程序14-10 预处理及调用c l o s e
3 B8 A, b+ T# a% _2 }
1 i5 l  ~, K3 ~  U* v' }! Vbool closest(Point1 X[], int n, Point1&amp; a, Point1&amp; b, float&amp; d)
" v- V) [' p+ d% h
- a  {" T5 Y+ l{// 在n &gt;= 2 个点中寻找最近点对
' T* X3 R% ?1 f9 M0 d& r" o% a
// 如果少于2个点，则返回f a l s e
3 P# B5 Y* a, I4 n" Z7 J3 Z
9 S2 Q# b: \: U$ ]( f2 J// 否则，在a 和b中返回距离最近的两个点
; X0 E* ~5 ~) r# p8 @/ y4 i- X
; B! |; f" p% z0 E  Xif (n &lt; 2) return false;

# G% T3 U. z+ \6 l// 按x坐标排序
# c$ a  a0 i' Y* d
M e r g e S o r t ( X , n ) ;

' @" i" {  N: W3 S8 q" V// 创建一个按y坐标排序的点数组
7 R$ [. e/ N; e
, T/ ]& `& G8 T/ XPoint2 *Y = new Point2 [n];

. X& b2 E/ e  Y9 ^3 i3 Yfor (int i = 0; i &lt; n; i++) {

// 将点i 从X 复制到Y

Y.p = i;
0 q& d' _# t! R* s( T$ b8 J
Y.x = X.x;
& H" V5 V% g: D
& y. k) U. T0 f. H' ?Y.y = X.y;

}
# F- {# M! \3 O: Z: @6 X. Q
M e r g e S o r t ( Y,n); // 按y坐标排序
; ~6 z4 B, a& B7 p
/ f" c7 Y; H5 B  S// 创建临时数组
! G( C/ Z- \% p& s' C9 f: v, F; \
Point2 *Z = new Point2 [n];

// 寻找最近点对

c l o s e ( X , Y, Z , 0 , n - 1 , a , b , d ) ;
7 m  P( O- Q& d
# W9 K3 ]) Q; |  ?3 W) H2 S+ B// 删除数组并返回
3 B8 i& l6 _6 R6 Z( L6 e0 k
% |2 L5 D! f0 f: xdelete [] Y;
- ^& f2 N$ a+ B' ^# o
delete [] Z;

* a, `% E# M% K' l, x, Dreturn true;

}

程序1 4 - 11 计算最近点对

4 @! v- \+ y7 a' G5 zvoid close(Point1 X[], Point2 Y[], Point2 Z[], int l, int r, Point1&amp; a, Point1&amp; b, float&amp; d)
+ c+ `) ]8 G# ^0 o
{//X[l:r] 按x坐标排序

//Y[l:r] 按y坐标排序

7 _9 ]3 {; @6 g7 c; g) |7 cif (r-l == 1) {// 两个点

5 I( K' X( I6 d6 ^+ _! x, Fa = X[l];
8 v4 ?/ S' p2 C0 ^
b = X[r];

d = dist(X[l], X[r]);
# P* Q9 h8 {' r% _
9 V9 H$ `* ~( j( u4 K5 F* H) wr e t u r n ; }

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

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

/ Y# c; O8 P+ \* d0 e: Tfloat d1 = dist(X[l], X[l+1]);
3 w9 u+ N0 N+ O8 C) G2 ?3 c
* r7 V; c$ ?3 c2 L9 ]( S  |float d2 = dist(X[l+1], X[r]);
0 v. m- e. ?# ]
float d3 = dist(X[l], X[r]);
5 i; s, U9 a( v0 C3 L; y* ~. o# H
// 寻找最近点对
1 t9 P& o$ g, ^1 q/ m
6 K- A2 |# x; ?' Gif (d1 &lt;= d2 &amp;&amp; d1 &lt;= d3) {

( u( y5 L) ]$ y' s8 k& pa = X[l];

b = X[l+1];

d = d1;
. [+ A: O/ V) w
- G: T# `3 _9 h* x! Er e t u r n ; }
# }5 Q4 u! D* X; q& }
if (d2 &lt;= d3) {a = X[l+1];
/ [, Y! C% p" s! o* f$ f0 a: h* J
7 h. d9 l/ x( @8 P; h1 Zb = X[r];

d = d2;}
% p4 h  o  z9 C7 R$ q. v' Q5 ]
( o  s- u: e, `else {a = X[l];

b = X[r];
2 U6 I$ [& Z7 }0 F7 Y  s3 m4 @
d = d3;}
+ ], R# O& }, |
% q! `/ M6 ]# O0 rr e t u r n ; }

4 s6 j4 U+ _# E4 N6 w/ /多于三个点，划分为两部分

' v) p& j) y% C( W/ s; g* Aint m = (l+r)/2; // X[l:m] 在A中，余下的在B中

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

4 J+ J; `' J' z5 Jint f = l, // Z[l:m]的游标

# ]0 W: K/ D! {; C% L5 E* b2 ]% _! x$ ?g = m+1; // Z[m+1:r]的游标
+ Z+ r; z( N- b* D+ @, u; l3 c) a$ Z
for (int i = l; i &lt;= r; i++)

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

else Z[f++] = Y;

, ^, o) ^7 w- w! y& F% K7 n2 R// 对以上两个部分进行求解
$ S( Y& p- d4 l8 I5 u) z* T
2 n- G$ N0 F( Z/ rc l o s e ( X , Z , Y, l , m , a , b , d ) ;

float dr;
  ^1 T' s* I  L8 n7 w/ v4 |
Point1 ar, br;

, C' I  ?% f+ N& \c l o s e ( X , Z , Y, m + 1 , r, a r, b r, d r ) ;
9 V5 T' z( z5 t! x' M) Q# G- V
// (a,b) 是两者中较近的点对

if (dr &lt; d) {a = ar;
3 M6 u+ O. n9 M) I/ B0 M
+ s8 v3 i2 x' U& E0 Y% Jb = br;

d = dr;}

M e r g e ( Z , Y,l,m,r);// 重构Y
9 O) `7 S* L- n) [) K, m- L
  I4 ?# z, M) d/ /距离小于d的点放入Z
/ y; v% G- W# q! D9 }& N+ a2 z- Y
int k = l; // Z的游标
0 ]' V. u. c1 E* a7 I8 q7 |
for (i = l; i &lt;= r; i++)
) h5 }- |8 s& \1 O' f1 l* U. F9 R
) t% }2 _7 t9 f5 m( l8 G+ Hif (fabs(Y[m].x - Y.x) &lt; d) Z[k++] = Y;
" f7 T; R( F4 v
. r/ Q( V6 O. [- E- _  s& v// 通过检查Z [ l : k - 1 ]中的所有点对，寻找较近的点对

for (i = l; i &lt; k; i++){

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

7 M# T  D4 {9 i! l( A. w/ U4 l) cj + + ) {

: C5 A, I) l2 V" T6 Q. p1 ^float dp = dist(Z, Z[j]);

7 d1 m- ?% C9 v; r0 L+ x% h/ uif (dp &lt; d) {// 较近的点对
1 J  I) b( U5 M- I2 R5 q
d = dp;

* w9 B+ h) N, U4 x: Sa = X[Z.p];
2 n, O. D/ L% E! X( d
b = X[Z[j].p];}

}
" g9 w! n  A3 ~* K0 [( s
" o  t+ Y. s( d; x$ P}
0 k# J: G( T! t6 S8 N9 Y6 u6 ~
( e% M- X0 x( N) f}

3 [, f% t' u6 ~* O4 b函数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 ]。
' q& e# c" x- a
, o0 D, _6 u: u. 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 ]最多有七个。
3 c. O7 M1 V( S* c5 V7 s; N
2. 复杂性分析

- j- o# g. O* }; t. t1 L6 a3 K2 z: C令t (n) 代表处理n 个点时，函数close 所需要的时间。当n＜4时，t (n) 等于某个常数d。当n≥4时，需花费(n) 时间来完成以下工作：将点集划分为两个部分，两次递归调用后重构Y，淘汰距分割线很远的点，寻找更好的第三类点对。两次递归调用需分别耗时t (「n /2ù」和t (?n /2?).

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