距离最近的点对

韩冰

<>给定n 个点（xi，yi）（1≤i≤n），要求找出其中距离最近的两个点。
: {9 o8 l4 W9 O- ?! ?) M( ^
; W: M0 M- q- g2 O2 b4 ^例14-7 假设在一片金属上钻n 个大小一样的洞，如果洞太近，金属可能会断。若知道任意两个洞的最小距离，可估计金属断裂的概率。这种最小距离问题实际上也就是距离最近的点对问题。

2 L. K5 R+ m9 z$ x8 x0 t7 t- 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进行递归求解。

/ X8 d% E4 D2 T/ g! v( w
if (n较小) {用直接法寻找最近点对

R e t u r n ; }
9 \) A# L& I3 Q! V! O; x
2 {8 ?  {, q! n. h( G// n较大

0 p2 d" |& k7 @8 s, C将点集分成大致相等的两个部分A和B
/ P% v& F4 l4 q1 S) ^% H) e5 W
确定A和B中的最近点对

) Q% m7 U- A+ }) G确定一点在A中、另一点在B中的最近点对

- U/ r' u% ~! h* n: o! ]从上面得到的三对点中，找出距离最小的一对点

8 F8 c* E! R: j+ K" E图14-13 寻找最近的点对

  R8 f; G1 Y4 b4 v3 W1 F( _. s
9 W: R* H# p3 t3 A7 v8 X5 ]3 f) `) Z/ Q为了确定第三种情况下的最近点对，需要采用一种不同的方法。这种方法取决于点集是如何被划分成A、B的。一个合理的划分方法是从xi（中间值）处划一条垂线，线左边的点属于A，线右边的点属于B。位于垂线上的点可在A和B之间分配，以便满足A、B的大小。

& {, M; q6 D+ S& {例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。
3 A! @* A$ p5 p2 }" I
设是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) 是最近的点对。
) g: p8 c8 m4 q5 h9 D
为了确定一个距离更小的第三类点，RA 中的每个点最多只需和RB 中的6个点比较，如图1 4 - 1 6所示。
+ `+ r% R7 r( P2 U3 a$ P: h
! c' y2 p7 ?% G$ S! ~- {9 T1. 选择数据结构
# n. {6 c0 Z" c% [
为了实现图1 4 - 1 3的分而治之算法，需要确定什么是“小问题”以及如何表示点。由于集合中少于两点时不存在最近点对，因此必须保证分解过程不会产生少于两点的点集。如果将少于四点的点集做为“小问题”，就可以避免产生少于两点的点集。
& J/ {# i% n0 c# H
. a8 O' R1 d4 C, z8 ?每个点可有三个参数：标号， x 坐标，y 坐标。假设标号为整数，每个点可用P o i n t l类（见程序1 4 - 8）来表示。为了便于按x 坐标对各个点排序，可重载操作符＜=。归并排序程序如1 4 -3所示。
3 ]4 S6 D* L% U+ ]
程序14-8 点类
; [* z1 m% Z- t4 B; r) `( h
class Point1 {
2 y/ R: ?: F4 ?9 n. q
3 n+ c; w" E9 e4 L' B4 Cfriend float dist(const Point1&amp;, const Point1&amp;);

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

friend bool closest(Point1 *, int, Point1&amp;, Point1&amp;,float&amp;);
1 z# ?& \$ P( S
friend void main();
' K" `5 [  d+ e1 b# E" X  A8 y( q
p u b l i c :
9 }$ {  {* E5 R- ?, {% W9 ^
int operator&lt;=(Point1 a) const
+ r/ E$ L  N3 `# S' _. M9 o" _1 Q
& U9 s# V( w  a0 m# R$ D{return (x &lt;= a.x);}
" r0 Q- p7 H0 C  ~* l
p r i v a t e :
, F) q/ u5 ]9 W3 X  K) r* P
int ID; // 点的编号
' m  w7 g: R; P& }# O" e" I) t
5 t/ ?' l4 b: |( a" zfloat x, y; // 点坐标
0 @. r; `, h# B0 u
- w# h+ V- [3 Q. e0 n9 d; Q} ;
, k" s' `) m  ?- S: s/ T
class Point2 {
6 n6 t1 d5 R" V( [
0 a1 Q4 x% F# e% [; T5 ]  L, |& ~friend float dist(const Point2&amp;, const Point2&amp;);
0 [6 r7 n; Q/ M+ a9 E
friend void close(Point1 *, Point2 *, Point2 *, int, int, Point1&amp;, Point1&amp;, float&amp;);
) S: ]! p+ w- i
friend bool closest(Point1 *, int, Point1&amp;, Point1&amp;, float&amp;);

friend void main();

p u b l i c :

int operator&lt;=(Point2 a) const
/ I" ^4 h; C% Z) ~# ~. F7 \
{return (y &lt;= a.y);}

p r i v a t e :
0 t0 ~& d# U! @& ~+ g: s# r
+ r2 B+ j3 f5 ?" _  A) yint p; // 数组X中相同点的索引

float x, y; // 点坐标
' b0 O6 @* d- l7 T7 g
9 l  U/ X: i8 V( k' i1 p} ;

所输入的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 R: p$ w- Q% n& |  {1 K
. P+ i& b. o, U. `6 a( @+ M/ x程序14-9 计算两点距离

) r+ o( I7 U) a% c6 i: Ztemplate<CLASS T>

inline float dist(const T&amp; u, const T&amp; v)
2 \3 q9 f. W) g& S
* u1 M2 v' y8 R, `" _{ / /计算点u 和v之间的距离
# e) u8 w! }+ V' R" a6 Y
float dx = u.x-v. x ;

' h; A2 T' i, W* f( Sfloat dy = u.y-v. y ;
$ {6 }' W1 o) T9 A. X0 |
4 @7 Z! X$ U. q! `6 {! sreturn sqrt(dx * dx + dy * dy);

( A+ ~) ]( k0 B# z! G}

8 z9 P9 V9 v# \; j) n; L如果点的数目少于两个，则函数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 )，该函数实际求解最近点对。
$ F0 V1 D  W2 a1 m
程序14-10 预处理及调用c l o s e
( H0 l/ F' g! h( J' `$ L
bool closest(Point1 X[], int n, Point1&amp; a, Point1&amp; b, float&amp; d)

( c# ?, F, M+ R6 j: }4 S: p{// 在n &gt;= 2 个点中寻找最近点对

" e- {3 N# I7 U3 u+ l3 r9 T+ c// 如果少于2个点，则返回f a l s e

// 否则，在a 和b中返回距离最近的两个点

: B* T5 B: K# r  `7 X7 `if (n &lt; 2) return false;
) \) z' l/ t: [& B+ B) E# c
// 按x坐标排序
! e) j- k% s4 A2 F( U9 j4 ?
M e r g e S o r t ( X , n ) ;

// 创建一个按y坐标排序的点数组
5 C/ R; p; t) @" t
Point2 *Y = new Point2 [n];
/ s7 P' z, [9 t/ }
+ z/ b# \& `" T- H# O8 T) ?" p0 @for (int i = 0; i &lt; n; i++) {
7 E. v9 ]' G9 [& G
// 将点i 从X 复制到Y
( \; O: X$ T3 d: q; O0 n
Y.p = i;

Y.x = X.x;
) g5 a4 [% S# f+ ~3 N$ }( \5 U
' k: K; q% |# hY.y = X.y;

}

M e r g e S o r t ( Y,n); // 按y坐标排序
) L7 Q0 l0 O; l$ ?
7 h! N$ j; q7 U6 i// 创建临时数组

& r; v/ x1 h0 }) Y& R4 ]- h) KPoint2 *Z = new Point2 [n];

& |$ T, _. I+ s4 S// 寻找最近点对
9 H% T+ d4 L; v
" ?% b) |0 p" R) U8 y2 Lc l o s e ( X , Y, Z , 0 , n - 1 , a , b , d ) ;

1 i% A- a' @3 b  X) b// 删除数组并返回
! p! C3 p. X8 E% m: D2 q" ]
1 l6 j! X* r/ G7 k* ?8 a! Sdelete [] Y;

+ Y0 B- @" r: o$ Bdelete [] Z;
6 s$ d8 i2 T4 Z4 {# h8 i5 V+ [6 }% n
return true;
# \, u1 X$ C1 a, F$ ]
}
: X5 w# \5 h8 N1 [
程序1 4 - 11 计算最近点对
4 s5 M& A! u; [3 G$ h* A% t
  m" V' ]/ L' T' j4 zvoid close(Point1 X[], Point2 Y[], Point2 Z[], int l, int r, Point1&amp; a, Point1&amp; b, float&amp; d)

- Y  ^3 U0 Z; U7 C{//X[l:r] 按x坐标排序

//Y[l:r] 按y坐标排序
# L/ Y8 u; W& O8 c6 o8 r, b
if (r-l == 1) {// 两个点

a = X[l];

7 j% j( E& [$ D1 P; db = X[r];

& x9 `. }  e4 U; \' Ad = dist(X[l], X[r]);

1 `' X  L; X  c- ~1 }' s0 wr e t u r n ; }
  g, d) Q# H8 \$ W) ]
# q9 ^9 X, A1 @% Lif (r-l == 2) {// 三个点

! H9 q% K* z7 Q7 o2 d$ J// 计算所有点对之间的距离
' ?5 x2 T$ L5 _, F
float d1 = dist(X[l], X[l+1]);
- m0 J) }& x0 j# q% k% u) M9 M  j5 n
float d2 = dist(X[l+1], X[r]);
$ v8 G: A6 s$ F
$ P! g2 b& @/ M0 [float d3 = dist(X[l], X[r]);
# b+ Y) V8 E9 x
// 寻找最近点对

if (d1 &lt;= d2 &amp;&amp; d1 &lt;= d3) {

7 j7 l% \6 J, i' f) ^$ a) Pa = X[l];
9 z4 C9 C, n7 z) ~
# N# {) o' g( k# I$ j( Xb = X[l+1];
& `( p" y: H0 ]
d = d1;

: K* K1 s6 R4 w) br e t u r n ; }

. x( [7 z6 w! R! d2 h$ l( pif (d2 &lt;= d3) {a = X[l+1];
8 o) v( z$ c: k3 ~' V& P/ }$ O/ o
b = X[r];

d = d2;}

5 E# f  w+ ^) c2 L. ~else {a = X[l];
) R6 H) b3 ]6 Y, i" h- Z8 s1 J- \
7 [4 u9 ^; r* W; |b = X[r];

$ {" m( `3 r3 }$ U( M* n* K3 ?d = d3;}

2 h2 e+ e6 V# i* Vr e t u r n ; }

6 Z; [& f  V3 A/ /多于三个点，划分为两部分

int m = (l+r)/2; // X[l:m] 在A中，余下的在B中
; w0 g8 O8 g% s# D" W
. Y" n1 C* V7 ]8 [- @# f7 a// 在Z[l:m] 和Z [ m + 1 : r ]中创建按y排序的表

8 |% i: Q$ p& Y" c* wint f = l, // Z[l:m]的游标

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

# K8 ~$ r9 x, D# E, Q$ tfor (int i = l; i &lt;= r; i++)

if (Y.p &gt; m) Z[g++] = Y;
  E: l5 e# w7 W% Q& ^& r* G
else Z[f++] = Y;

' H0 c' o8 x4 ^0 U' d$ o// 对以上两个部分进行求解
5 m: Q8 O- {' a9 |; S# ]
+ V3 n2 n, Y% }7 V8 t# Q9 Rc l o s e ( X , Z , Y, l , m , a , b , d ) ;

% Z5 K$ p( m% a" T1 ?  S" Efloat dr;
: @" e2 R1 ~/ E% }9 U% f) b) c4 ^
Point1 ar, br;

c l o s e ( X , Z , Y, m + 1 , r, a r, b r, d r ) ;
/ F( |8 ^$ S$ k8 `* T
$ s, m6 M! f# V3 k0 p1 {7 Z/ C& p// (a,b) 是两者中较近的点对

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

b = br;

1 w0 w  b5 C) }" Vd = dr;}

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

: J$ M6 S3 Y6 B8 T/ /距离小于d的点放入Z
* |& }1 P3 F9 F/ x6 w
int k = l; // Z的游标

for (i = l; i &lt;= r; i++)

if (fabs(Y[m].x - Y.x) &lt; d) Z[k++] = Y;

// 通过检查Z [ l : k - 1 ]中的所有点对，寻找较近的点对
: I; b' d! O0 {, _
for (i = l; i &lt; k; i++){
% f% x* p# A" S5 H" |0 W1 F8 W" ]6 P( u
for (int j = i+1; j &lt; k &amp;&amp; Z[j].y - Z.y &lt; d;

9 u; ?& F- J4 k0 ?3 a4 }j + + ) {

float dp = dist(Z, Z[j]);

if (dp &lt; d) {// 较近的点对

d = dp;
9 z3 }; d1 X* R' Z% L+ F
a = X[Z.p];
/ r3 a& G1 S4 q# U$ H
" D- `8 D. \- h$ j+ U4 @b = X[Z[j].p];}

}

}
* \; Z8 \! W& j0 k
3 g) Y8 @. Y; B1 W}
( U, Z3 r! V6 d+ I
函数c l o s e（见程序1 4 - 11）用来确定X[1:r] 中的最近点对。假定这些点按x 坐标排序。在Y [ 1 : r ]中对这些点按y 坐标排序。Z[ 1 : r ]用来存放中间结果。找到最近点对以后，将在a, b中返回最近点对，在d 中返回距离，数组Y被恢复为输入状态。函数并未修改数组X。
" {% n$ B2 N* o- o
首先考察“小问题”，即少于四个点的点集。因为分割过程不会产生少于两点的数组，因此只需要处理两点和三点的情形。对于这两种情形，可以尝试所有的可能性。当点数超过三个时，通过计算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 ]。

为实现图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
" z! _; E, T( D# H0 [/ M7 H
5 Z. q5 Y' d' W) E还是在RB中）与Z[j] 配对来实现，其中i＜j 且Z [ j ] . y - Z [ i ] . y＜ 。对每一个Z [ i ]，在2 × 区域内所检查的点如图1 4 - 1 7所示。由于在每个2 × 子区域内的点至少相距。因此每一个子区域中的点数不会超过四个，所以与Z [ i ]配对的点Z [ j ]最多有七个。
$ A" L; w: L4 K* l, C, p
2 S/ h, L: n7 {% x1 b! \5 s0 p2. 复杂性分析
6 r9 L; H# D: l8 Q" m( U% D; G" s
. _# a. a" M+ d% s令t (n) 代表处理n 个点时，函数close 所需要的时间。当n＜4时，t (n) 等于某个常数d。当n≥4时，需花费(n) 时间来完成以下工作：将点集划分为两个部分，两次递归调用后重构Y，淘汰距分割线很远的点，寻找更好的第三类点对。两次递归调用需分别耗时t (「n /2ù」和t (?n /2?).
/ O! M) y6 i& R8 y9 k" X, Z
这个递归式与归并排序的递归式完全一样，其结果为t (n) = (nl o gn)。另外，函数c l o s e s t还需耗时(nl o gn)来完成如下额外工作：对X进行排序，创建Y和Z，对Y进行排序。因此分而治之最近点对求解算法的时间复杂性为(nl o gn)。</P>