距离最近的点对

韩冰

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

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

通过检查所有的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进行递归求解。

" [1 O* Y- w; q
( z7 X9 ?9 W. d/ Sif (n较小) {用直接法寻找最近点对

0 \$ J1 R* x/ [; j1 h, ^. FR e t u r n ; }

' t2 L+ i* n! N( K  r/ g// n较大

将点集分成大致相等的两个部分A和B
1 M5 g: R* J0 p
0 a* F- m  G! ?$ J4 v) I确定A和B中的最近点对
$ |+ z4 |0 f# k, K
确定一点在A中、另一点在B中的最近点对
" G5 k' \( Q+ E4 h% U
从上面得到的三对点中，找出距离最小的一对点

图14-13 寻找最近的点对

9 }* [7 ]9 Q& z7 e, d% Q" _+ i
为了确定第三种情况下的最近点对，需要采用一种不同的方法。这种方法取决于点集是如何被划分成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。
" F( e/ M9 b: A; w( H! X
( u  L' y6 @, `设是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）。
5 X* x" s! r) a2 \' H
0 R7 l) U2 i' b/ E; C& Q* [, E6 k7 M例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$ Y, j# ~% {3 c
为了确定一个距离更小的第三类点，RA 中的每个点最多只需和RB 中的6个点比较，如图1 4 - 1 6所示。

7 c2 S8 l% g  D. y" Y! i! {* G1. 选择数据结构
, t- {' y5 S- T* M: `) l  E1 F
" z& O6 h2 P6 M% B2 d为了实现图1 4 - 1 3的分而治之算法，需要确定什么是“小问题”以及如何表示点。由于集合中少于两点时不存在最近点对，因此必须保证分解过程不会产生少于两点的点集。如果将少于四点的点集做为“小问题”，就可以避免产生少于两点的点集。

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

; V) j9 T4 b- T3 d# N2 {程序14-8 点类
; X6 A+ l) J+ w: L# a
! @8 t6 W7 U/ t& T. y2 b. l: H# Nclass Point1 {
: P  L) O, r& l" P( e% D  C9 f
friend float dist(const Point1&amp;, const Point1&amp;);
' h& \# G/ A' N* o
friend void close(Point1 *, Point2 *, Point2 *, int, int, Point1&amp;, Point1&amp;, float&amp;);
% O. i: B1 ]0 P- x- {* Q
" O) K9 u6 W1 e1 C/ p( Q  Rfriend bool closest(Point1 *, int, Point1&amp;, Point1&amp;,float&amp;);
' ]% f: O& y* ^# a4 F
: }! {2 C- W1 V. U; @3 I0 d8 Bfriend void main();

p u b l i c :
0 ^' i2 D) t! S4 q1 R0 W
int operator&lt;=(Point1 a) const

! P+ X  ~$ x; z6 Y4 I{return (x &lt;= a.x);}
9 |# R4 A0 H+ N3 w8 G
; c: m% g. d7 ?/ m  E& c8 ~p r i v a t e :

! V9 ~& v# \! W# k+ u# ?int ID; // 点的编号

float x, y; // 点坐标

} ;

class Point2 {
6 b- T/ j; G5 d4 N
, x" M9 [, b9 j6 J6 q% k& lfriend float dist(const Point2&amp;, const Point2&amp;);

1 m0 u+ y: R0 ^, b/ {friend void close(Point1 *, Point2 *, Point2 *, int, int, Point1&amp;, Point1&amp;, float&amp;);
0 @$ d3 t6 {. E; `2 i& Q# q3 u
& C6 n  `5 D  B9 A4 d7 p! Sfriend bool closest(Point1 *, int, Point1&amp;, Point1&amp;, float&amp;);

, y3 k3 V. h; |! s# v* Q* s7 ffriend void main();
: c2 V+ @3 [' \, j; @3 x$ y5 {
; I) s# Q$ S8 u7 @; f/ I3 T* cp u b l i c :

int operator&lt;=(Point2 a) const
4 g) Y# s9 s' j5 D5 ~+ A' e
% ], f3 i$ Y5 U' ]6 ~{return (y &lt;= a.y);}
4 j5 Y! r; r' ?' n
p r i v a t e :

9 k; y5 p5 o! A5 rint p; // 数组X中相同点的索引
' x, T+ F1 L( [. y& e/ L4 g8 V
2 |5 M, j+ \: o6 Rfloat x, y; // 点坐标
; D  d; j4 K  |
} ;

% M2 h0 m$ G9 G+ ~  g所输入的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中的对应点。

* H8 _' ^3 a! V3 f3 q确定了必要的数据结构之后，再来看看所要产生的代码。首先定义一个模板函数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 计算两点距离
, f( y! D+ t+ W3 K& b
2 {6 ~8 R9 c6 Stemplate<CLASS T>

# I0 \/ u3 G0 Finline float dist(const T&amp; u, const T&amp; v)

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

float dx = u.x-v. x ;
* N* l- d. c2 {
! h! o" r: j0 S& d$ hfloat dy = u.y-v. y ;

return sqrt(dx * dx + dy * dy);

* `1 G( x$ I( D2 H2 K}
6 s* w$ }. Z/ s
$ L& I6 ]0 f9 o" k3 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 )，该函数实际求解最近点对。

' Q+ t1 ?. j8 ], H程序14-10 预处理及调用c l o s e

8 f  u) @* b, P6 i+ a9 s+ Ubool closest(Point1 X[], int n, Point1&amp; a, Point1&amp; b, float&amp; d)
; R3 S  w5 R# Y7 j# ^
{// 在n &gt;= 2 个点中寻找最近点对

+ g* \5 N) E$ [2 d5 O// 如果少于2个点，则返回f a l s e
) x- v7 R/ ]; n* F& W0 D; s
; y, b0 l: i- N& I// 否则，在a 和b中返回距离最近的两个点
; O/ V! w' Y! K' M
: F: @  X- A7 ~9 ?: Lif (n &lt; 2) return false;

// 按x坐标排序
* g9 z" Z9 l! K' h
M e r g e S o r t ( X , n ) ;

: _" x5 C* f, z6 z% i  V% z: C! y+ e// 创建一个按y坐标排序的点数组
1 p5 {4 Z% ^* @- K, a' k* v
Point2 *Y = new Point2 [n];

5 c+ U7 \) [# j! Ffor (int i = 0; i &lt; n; i++) {
0 a8 Q6 ?4 Q5 r1 L. I+ N0 T, ?
1 W; b$ P3 z- V' S$ d: w// 将点i 从X 复制到Y
9 Q, N; ]! s' w2 \. o: `. j* k; O
' H4 N1 \3 W$ X( K' s0 V# RY.p = i;

Y.x = X.x;

5 j; P1 `1 \% zY.y = X.y;

}
  x, B9 a! k6 ?: c0 N& g! H9 b
2 z7 K& X& o: k* R/ r  n2 cM e r g e S o r t ( Y,n); // 按y坐标排序

// 创建临时数组
1 d$ i1 L5 a- X* w  @
9 ]% s- j5 Z' |. l& ~0 FPoint2 *Z = new Point2 [n];
- q8 a6 ]2 L1 }" ?
3 E3 m4 n* u' j* X4 ]// 寻找最近点对

c l o s e ( X , Y, Z , 0 , n - 1 , a , b , d ) ;
0 y+ x0 D$ ^- p. C& i3 i0 [
// 删除数组并返回
$ ?! d- W( b' F! q1 S' U  I' m
* G6 z5 m8 Y& O6 [delete [] Y;
: u. F9 J8 R: ~$ `% _/ y
9 ^1 b5 h! t3 pdelete [] Z;
5 T' d  R+ u7 e
- f( @6 S! X* [( j8 wreturn true;

) E* M6 ~. e5 ~}

程序1 4 - 11 计算最近点对
# j1 E& c8 ~0 T' Q( L$ A! ~
void close(Point1 X[], Point2 Y[], Point2 Z[], int l, int r, Point1&amp; a, Point1&amp; b, float&amp; d)

{//X[l:r] 按x坐标排序

2 l1 Z) g7 Z+ \/ M# ~//Y[l:r] 按y坐标排序
2 g- B4 @! v6 z* L6 r% G
  E7 f) b9 [' X4 k- A- c! b- g( Cif (r-l == 1) {// 两个点
2 J$ F5 S8 d. C4 w( z; g! ^4 j
a = X[l];
6 E  w) i3 p2 f% N* z3 h1 Z$ C
5 A3 v9 ~. `2 l2 M# Ab = X[r];
( Z4 g: y+ |+ }4 u
d = dist(X[l], X[r]);

r e t u r n ; }
* ^5 Q# Z! y' J: n5 S3 t
if (r-l == 2) {// 三个点

2 _" m9 T4 ]5 ]* [/ L  F/ U  c; B( j// 计算所有点对之间的距离
) h/ R. K6 X0 Y( ^1 |) b( ]
/ }0 \6 |! d; O6 Nfloat d1 = dist(X[l], X[l+1]);

3 y5 i& j3 v& S! f! ]  S4 W- W4 H5 pfloat d2 = dist(X[l+1], X[r]);
" S) ?( K4 W0 N6 C% Q
float d3 = dist(X[l], X[r]);
. }5 {# l# }! ?
// 寻找最近点对

9 ]( v5 ^% R) |/ ]8 dif (d1 &lt;= d2 &amp;&amp; d1 &lt;= d3) {

a = X[l];

# t9 q/ [) a3 l! K( gb = X[l+1];

d = d1;
( r7 X1 c- I. i0 u4 s
r e t u r n ; }
8 }: x7 ?$ o7 ]/ Q5 S
if (d2 &lt;= d3) {a = X[l+1];

b = X[r];
+ A* c  m3 S% d
' [5 V: z2 J: l: v& H* f& B, k+ D& \d = d2;}

else {a = X[l];
" J/ I$ z) x0 n+ g) p
b = X[r];
4 e! N, \. s, W8 ]7 p0 M6 i! B
d = d3;}

r e t u r n ; }

) E% a8 K" y! I: l+ c) X: r/ /多于三个点，划分为两部分

4 M! x6 g3 B  x9 D) q3 p  dint m = (l+r)/2; // X[l:m] 在A中，余下的在B中
# Q1 Q. P2 g' Q
// 在Z[l:m] 和Z [ m + 1 : r ]中创建按y排序的表
1 q6 V" {0 g2 W; z) M' M
int f = l, // Z[l:m]的游标

0 O$ J; ]  h8 Q* i, s; ug = m+1; // Z[m+1:r]的游标

" J! T0 X: I- s& Z( f) Gfor (int i = l; i &lt;= r; i++)
6 H1 L: I8 l. k. \
% Y2 k! {1 {. Q3 Jif (Y.p &gt; m) Z[g++] = Y;

else Z[f++] = Y;

, E% M# V- t2 \+ ]3 P2 p0 U// 对以上两个部分进行求解

c l o s e ( X , Z , Y, l , m , a , b , d ) ;
# S6 O8 C0 i  s: m0 m
( g% @# }' l' afloat dr;
9 D& w0 C! t# @' @: l
4 t% c( {2 D* z+ _. s6 mPoint1 ar, br;
0 o- X7 H% R. S0 d
c l o s e ( X , Z , Y, m + 1 , r, a r, b r, d r ) ;
9 X& ~: ]) @7 s+ [8 g$ p: q7 U% {0 u' Y
& k7 {4 M. h' g- i* T* h// (a,b) 是两者中较近的点对

1 U1 g7 x1 @7 D' `& Y, Bif (dr &lt; d) {a = ar;

- d9 V3 e  _1 J- ^$ Kb = br;
6 t# c$ _  c# _$ p/ J
d = dr;}

M e r g e ( Z , Y,l,m,r);// 重构Y
5 M& v) W8 L- o! ^( }# t! ?
  h" J( K  b& |' q/ /距离小于d的点放入Z

int k = l; // Z的游标

! z' ^* o8 I" m5 L1 y1 Y) [5 p$ |for (i = l; i &lt;= r; i++)
5 @0 ]. v" y+ z4 K+ k, t! `
& e, \8 V7 g( M7 G* r9 |if (fabs(Y[m].x - Y.x) &lt; d) Z[k++] = Y;
+ M3 c) T6 t% k& N( Y3 f0 N
, h# \! l/ L; Q( I// 通过检查Z [ l : k - 1 ]中的所有点对，寻找较近的点对

for (i = l; i &lt; k; i++){
8 j( k2 I& S  \0 @
for (int j = i+1; j &lt; k &amp;&amp; Z[j].y - Z.y &lt; d;

9 q; z' ]( a8 m& R" R  {j + + ) {

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

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

d = dp;

5 K- n1 [6 J0 p4 A- C  P2 Aa = X[Z.p];

8 d4 U% Y) J  p; ~) Yb = X[Z[j].p];}

}

4 \0 K3 {. J* p8 v2 \, f5 s}
9 j+ j, r6 o- A" B( C
}

函数c l o s e（见程序1 4 - 11）用来确定X[1:r] 中的最近点对。假定这些点按x 坐标排序。在Y [ 1 : r ]中对这些点按y 坐标排序。Z[ 1 : r ]用来存放中间结果。找到最近点对以后，将在a, b中返回最近点对，在d 中返回距离，数组Y被恢复为输入状态。函数并未修改数组X。
, U8 a8 n0 ?7 x0 T2 v/ I
6 k3 K  ?% s  s/ Q. u+ i2 ^首先考察“小问题”，即少于四个点的点集。因为分割过程不会产生少于两点的数组，因此只需要处理两点和三点的情形。对于这两种情形，可以尝试所有的可能性。当点数超过三个时，通过计算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 ]。
  x% P$ g* T$ O9 Z6 O) }
6 k: I+ q0 d  k7 |# t+ G& t# 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

还是在RB中）与Z[j] 配对来实现，其中i＜j 且Z [ j ] . y - Z [ i ] . y＜ 。对每一个Z [ i ]，在2 × 区域内所检查的点如图1 4 - 1 7所示。由于在每个2 × 子区域内的点至少相距。因此每一个子区域中的点数不会超过四个，所以与Z [ i ]配对的点Z [ j ]最多有七个。

4 Q" b$ `) ?/ u( l; ~2. 复杂性分析

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

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