距离最近的点对

韩冰

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

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

6 y& |5 ~1 }! D0 e# h3 S& d, t通过检查所有的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进行递归求解。

+ B5 y! w! n3 r4 f- K( h9 I
if (n较小) {用直接法寻找最近点对
8 d" ~3 I3 x. y$ ?! I5 _0 o+ T
" o) Z0 Q4 O+ n/ I6 @& FR e t u r n ; }

% g0 w: V! H4 l/ ]' B- A& X; {  V: P+ y// n较大

+ d9 U: |+ r) w4 ]将点集分成大致相等的两个部分A和B

确定A和B中的最近点对

7 v  I- A6 O& `6 F$ ?确定一点在A中、另一点在B中的最近点对

从上面得到的三对点中，找出距离最小的一对点
5 P# w: U! _1 {
/ C3 f* F( J" t- D2 P/ {1 j图14-13 寻找最近的点对

- s4 q( ^. ]9 i# a" F& C
, k7 M0 l* ^  c( i8 {为了确定第三种情况下的最近点对，需要采用一种不同的方法。这种方法取决于点集是如何被划分成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。
& S' m0 h/ U( e  {
设是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）。
( p3 i; B6 j  {' P  u! d; m: f( Q
: `+ o3 m/ m" v7 L% b; D! A8 S例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所示。
6 t+ x$ C9 [% h2 j! w& B, @
1. 选择数据结构

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

  T, @6 V* d; o, `3 q8 \% e- f每个点可有三个参数：标号， x 坐标，y 坐标。假设标号为整数，每个点可用P o i n t l类（见程序1 4 - 8）来表示。为了便于按x 坐标对各个点排序，可重载操作符＜=。归并排序程序如1 4 -3所示。

! Z; h3 {. Z4 ?. e. r程序14-8 点类
, U7 Q" T  \4 o
, a% c+ m; H. Pclass Point1 {

! p  N8 `2 ~( i# A6 T! ?0 W1 ufriend float dist(const Point1&amp;, const Point1&amp;);

8 _6 W/ T! I5 ~) ifriend void close(Point1 *, Point2 *, Point2 *, int, int, Point1&amp;, Point1&amp;, float&amp;);
% N* v, S/ v) `8 P
friend bool closest(Point1 *, int, Point1&amp;, Point1&amp;,float&amp;);

friend void main();
3 x  n, l( V9 ]
  t9 S, ]) V7 N$ l+ {! Np u b l i c :
0 l+ o2 u: q7 K2 ?. `' U7 a- Y
% ]- w6 \- \. X. V! K4 A& Cint operator&lt;=(Point1 a) const
+ u, F, U6 J8 q, L$ G( o. G
7 }/ G4 L; I: L0 ]{return (x &lt;= a.x);}

p r i v a t e :

int ID; // 点的编号

float x, y; // 点坐标

} ;
& x* R: R0 z. m3 l. t9 r
class Point2 {

friend float dist(const Point2&amp;, const Point2&amp;);

8 M4 W, A5 r$ y. W# D3 @friend void close(Point1 *, Point2 *, Point2 *, int, int, Point1&amp;, Point1&amp;, float&amp;);
. u5 o! |: H9 X  k
, u5 K0 L0 y% g/ y0 J: H! Yfriend bool closest(Point1 *, int, Point1&amp;, Point1&amp;, float&amp;);

/ O, {5 b1 S$ mfriend void main();
' J9 l# |$ A8 [, X
p u b l i c :
8 D1 m2 d7 [: J7 F6 c$ E
int operator&lt;=(Point2 a) const
0 S) o& L, `1 `) a
{return (y &lt;= a.y);}

. H) T7 ?8 P1 h8 v" D" h4 np r i v a t e :
1 d, `" p- z& ~) Y( e) ]8 R/ K. ^
: O1 O2 F$ z* g4 b1 Lint p; // 数组X中相同点的索引
* G9 f5 _( x0 L# O
float x, y; // 点坐标

} ;
4 g- |/ A  i' v( 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中的对应点。
. k+ s9 Y) Z% W. v. d& d
/ r' _' r& Y% `( u* S' C" Z确定了必要的数据结构之后，再来看看所要产生的代码。首先定义一个模板函数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类的友元。
# E) ~1 M& L% s/ [9 o( w4 ~
7 r* }' v  W5 [9 L9 F程序14-9 计算两点距离

/ S4 N6 B1 h6 z' C3 ~template<CLASS T>

1 T$ @. h7 w' V& q9 O2 ]inline float dist(const T&amp; u, const T&amp; v)

{ / /计算点u 和v之间的距离

float dx = u.x-v. x ;

float dy = u.y-v. y ;
0 V, ~* H. B4 }, f2 w( I+ s: h
return sqrt(dx * dx + dy * dy);

8 J- p/ M1 a: g+ w}

+ C& \* q  o4 W  z  |9 |( ?0 Q9 B; o如果点的数目少于两个，则函数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 )，该函数实际求解最近点对。
  k+ {6 j3 Q- \
程序14-10 预处理及调用c l o s e

1 I* p9 b( }' f8 s5 W3 i6 y) A0 `3 wbool closest(Point1 X[], int n, Point1&amp; a, Point1&amp; b, float&amp; d)
# ?5 d) A4 o" G" `% _  Z
{// 在n &gt;= 2 个点中寻找最近点对

, R) v$ B0 s& P+ V9 i// 如果少于2个点，则返回f a l s e
1 Z1 @% R+ w4 G& O" O, S3 N' F9 [3 A
// 否则，在a 和b中返回距离最近的两个点

if (n &lt; 2) return false;
  p3 |- X& T& [- W
8 Q4 u8 Z& k' Z' y# p$ {// 按x坐标排序
4 A! r2 |2 W$ P+ I4 B0 I6 K; t6 ]
M e r g e S o r t ( X , n ) ;
$ I2 |2 h7 s0 c$ i2 p6 @4 R/ A  u
// 创建一个按y坐标排序的点数组
) V1 j6 t8 Y. e" Q6 w) c0 z. G: x4 q
Point2 *Y = new Point2 [n];
. r6 V4 w3 D' K
5 p+ R* s* d( Y; R$ tfor (int i = 0; i &lt; n; i++) {
  J# E: l" d/ ~+ E( ]7 L/ X
5 V4 V3 _, r  G: j# o. K  e. b5 Z# m// 将点i 从X 复制到Y

4 C! g& O! ?8 k0 Y8 zY.p = i;

+ M. g" `( |# A+ l1 m/ LY.x = X.x;

Y.y = X.y;
9 G' L( r/ D0 \; k, z. _+ l# q$ ~0 s7 s
1 f, H4 Q7 ^) D& ?}

+ @' A* _$ p) y2 j# a  A3 I3 iM e r g e S o r t ( Y,n); // 按y坐标排序

5 W/ D1 u, N- w9 w+ U6 C// 创建临时数组

Point2 *Z = new Point2 [n];
3 W) A) s9 U8 T+ h1 x! }
* L: b& C& @7 g% j7 e/ T// 寻找最近点对
: H5 B& @- @% q$ N
3 @" B: Z1 s# Q  Dc l o s e ( X , Y, Z , 0 , n - 1 , a , b , d ) ;

  l' `1 Y& N7 l6 q- P// 删除数组并返回
& n; T* D& t$ @$ R" B2 D
6 f5 M6 R/ T* h8 Q, E5 H5 F! Sdelete [] Y;

delete [] Z;
% l0 i" x- T0 U2 J9 b6 s
2 J& Z& o* E7 Q. B* u* zreturn true;

}
% ]% }3 W! X: p
程序1 4 - 11 计算最近点对
/ `7 o6 L# [* v0 v
; M+ X! m  P3 A$ K5 Jvoid close(Point1 X[], Point2 Y[], Point2 Z[], int l, int r, Point1&amp; a, Point1&amp; b, float&amp; d)

; L# S, c- [6 n, L$ N, ?2 i3 s{//X[l:r] 按x坐标排序
) B, Z' H9 k3 q! j& p
4 i3 P; r0 H- S+ r6 ~//Y[l:r] 按y坐标排序
4 s5 K4 r+ J: O" O
if (r-l == 1) {// 两个点

a = X[l];

b = X[r];

+ C8 ^$ E5 j& v1 a+ Id = dist(X[l], X[r]);

r e t u r n ; }
% [6 }/ E. d' V  q9 w1 L) N
if (r-l == 2) {// 三个点
" @2 D$ H, L5 h3 [; _" x, J& L
. x1 f* v- V5 C* Q/ h$ C' S% p// 计算所有点对之间的距离

% O. ?- p! w! Sfloat d1 = dist(X[l], X[l+1]);

float d2 = dist(X[l+1], X[r]);
$ n# s3 b6 ?2 R* T+ g( _
. `; f. w! v8 p8 g; S* J2 Zfloat d3 = dist(X[l], X[r]);

// 寻找最近点对

3 H) L+ ^( }( u5 bif (d1 &lt;= d2 &amp;&amp; d1 &lt;= d3) {
, V6 k& ^+ r. c9 F; t" t) O
/ b# I- C' I5 z7 d/ C- Ia = X[l];
9 a+ ?# ], t  r3 D0 z: b
b = X[l+1];

: O  W+ r6 w1 D  @- zd = d1;
2 G& t% A/ E1 n. w
r e t u r n ; }

if (d2 &lt;= d3) {a = X[l+1];
& c+ A' `; N5 Y; f* `
b = X[r];

d = d2;}
: V; T7 U$ S, o) ]% R8 a& y
2 f  N% N0 R2 B" Y- I' c; t+ b0 E! `else {a = X[l];
5 |/ {- k! g+ R$ Y' m9 N0 V6 L
/ W0 J+ j1 k, Y" I, q2 H: I$ ab = X[r];

d = d3;}

# _% K4 J+ ?& }  _( kr e t u r n ; }
/ _) ^2 ?2 K1 K% C  N! @, \
4 s/ P$ k; p2 O/ /多于三个点，划分为两部分
; C: Y$ `+ }! C3 a  G
int m = (l+r)/2; // X[l:m] 在A中，余下的在B中

2 |8 H( ^- }# t6 Y, {// 在Z[l:m] 和Z [ m + 1 : r ]中创建按y排序的表
# c/ B; e  E7 I, O* d
& U; Z. S+ @- k, wint f = l, // Z[l:m]的游标
$ G* |+ S6 i; H3 q/ x  e7 w- C3 s
4 D9 m9 V( z4 T/ y% rg = m+1; // Z[m+1:r]的游标

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

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

& A# `6 f9 c# X/ o8 ~4 `" melse Z[f++] = Y;
- [9 S; d: U6 l0 Y" G1 f( [0 c. `
// 对以上两个部分进行求解
5 b# A% z2 X7 M5 F
$ c! z2 T0 V5 p9 S1 Zc l o s e ( X , Z , Y, l , m , a , b , d ) ;
8 h( [+ W( \, M; P1 t6 q! Y8 y
float dr;

4 Q: W/ ]6 D+ yPoint1 ar, br;

5 ~/ W7 C, v) W4 ?7 Q; [% Kc l o s e ( X , Z , Y, m + 1 , r, a r, b r, d r ) ;

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

) G, [: F4 w" k7 kif (dr &lt; d) {a = ar;

8 T3 \$ w9 [$ d* V# o- tb = br;

6 y0 R+ m! W9 F' I1 id = dr;}

M e r g e ( Z , Y,l,m,r);// 重构Y
/ a" P- o- ?/ Q5 t! O4 e
/ /距离小于d的点放入Z
  `( T" }6 [9 h, K3 D$ x
* D, [9 t/ s5 K) H% gint k = l; // Z的游标

for (i = l; i &lt;= r; i++)
+ u- J, h- O4 O
if (fabs(Y[m].x - Y.x) &lt; d) Z[k++] = Y;
% M$ ~# M9 ]1 R' {
// 通过检查Z [ l : k - 1 ]中的所有点对，寻找较近的点对
/ _  b. n- V3 Q. l
% @# i& z" S4 S% U6 Q% Afor (i = l; i &lt; k; i++){
* x1 I! o3 n4 X, k# T& K/ W
for (int j = i+1; j &lt; k &amp;&amp; Z[j].y - Z.y &lt; d;

j + + ) {
( }5 ^+ ~1 I9 e9 l* \3 u
) k, {( Z8 f$ X, Lfloat dp = dist(Z, Z[j]);
2 O9 h0 F$ X6 F/ E) O5 S9 Y
3 c' }1 _) x) j0 iif (dp &lt; d) {// 较近的点对

5 ~# h9 w3 V7 ]& u8 A. d# O+ l: ]d = dp;
) ?! Y2 Z4 k7 N. e  u, B# r5 b
$ {, b/ w# q0 v$ \4 Ga = X[Z.p];

1 H' s' E+ V* n+ @5 Nb = X[Z[j].p];}
+ j( s8 b9 ^- ~1 M0 F, z
3 f6 V5 Q+ P6 Y}

# Z( E9 y) U7 u5 X" r8 X}
; t1 B. ~/ p! u+ z, z' M' ]8 h0 M
}

: f& C2 F# V% n$ K2 g8 D; U函数c l o s e（见程序1 4 - 11）用来确定X[1:r] 中的最近点对。假定这些点按x 坐标排序。在Y [ 1 : r ]中对这些点按y 坐标排序。Z[ 1 : r ]用来存放中间结果。找到最近点对以后，将在a, b中返回最近点对，在d 中返回距离，数组Y被恢复为输入状态。函数并未修改数组X。
' h. D* e1 o5 R  L  {
3 c- s1 i. L6 v) Q) U首先考察“小问题”，即少于四个点的点集。因为分割过程不会产生少于两点的数组，因此只需要处理两点和三点的情形。对于这两种情形，可以尝试所有的可能性。当点数超过三个时，通过计算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 |$ f! _1 n' h! t7 d0 H
$ t& q, V+ t* g7 o- s为实现图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

4 }/ ]5 ~# p! N& j" L还是在RB中）与Z[j] 配对来实现，其中i＜j 且Z [ j ] . y - Z [ i ] . y＜ 。对每一个Z [ i ]，在2 × 区域内所检查的点如图1 4 - 1 7所示。由于在每个2 × 子区域内的点至少相距。因此每一个子区域中的点数不会超过四个，所以与Z [ i ]配对的点Z [ j ]最多有七个。
3 T/ w; Z! s! Y- k( U
2. 复杂性分析
( ^1 u5 T- U2 O  K$ @3 o
* k8 A7 e. `- r2 X0 Q; f6 J令t (n) 代表处理n 个点时，函数close 所需要的时间。当n＜4时，t (n) 等于某个常数d。当n≥4时，需花费(n) 时间来完成以下工作：将点集划分为两个部分，两次递归调用后重构Y，淘汰距分割线很远的点，寻找更好的第三类点对。两次递归调用需分别耗时t (「n /2ù」和t (?n /2?).
8 l7 a9 a& `! E/ j
这个递归式与归并排序的递归式完全一样，其结果为t (n) = (nl o gn)。另外，函数c l o s e s t还需耗时(nl o gn)来完成如下额外工作：对X进行排序，创建Y和Z，对Y进行排序。因此分而治之最近点对求解算法的时间复杂性为(nl o gn)。</P>