距离最近的点对

韩冰

<>给定n 个点（xi，yi）（1≤i≤n），要求找出其中距离最近的两个点。
$ M+ w: u$ L7 X" n5 Q
& U) V% `' E6 f5 i例14-7 假设在一片金属上钻n 个大小一样的洞，如果洞太近，金属可能会断。若知道任意两个洞的最小距离，可估计金属断裂的概率。这种最小距离问题实际上也就是距离最近的点对问题。
* U) c8 X2 z+ R: F$ G
' t7 Q. y9 H' ~7 a% z; v% J6 A通过检查所有的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进行递归求解。

$ Q! Y- b, N9 Y$ v( r
if (n较小) {用直接法寻找最近点对
7 S$ I0 K1 n  q9 ]; }3 L5 S
R e t u r n ; }
+ n. ?2 H% R+ M4 a- p$ b  r) D
+ y; W% L7 V2 u- V5 Y. T! e// n较大

& D' e5 C. b5 [将点集分成大致相等的两个部分A和B
3 D4 g, V' H; n0 b6 P* J
: H* [9 T9 Q( y3 ~: s确定A和B中的最近点对

确定一点在A中、另一点在B中的最近点对
! I* u0 v$ k8 y# T+ m9 {) X/ Z; `+ \
从上面得到的三对点中，找出距离最小的一对点

, w8 |3 c, A* k" y% x图14-13 寻找最近的点对
: u4 O4 j7 d& z7 _3 V
4 ?: r- B( S. u! X; B
' Z+ G$ G5 B6 o% \1 w1 s+ F- C为了确定第三种情况下的最近点对，需要采用一种不同的方法。这种方法取决于点集是如何被划分成A、B的。一个合理的划分方法是从xi（中间值）处划一条垂线，线左边的点属于A，线右边的点属于B。位于垂线上的点可在A和B之间分配，以便满足A、B的大小。
4 v$ e* R, l+ }) v
& F. |( B. f. q. ~例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 S5 H2 `% W  {
设是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）。
* q. K0 [5 c& s0 D1 j6 a: z
9 K: ]4 ~% _  h9 Y: n/ z$ N3 F例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所示。
5 U9 _$ u4 E9 o0 h$ ~( f
# F1 v) D% B5 H1. 选择数据结构

为了实现图1 4 - 1 3的分而治之算法，需要确定什么是“小问题”以及如何表示点。由于集合中少于两点时不存在最近点对，因此必须保证分解过程不会产生少于两点的点集。如果将少于四点的点集做为“小问题”，就可以避免产生少于两点的点集。

每个点可有三个参数：标号， x 坐标，y 坐标。假设标号为整数，每个点可用P o i n t l类（见程序1 4 - 8）来表示。为了便于按x 坐标对各个点排序，可重载操作符＜=。归并排序程序如1 4 -3所示。
2 l: w0 Y  s; o9 ~
8 Y7 H+ B7 p+ R* _程序14-8 点类
; L9 |* Q$ V6 q; K$ n4 G0 v
class Point1 {

* r+ s! v  }- p" f# s4 ffriend float dist(const Point1&amp;, const Point1&amp;);
% V/ L% t6 J1 \# l" E$ \3 S( i
friend void close(Point1 *, Point2 *, Point2 *, int, int, Point1&amp;, Point1&amp;, float&amp;);

! |$ J# m* @, A: Sfriend bool closest(Point1 *, int, Point1&amp;, Point1&amp;,float&amp;);

friend void main();
) f8 j+ k3 h7 M- `/ q5 E$ `2 f! n6 S
8 c* H" P8 s8 L) u7 yp u b l i c :
' t) q6 Y4 P& m3 x0 P
" \1 f# L% ]8 D; H4 J  \/ uint operator&lt;=(Point1 a) const
: H  p7 n; N; ]+ A# y
{return (x &lt;= a.x);}
/ D" i6 Y+ B% M
p r i v a t e :
4 c8 \7 s* F0 {( K
int ID; // 点的编号

float x, y; // 点坐标
: V, r$ l2 w& R; B
} ;

: o2 R2 W$ x2 P! s) xclass Point2 {
/ ^2 {, F3 u6 o( F9 y' y" D
friend float dist(const Point2&amp;, const Point2&amp;);
' i' t* z+ ^( e# n8 z5 j
2 [+ ]7 V% V; Ufriend void close(Point1 *, Point2 *, Point2 *, int, int, Point1&amp;, Point1&amp;, float&amp;);
7 |, F$ C* r# c& Z, M, D" I$ Y
% L2 R1 G6 o& P% _( q2 `friend bool closest(Point1 *, int, Point1&amp;, Point1&amp;, float&amp;);

friend void main();

3 P' }, |  N9 P+ o% Q5 fp u b l i c :
+ z* j' O' I( t  N; ?
. G8 R) R1 r3 h% o% @% t9 Sint operator&lt;=(Point2 a) const

{return (y &lt;= a.y);}
. Z, [5 V: _) c/ d; g4 z! W9 P
p r i v a t e :
9 M" Z3 m# ^& }* V/ Q
. D3 \; q. ]  U- `8 Pint p; // 数组X中相同点的索引
. i! }+ p: z3 e- H% \4 o' E" A
float x, 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中的对应点。
- F  l: D4 C3 m) y2 r3 D
确定了必要的数据结构之后，再来看看所要产生的代码。首先定义一个模板函数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 计算两点距离

( y; i" t/ c. S- T, r5 Ntemplate<CLASS T>
6 ?1 S( c8 x1 |* p4 k- K# R% [
inline float dist(const T&amp; u, const T&amp; v)
0 c, I2 l( w% |; o4 ]: v* U
" H& Q2 N5 k9 _: F{ / /计算点u 和v之间的距离
6 d! _3 N! D2 r6 ^( C3 o
* D/ N8 h4 e) y3 U) u8 q% @) `float dx = u.x-v. x ;

float dy = u.y-v. y ;
/ B7 Y/ ?, z. W) H4 a) i* {( j
/ `" T5 j) x7 V6 ]return sqrt(dx * dx + dy * dy);
2 A: R5 u' _, X/ V" ?9 I
}

6 e9 L5 x( z9 l8 c" e0 j/ S如果点的数目少于两个，则函数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 )，该函数实际求解最近点对。
" N( c6 B/ A# G2 V6 F( `& f% f
程序14-10 预处理及调用c l o s e
9 H) i4 l0 {. J5 L9 E( J
7 ^7 K( a8 h% S+ N$ c4 {; H, e6 Z3 Fbool closest(Point1 X[], int n, Point1&amp; a, Point1&amp; b, float&amp; d)
& Q& g4 a/ x: i- o4 S" W0 I7 c$ _
. V0 q1 r" G8 x# C1 _{// 在n &gt;= 2 个点中寻找最近点对

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

0 U' k  l! Q7 F8 S. J. \// 否则，在a 和b中返回距离最近的两个点
3 A" H/ f+ l6 F$ \7 ~
if (n &lt; 2) return false;
9 r. }! Q) R4 l
// 按x坐标排序
+ V4 `! m7 Q2 c
, g4 n  {* K. q3 V5 C" ~: j' i4 LM e r g e S o r t ( X , n ) ;
8 v: i5 j' ~7 c" m, |
// 创建一个按y坐标排序的点数组

Point2 *Y = new Point2 [n];
2 \+ P  `7 W9 I- y: i. v/ w# F5 c
for (int i = 0; i &lt; n; i++) {
9 u, _6 X- F% o5 p; c; A
// 将点i 从X 复制到Y
. o: ]9 D& q9 I# a6 Z0 K( _
Y.p = i;

; C% }& }( q% LY.x = X.x;
' G7 a8 ^3 N. c% x2 I" I
, I; p/ G! q) w% OY.y = X.y;

' o6 B) w5 z9 l" y/ z}

M e r g e S o r t ( Y,n); // 按y坐标排序
" b; N( Z- z, u* F: }8 `. [+ V
// 创建临时数组
8 A2 z3 R. k( c/ z4 s, v& Y$ _
Point2 *Z = new Point2 [n];

; E2 C- [1 P# Y) a0 J6 s5 J% p2 X0 O// 寻找最近点对
, Z( ?$ g4 M9 G3 b
" ?/ M5 g4 K. T( sc l o s e ( X , Y, Z , 0 , n - 1 , a , b , d ) ;
, i2 o, \0 I. v
// 删除数组并返回
- m2 o4 s4 w; Q( f3 B1 E
3 o' Z8 q' Y( w5 t1 @delete [] Y;

delete [] Z;

& r! T) h' x6 w2 `+ _return true;
. a# s' O* I7 _& p: }$ e8 G% j( g
}

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

void close(Point1 X[], Point2 Y[], Point2 Z[], int l, int r, Point1&amp; a, Point1&amp; b, float&amp; d)
, |. e9 k: C3 c2 C, g
{//X[l:r] 按x坐标排序
& j/ g3 u' Q; a" k- U# ~
//Y[l:r] 按y坐标排序

6 X7 {8 a6 M9 _if (r-l == 1) {// 两个点

! J! i$ Z! r1 v" Y, Va = X[l];
9 w7 L* s8 M/ D7 X6 }" t6 E
4 |7 N" `6 X7 l- C9 Bb = X[r];

$ [) J5 k. d& hd = dist(X[l], X[r]);

+ [6 u; [7 l5 }! Tr e t u r n ; }

+ ~" u% O; A2 b6 t) {& J+ J- hif (r-l == 2) {// 三个点
6 r$ P3 F8 o/ G# c7 {7 _  j
" X( \  P9 d" d8 \: R// 计算所有点对之间的距离

( ]& |4 \4 }6 W/ _3 v$ P- l2 Bfloat d1 = dist(X[l], X[l+1]);
5 e6 n& X! e0 [7 ]
float d2 = dist(X[l+1], X[r]);

+ a$ N& K! z3 m  Ofloat d3 = dist(X[l], X[r]);

! ]# {/ z( z$ d7 N5 Z9 L// 寻找最近点对
2 q. a# p9 c+ y& b. F) W1 |2 W
3 u! S+ S7 |6 _7 Oif (d1 &lt;= d2 &amp;&amp; d1 &lt;= d3) {

a = X[l];
0 f# d- B* t" v7 w# d) o9 h2 T
b = X[l+1];
, m4 L: C% d- q7 W: W, o
' W; [. I  T$ C6 Ed = d1;

: `( U; a. V( X! I, @) e/ u3 m  c1 ?& Jr e t u r n ; }

+ R/ [$ V8 n; J! Wif (d2 &lt;= d3) {a = X[l+1];
0 K$ S( r' {7 ~9 ~+ r" J, J
7 h; e& q1 E" P( u- q: T( Vb = X[r];

d = d2;}
/ \; P. ~1 B% f  N
; w. i* _7 Z6 nelse {a = X[l];
2 c5 S0 M  }( ~* l6 n
b = X[r];

d = d3;}

4 h7 d0 H, k1 T% w. \! nr e t u r n ; }

/ /多于三个点，划分为两部分

int m = (l+r)/2; // X[l:m] 在A中，余下的在B中
; a" S3 I3 D1 d' p, o9 e
// 在Z[l:m] 和Z [ m + 1 : r ]中创建按y排序的表

( w8 G" C( t# T! r+ pint f = l, // Z[l:m]的游标

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

for (int i = l; i &lt;= r; i++)
+ y& B3 V- ^/ g: @0 r
4 s: n# Z0 h# V, ?if (Y.p &gt; m) Z[g++] = Y;

else Z[f++] = Y;

// 对以上两个部分进行求解

/ x7 Z0 E8 }: Rc l o s e ( X , Z , Y, l , m , a , b , d ) ;
  g% V6 B8 Y; `# u# B' h, o! u
/ {0 F' C& R# h9 w) Q- i, b7 ~float dr;

. L* S- t5 m- z$ n7 QPoint1 ar, br;
# Q4 [/ Z" V3 o# k
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;

% ?9 C$ |! V. ab = br;

7 F! S. B/ Y& @d = dr;}
$ K9 _! E9 x% ^. p7 r& q* K9 R
M e r g e ( Z , Y,l,m,r);// 重构Y

/ o) [6 t9 p  b/ y7 [/ /距离小于d的点放入Z
. p) V6 n& o* Z# U0 |! W+ m8 G0 ~% U
8 |* H" o+ i( S  K; j  Cint k = l; // Z的游标
- @3 X0 `/ G3 H- r( ?, b5 x
) ]! y0 `9 d4 |4 y6 r5 x6 sfor (i = l; i &lt;= r; i++)
8 S0 A8 k7 E# v& D+ b5 d3 v  [
if (fabs(Y[m].x - Y.x) &lt; d) Z[k++] = Y;
6 Z2 j: v( Y. w1 U4 K
// 通过检查Z [ l : k - 1 ]中的所有点对，寻找较近的点对
! {& K" b$ T! l( W/ f
for (i = l; i &lt; k; i++){
' x4 {9 i+ X) F, T" J( K# V
4 @3 y' R6 f! o6 Q, x9 ]. H* `9 Mfor (int j = i+1; j &lt; k &amp;&amp; Z[j].y - Z.y &lt; d;

5 O% P5 b0 u0 E2 y5 U/ bj + + ) {
# Q1 r/ d& f3 P% s
float dp = dist(Z, Z[j]);

8 y  A* D: l! H: A  tif (dp &lt; d) {// 较近的点对

. b6 N3 o! b' f6 f" qd = dp;

/ y& U5 b' l5 l0 |4 ^- E6 B2 \- ?a = X[Z.p];

b = X[Z[j].p];}
& C. H! ^. s9 |- g( `. |) d( W
}
" v; p! E0 H% k; ]9 D5 Q( l
1 B; H( n6 T% ]}
, s; E2 ?: Q. U4 _+ m
}
, o7 y' {5 M: x. j, f) q& e, @8 @
函数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 ]。
6 {& f. @: Z3 ~5 ?9 r; z+ c' a
  H3 A7 [% f" {4 @为实现图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

+ H( \5 |8 M0 e; I; r还是在RB中）与Z[j] 配对来实现，其中i＜j 且Z [ j ] . y - Z [ i ] . y＜ 。对每一个Z [ i ]，在2 × 区域内所检查的点如图1 4 - 1 7所示。由于在每个2 × 子区域内的点至少相距。因此每一个子区域中的点数不会超过四个，所以与Z [ i ]配对的点Z [ j ]最多有七个。
7 |* H# W; s! c7 c( z
! i, c+ t# _  N4 ]2. 复杂性分析

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

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