距离最近的点对

韩冰

<>给定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进行递归求解。

, Q8 h* O& @. p& V5 X- r
# Z% ^7 z4 s/ \if (n较小) {用直接法寻找最近点对
& B! y  D8 z6 w9 u  a% b# G+ |# a
8 b. y4 x/ x% l) e3 iR e t u r n ; }
9 H! b$ o7 q8 p! m
// n较大
4 f) W. L( s. }% ?5 s
; Z8 X' a6 P' \. A) W& d6 ^将点集分成大致相等的两个部分A和B

# @" t! f( f  a4 h5 K3 k2 W确定A和B中的最近点对
( P  V6 o1 `7 ]. U
确定一点在A中、另一点在B中的最近点对
0 t: |5 L* M5 ]
2 p, L- c: ^( X8 r# P- W从上面得到的三对点中，找出距离最小的一对点
5 z! g, ]2 i' m0 s( @$ J; m; {
. \! L6 _0 ?( p( T图14-13 寻找最近的点对
6 r" r9 {' S' u9 \
6 ?+ r; ]( W, n: x: L
为了确定第三种情况下的最近点对，需要采用一种不同的方法。这种方法取决于点集是如何被划分成A、B的。一个合理的划分方法是从xi（中间值）处划一条垂线，线左边的点属于A，线右边的点属于B。位于垂线上的点可在A和B之间分配，以便满足A、B的大小。
7 V& V4 O; v: j# Z' X7 U! 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。
4 ?: r6 Q9 A# Z
设是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）。

* f& {7 Z' V5 o* T/ u& `例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) 是最近的点对。

! L! v2 s' ~5 Y/ c1 A为了确定一个距离更小的第三类点，RA 中的每个点最多只需和RB 中的6个点比较，如图1 4 - 1 6所示。
3 [3 @$ y8 J7 X- N
- u# O3 _: |" U2 X1. 选择数据结构

2 {# [" `+ W  `4 t! _' |为了实现图1 4 - 1 3的分而治之算法，需要确定什么是“小问题”以及如何表示点。由于集合中少于两点时不存在最近点对，因此必须保证分解过程不会产生少于两点的点集。如果将少于四点的点集做为“小问题”，就可以避免产生少于两点的点集。
% l; e: V9 W5 T, G+ ~! g, O$ f
每个点可有三个参数：标号， x 坐标，y 坐标。假设标号为整数，每个点可用P o i n t l类（见程序1 4 - 8）来表示。为了便于按x 坐标对各个点排序，可重载操作符＜=。归并排序程序如1 4 -3所示。

0 h/ [0 ^1 N7 @" }3 Z程序14-8 点类
+ i2 A* T& T' o1 b# ^
. n! R* [" Y( R7 M4 wclass Point1 {

friend float dist(const Point1&amp;, const Point1&amp;);
( o/ t, F. O) {/ h8 i, L8 P
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;);
9 _9 M( v, j$ F6 f; N
7 d( `# B! E8 m) S1 |+ X' G' ?friend void main();
# Q. M6 p- m) ^- h
/ j! K+ L& j1 C4 Kp u b l i c :

/ ^9 o: c: |( `& J# Tint operator&lt;=(Point1 a) const

# ~4 {8 W, d. b, \+ s! X{return (x &lt;= a.x);}

p r i v a t e :
. G- N1 b# A. Y4 Y1 L, P
int ID; // 点的编号

float x, y; // 点坐标

2 u7 N, h% R1 k4 F: \" U* E1 o3 N} ;

! J& e2 p- T6 G7 }4 M' Q) xclass Point2 {
  M/ ]( Q; r. _5 b
. F: o2 ^: y  @  Z! q3 T& l0 Nfriend float dist(const Point2&amp;, const Point2&amp;);

friend void close(Point1 *, Point2 *, Point2 *, int, int, Point1&amp;, Point1&amp;, float&amp;);
3 J1 [0 S7 ^7 A4 S! C9 U
# {# m6 M8 k1 k& n6 Sfriend 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

8 J2 f4 O7 C9 F7 `  S8 }7 R{return (y &lt;= a.y);}
- J: a+ c$ g7 t: V' Y
) i! h6 L) a  E5 x& Y: X! ~; xp r i v a t e :
) v- y% G; [+ a; T, {* |& l
int p; // 数组X中相同点的索引
; X8 I( e7 Q9 P* N2 X. f0 K
9 w+ s9 a8 N' Pfloat x, y; // 点坐标
$ Z1 O' _: K9 q% Y- W$ E: Z9 K* X
5 K) a- q( K- P/ o' Q4 X$ s9 V} ;

+ r$ O! @4 b2 b# Y. i所输入的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中的对应点。

+ g, h: ]' `1 s) u确定了必要的数据结构之后，再来看看所要产生的代码。首先定义一个模板函数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 计算两点距离

" A8 i$ \9 |4 z  G' Y1 rtemplate<CLASS T>
( y& p5 e- f0 H( x$ c
inline float dist(const T&amp; u, const T&amp; v)

{ / /计算点u 和v之间的距离
0 |( W; e0 F" ~& G3 ~
) [6 ]! F3 x8 H1 h3 _# K0 ?( ]! W/ ]float dx = u.x-v. x ;

float dy = u.y-v. y ;
& ?$ w% Q6 }  t7 [% }' c6 [
5 q' q5 n3 A) \: @/ o) z8 V1 T, l1 T! Qreturn sqrt(dx * dx + dy * dy);

}
' x" r  S3 Q- p/ u4 z0 N9 B3 x
/ m7 j4 @$ s, [9 Y% C5 t如果点的数目少于两个，则函数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 )，该函数实际求解最近点对。
( G+ X+ c& p. x6 J
3 k2 K1 J# ]+ {0 ?/ u: s& ~; @6 y# d  ^程序14-10 预处理及调用c l o s e
  p, C" ]+ K8 i: _1 p) O
bool closest(Point1 X[], int n, Point1&amp; a, Point1&amp; b, float&amp; d)
$ p1 i' @6 H% x. v+ G" T4 R
! u1 k. H, w. g8 i4 O2 o7 v' h{// 在n &gt;= 2 个点中寻找最近点对
% y3 D& L; f8 f4 F6 a9 |
! q. H+ W. R; y. {// 如果少于2个点，则返回f a l s e
1 J) ^% u$ N& U$ l! v: U
// 否则，在a 和b中返回距离最近的两个点

+ x9 t1 `: t8 c8 Cif (n &lt; 2) return false;
0 A! S/ P, |- c- v0 W
/ J. R. \# V  }# W3 ~// 按x坐标排序
* M7 w/ T' m- E
# K1 d8 b$ w1 |; P5 ]5 PM e r g e S o r t ( X , n ) ;

" j6 k" Q0 s" t0 {" J5 a9 C, P// 创建一个按y坐标排序的点数组

Point2 *Y = new Point2 [n];

for (int i = 0; i &lt; n; i++) {

9 ?$ q  \1 A1 D' t4 r* t// 将点i 从X 复制到Y
! Z& }/ ^* v" e- _- J
Y.p = i;
% i# G9 @! {" C  ^
Y.x = X.x;
9 ~2 M7 K/ T; k3 T1 B( A: h( ]4 E
Y.y = X.y;
% [2 X1 U- A! S8 K( q( h
0 b- \3 K  U* X5 `1 l3 j; K}
; m9 u3 c, ]6 W/ F
, K4 Q0 u! C) r' r+ _# E. O1 gM e r g e S o r t ( Y,n); // 按y坐标排序

// 创建临时数组

Point2 *Z = new Point2 [n];

// 寻找最近点对
3 C+ ]7 w: E$ {0 ^& e
c l o s e ( X , Y, Z , 0 , n - 1 , a , b , d ) ;

% `6 D8 S. n: S! P// 删除数组并返回

delete [] Y;
) F) d/ D0 [; g! c
delete [] Z;

# j/ f1 y# b$ m4 ]return true;
! n! X" Y2 I6 G( ^
}

+ W" [  `5 d) s$ q程序1 4 - 11 计算最近点对
1 k* z; f: ]) K3 k
+ W& s+ P( [3 @void close(Point1 X[], Point2 Y[], Point2 Z[], int l, int r, Point1&amp; a, Point1&amp; b, float&amp; d)

# V3 ?) j9 j/ {% \{//X[l:r] 按x坐标排序

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

/ q* Z1 j7 {: ?0 `9 }if (r-l == 1) {// 两个点

/ Y: q- M8 Q' ^: d% aa = X[l];
3 B# g8 g, p- V' c# R- V
$ v$ G- b# {8 l/ Hb = X[r];

d = dist(X[l], X[r]);
8 q2 D2 ^9 L* s8 n* x
r e t u r n ; }

if (r-l == 2) {// 三个点
; B* K  o$ B& }
6 J. F8 @$ X8 X" T  p$ n# Y// 计算所有点对之间的距离

; r( s" z* ]+ L5 z7 s% G& s+ b6 p, Yfloat d1 = dist(X[l], X[l+1]);

float d2 = dist(X[l+1], X[r]);

float d3 = dist(X[l], X[r]);

+ {# X; z9 z) r* K9 w// 寻找最近点对

if (d1 &lt;= d2 &amp;&amp; d1 &lt;= d3) {
) ]! F; r; v/ p  ?, A
  a0 M. d9 @- w8 t- z8 {a = X[l];
0 B5 E4 ^5 E7 K  x4 z& ?  o& ?8 G
b = X[l+1];

d = d1;
/ V# A% K2 G5 P  m' b, j3 i0 ], k! U
! j3 ^, L, l1 }1 B/ z- R2 wr e t u r n ; }

if (d2 &lt;= d3) {a = X[l+1];

; ^# u% }* i/ E# o5 h; P  k' r* nb = X[r];

d = d2;}
, }" t+ S2 W3 g2 y3 h% A8 d7 k1 B; v
else {a = X[l];

b = X[r];
# y7 R( L5 t& _6 `( \0 w+ M9 U
d = d3;}
& L- W# x+ z  E2 w  B
4 Y$ }$ {9 K3 x: _3 P. O% y. nr e t u r n ; }
  C. \) u  I& ]& o( z! q
/ /多于三个点，划分为两部分

' Y* Y$ d" h3 h: C, d& k) oint m = (l+r)/2; // X[l:m] 在A中，余下的在B中
  h/ r  @9 ?0 I( g; o
- C! n( Y  U6 I, O( m// 在Z[l:m] 和Z [ m + 1 : r ]中创建按y排序的表
0 U9 \: [/ `7 A/ |: L$ l. n- _
int f = l, // Z[l:m]的游标
8 x0 Z0 p, H; k# F7 s/ j4 i0 I9 i' e
g = m+1; // Z[m+1:r]的游标

* Y, y# x) B# U5 G: g4 _% Pfor (int i = l; i &lt;= r; i++)

if (Y.p &gt; m) Z[g++] = Y;
. Z& B: v& ]# n  P$ S4 g: A2 p
else Z[f++] = Y;

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

c l o s e ( X , Z , Y, l , m , a , b , d ) ;
  J: G! M$ G% x5 I8 \
float dr;
* [" T: p3 p8 Q/ P
" I. v9 Z( H# b4 k3 {Point1 ar, br;
# ]' Q; t6 Q3 W; E
( R& Z: p6 o3 @6 o4 L1 Zc l o s e ( X , Z , Y, m + 1 , r, a r, b r, d r ) ;

// (a,b) 是两者中较近的点对
% T, ~, d! c: C
if (dr &lt; d) {a = ar;

b = br;

d = dr;}
/ X9 \# y9 a; b$ L! c8 z$ l% g8 S
M e r g e ( Z , Y,l,m,r);// 重构Y
, e4 n% R6 F* {" x! n4 V. }6 P& d
/ /距离小于d的点放入Z
6 p% m( `3 U3 {! @) D8 Y' w
int k = l; // Z的游标
  O" s  O4 C0 U" [
for (i = l; i &lt;= r; i++)

if (fabs(Y[m].x - Y.x) &lt; d) Z[k++] = Y;
$ Y, W7 z! h8 O# F  F0 E$ U
// 通过检查Z [ l : k - 1 ]中的所有点对，寻找较近的点对
4 z, N; m3 Y4 {7 b) i, I
$ _( a7 J4 W9 N+ U4 p6 ?% s) F: M/ r) cfor (i = l; i &lt; k; i++){
5 Q; m+ V# [9 z* s
for (int j = i+1; j &lt; k &amp;&amp; Z[j].y - Z.y &lt; d;

! r% B* r; ]9 ]3 V: ^: vj + + ) {

float dp = dist(Z, Z[j]);
  Y2 s( p4 ]) Y- x3 ?$ R6 M
0 K& \  U5 y$ F0 }4 Lif (dp &lt; d) {// 较近的点对

d = dp;

( Q6 g, ~4 [4 _& a5 v- ]a = X[Z.p];
0 ~3 H  \) l9 P7 D) _
b = X[Z[j].p];}

}

}
- o. u, x4 [/ c4 u  T
8 X, Y$ T7 i: \) l}
: P+ n8 d4 F: t* p1 v' e
9 V; S9 }8 Z% P; ]+ [* k: r函数c l o s e（见程序1 4 - 11）用来确定X[1:r] 中的最近点对。假定这些点按x 坐标排序。在Y [ 1 : r ]中对这些点按y 坐标排序。Z[ 1 : r ]用来存放中间结果。找到最近点对以后，将在a, b中返回最近点对，在d 中返回距离，数组Y被恢复为输入状态。函数并未修改数组X。
$ [+ |$ L4 ], w/ a$ l: r
4 k' o) m2 |& v. L首先考察“小问题”，即少于四个点的点集。因为分割过程不会产生少于两点的数组，因此只需要处理两点和三点的情形。对于这两种情形，可以尝试所有的可能性。当点数超过三个时，通过计算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

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

/ V2 F) C2 |% [" y  x" n8 E2. 复杂性分析

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

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