距离最近的点对

韩冰

<>给定n 个点（xi，yi）（1≤i≤n），要求找出其中距离最近的两个点。
' L# s! P3 \0 Y1 O4 j/ |" z, J( A) w
9 o! S/ d% X* A7 _" M- ~7 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进行递归求解。

  t; [3 j, V/ a& h% D' ~
! v- a; P1 H* r4 ?8 j4 Wif (n较小) {用直接法寻找最近点对

$ ^* S0 k, q8 X& Y( u$ I% p2 p. B& QR e t u r n ; }
* Q7 r! w- H8 X, c! \/ ~' A: O$ ~
// n较大

: Z3 E* }5 I, \. b3 y将点集分成大致相等的两个部分A和B

% k/ z6 K3 R$ E1 _  \5 i确定A和B中的最近点对

确定一点在A中、另一点在B中的最近点对
( r$ l5 Z& ?" n7 R' j. J
从上面得到的三对点中，找出距离最小的一对点

1 ~# O. ^# Y2 f图14-13 寻找最近的点对
/ p; `$ C# Y. J2 D: n  w3 [. @" c
9 ]8 m% n0 o2 u' q$ K8 Q3 y/ i
' [  t; m0 B0 l: a9 d为了确定第三种情况下的最近点对，需要采用一种不同的方法。这种方法取决于点集是如何被划分成A、B的。一个合理的划分方法是从xi（中间值）处划一条垂线，线左边的点属于A，线右边的点属于B。位于垂线上的点可在A和B之间分配，以便满足A、B的大小。
) G. o0 i/ B8 @1 i
3 E4 g" _7 T( ?( v6 ?例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。

设是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）。
6 ^6 V4 w: ]/ x/ I7 \2 q
( Y* r  R6 i2 R9 k例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) 是最近的点对。
, j1 W/ q0 l0 g" R, s' r
为了确定一个距离更小的第三类点，RA 中的每个点最多只需和RB 中的6个点比较，如图1 4 - 1 6所示。

1. 选择数据结构
0 p; Y) {6 W4 A* S2 j1 r8 |  h
为了实现图1 4 - 1 3的分而治之算法，需要确定什么是“小问题”以及如何表示点。由于集合中少于两点时不存在最近点对，因此必须保证分解过程不会产生少于两点的点集。如果将少于四点的点集做为“小问题”，就可以避免产生少于两点的点集。

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

# [" l: l/ L! b) N4 _/ k. a程序14-8 点类

# d( K4 d3 ~- n1 H; D+ gclass Point1 {
, J4 |' B6 s- m! ^$ M# [
friend float dist(const Point1&amp;, const Point1&amp;);
* r- b" d$ o/ o* |8 l1 Z, ^
: x$ U1 o* H8 E) t) `, Wfriend void close(Point1 *, Point2 *, Point2 *, int, int, Point1&amp;, Point1&amp;, float&amp;);

6 _+ @5 U( ?; m8 V9 F8 Yfriend bool closest(Point1 *, int, Point1&amp;, Point1&amp;,float&amp;);

friend void main();

2 R6 x  j; F; X, x" B+ T3 yp u b l i c :
1 D  Q& i9 o% N
int operator&lt;=(Point1 a) const
% v9 z9 v% I! r( k3 D/ h9 P
{return (x &lt;= a.x);}
/ H0 y7 e. P2 o& j
p r i v a t e :

2 h6 V8 f( u: D2 C7 ^; nint ID; // 点的编号
' v% W  U( x$ \5 W$ r2 v& H" T# v
- V. f$ \1 p- Ofloat x, y; // 点坐标
0 Q" d. C9 g9 V; K3 \! j
& l( `! {( X. y/ W- i4 E5 V7 u} ;

class Point2 {
0 J( @/ D  m" f6 H( P* s! m
friend float dist(const Point2&amp;, const Point2&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;);

' S. |* W. S" f/ xfriend void main();
' Y+ q; c' m) T( |. C# ^
! e8 p5 T# Q8 n- sp u b l i c :

int operator&lt;=(Point2 a) const

' e2 ?7 B# @2 Q# A0 y% O+ b8 H{return (y &lt;= a.y);}
" t# ~" i$ |$ j$ w$ O# o; W
- f4 _. ^8 y) u) V' _3 fp r i v a t e :

int p; // 数组X中相同点的索引
5 Y! {$ G& I1 {; D, D
float x, y; // 点坐标
' P+ ~  Q- l2 x. V
7 C* k" o2 Q3 Q0 h& [1 H} ;
2 f, M  I! R  J  t# v
- ?$ Y/ H, |) _5 U所输入的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中的对应点。
* r; e7 a8 f: E
确定了必要的数据结构之后，再来看看所要产生的代码。首先定义一个模板函数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类的友元。
: ~0 s; [  M5 l5 ]; E
程序14-9 计算两点距离
  X& F( A( V  C
template<CLASS T>
; l! c5 ~$ k" _2 ^% Z9 n
6 w4 O9 D; P* w3 t* `inline float dist(const T&amp; u, const T&amp; v)

# ^+ ]: `8 d& o3 A7 ?% v. S{ / /计算点u 和v之间的距离
* |4 c$ j  M$ M' y; L) L, p
float dx = u.x-v. x ;
$ B4 }' r6 w' Q4 b; i5 m( B" s
" z, P1 |, ]; r: R& cfloat dy = u.y-v. y ;
6 i# b+ O7 L2 ?( N: [  y2 B) U
5 s8 X: O( q2 z( @return sqrt(dx * dx + dy * dy);

% |3 |+ k! {% \# B}
& b! [( n" q6 @" b
) ~3 j. z5 H2 I& x1 |9 z如果点的数目少于两个，则函数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 )，该函数实际求解最近点对。

+ ~0 E$ N& `  J/ B程序14-10 预处理及调用c l o s e

bool closest(Point1 X[], int n, Point1&amp; a, Point1&amp; b, float&amp; d)

+ x# W, d" Q: v{// 在n &gt;= 2 个点中寻找最近点对
6 H: [3 R9 r( Q3 X' M+ U
' {! r2 w( i1 v$ S) i3 s. E// 如果少于2个点，则返回f a l s e
9 b$ V4 G# x2 ?
// 否则，在a 和b中返回距离最近的两个点

if (n &lt; 2) return false;
& E7 W$ L) c  G7 @( _# m* p
// 按x坐标排序

3 @" v9 T* o  Y0 J" z5 }4 NM e r g e S o r t ( X , n ) ;

7 u! g; N9 ~3 U! ?: W6 Z// 创建一个按y坐标排序的点数组

6 @- h/ j$ o8 y* g! Y. c' b- `( ^Point2 *Y = new Point2 [n];

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

// 将点i 从X 复制到Y
/ p8 |6 C: L  V- e9 q
( S5 ]2 v( Q. {0 l7 YY.p = i;
- F2 q! Y+ K/ ^+ i# j+ [& |( Q6 {
* u% ^( B/ o: e- E& s; m6 j! ]Y.x = X.x;

9 `( R+ G8 E& |* C1 tY.y = X.y;

}

M e r g e S o r t ( Y,n); // 按y坐标排序
; L; u" M2 T" N6 P  c) V1 k! t
# B9 [  ^, l1 x// 创建临时数组

Point2 *Z = new Point2 [n];
) J4 T1 j8 u# P* I
- @1 Y! ]9 l5 H  V. v// 寻找最近点对
0 N+ p0 w/ B! T: g  }
c l o s e ( X , Y, Z , 0 , n - 1 , a , b , d ) ;
" J6 y, Q0 O8 e' i6 [4 X0 Y4 J$ o8 M7 k
& B  w5 q9 S  s) b4 c// 删除数组并返回
$ A) \! i, L. U# q7 ]
delete [] Y;

delete [] Z;
9 I+ Q; d4 A0 k. l. d
: j( C1 Q2 h. l8 I# y& p) u6 Nreturn true;
& v: b& g- I9 H$ N, [
}

4 n2 }( ^  M4 S; o, T程序1 4 - 11 计算最近点对

. s. y' X+ M3 z# n) i1 ^& ovoid close(Point1 X[], Point2 Y[], Point2 Z[], int l, int r, Point1&amp; a, Point1&amp; b, float&amp; d)

{//X[l:r] 按x坐标排序
4 ]0 b1 F  m/ y
" R' ^9 O2 S- e* i* E- ?9 v//Y[l:r] 按y坐标排序
7 C$ X# H0 C( j) _
if (r-l == 1) {// 两个点

a = X[l];
8 i: H' z; d1 `+ m( L% ^$ x
b = X[r];
( d$ ]) |" h0 j, B
d = dist(X[l], X[r]);
2 D2 @  u. T* P$ ^
  [1 P- y" {3 F# H) ?' M' Vr e t u r n ; }

if (r-l == 2) {// 三个点
5 z% N$ V& L7 N; d
' K" o6 T5 d8 ~// 计算所有点对之间的距离

: h' v. l+ j) R2 D" pfloat d1 = dist(X[l], X[l+1]);

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

! z( @3 l6 F5 S, C  W9 ~float d3 = dist(X[l], X[r]);

// 寻找最近点对
5 u' @' B* R7 W4 L
if (d1 &lt;= d2 &amp;&amp; d1 &lt;= d3) {
( I# i* Q. S. `9 ~$ v
a = X[l];
) _& a' c) E& `' P; }! e  y
- T: Z* l3 u4 p! {% o8 zb = X[l+1];

d = d1;
' Q- }5 p9 J8 r$ g; p4 {- ~; a: J
5 G) K* d3 z" O" s0 Wr e t u r n ; }
. u& g0 j4 D" v- B" W( h
9 ^% M9 {0 ]% e1 a* z* E. qif (d2 &lt;= d3) {a = X[l+1];
7 ?1 I* z  d; m; y4 ?. s1 @
8 b3 ~1 C4 j7 Wb = X[r];
9 O6 L+ Y0 [, [8 Y% I: p) _; i
d = d2;}
. c9 T/ S) B1 o
else {a = X[l];
1 L) i' V' L  Q! E2 g! K0 M! F
7 e  p7 `6 d; Y7 c" ub = X[r];
; q% Z. {( U) G- }7 d3 P
d = d3;}

1 D1 {8 X' ]. R2 V( r' lr e t u r n ; }

/ /多于三个点，划分为两部分
: h' i; I6 z; u6 P. W) \6 M1 V
3 x5 _- e+ U  vint m = (l+r)/2; // X[l:m] 在A中，余下的在B中
2 T+ H0 ?& A* ^! U
# i3 N$ p! m1 i// 在Z[l:m] 和Z [ m + 1 : r ]中创建按y排序的表

int f = l, // Z[l:m]的游标
" u- M& I3 [, h& H( m' J* `5 g$ b9 R
g = m+1; // Z[m+1:r]的游标
* C5 [# C4 ^! {/ C
for (int i = l; i &lt;= r; i++)
* h# T) d$ K, [# q5 H9 X
% c8 D7 f; S% r4 A4 i7 l/ uif (Y.p &gt; m) Z[g++] = Y;

else Z[f++] = Y;
( y+ H- V1 j, O4 E# ~* Z" v
7 q' L, g6 v0 Z1 K// 对以上两个部分进行求解

" u( {/ i- [& V" U+ b1 v- wc l o s e ( X , Z , Y, l , m , a , b , d ) ;

# ~, z  h4 h% e- afloat dr;

# K$ i: A) J$ J" U7 N4 X6 ePoint1 ar, br;

& s; v5 J7 H- f0 |: @! W: z% Hc l o s e ( X , Z , Y, m + 1 , r, a r, b r, d r ) ;
6 s8 G7 Z9 W$ z; i
% d, E5 K# K* a4 e# j- |// (a,b) 是两者中较近的点对

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

3 |4 A* I' l0 ?, Q+ {) vb = br;
$ R! N( Q5 B  t; I
d = dr;}
: L; |7 q& e& {% C, h; O3 W
M e r g e ( Z , Y,l,m,r);// 重构Y
9 m- Y# [0 r/ ^
/ /距离小于d的点放入Z
* V2 o" h+ Q4 y2 A  P6 u4 n3 S# o
int k = l; // Z的游标

$ I) b/ k; v# K# zfor (i = l; i &lt;= r; i++)

4 C+ B- F8 H# D* a- L0 M+ lif (fabs(Y[m].x - Y.x) &lt; d) Z[k++] = Y;

% [$ F7 B) E' i' h6 P7 _// 通过检查Z [ l : k - 1 ]中的所有点对，寻找较近的点对
& P4 }' z5 ~' A, c5 e
# v4 x* a7 {. m$ B# U7 jfor (i = l; i &lt; k; i++){
6 K% ~# q; f: S: D4 |
for (int j = i+1; j &lt; k &amp;&amp; Z[j].y - Z.y &lt; d;

j + + ) {

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

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

4 Q7 `4 }+ P* b) Nd = dp;
5 D0 O3 ~0 e: u  K& {
a = X[Z.p];
5 w2 E+ U+ M0 _3 H/ B
b = X[Z[j].p];}
5 \. x4 b0 _- N. e
; m4 I' ^) v3 N; d5 O! }+ p6 n) ~}

- ^( L! J) o: \+ Z}

7 V" c9 _" r" z8 V6 o% W: F: Q}
6 f2 D5 r& @  s  P" ^2 W
0 W7 J/ G/ D: b5 r0 Y9 d  p函数c l o s e（见程序1 4 - 11）用来确定X[1:r] 中的最近点对。假定这些点按x 坐标排序。在Y [ 1 : r ]中对这些点按y 坐标排序。Z[ 1 : r ]用来存放中间结果。找到最近点对以后，将在a, b中返回最近点对，在d 中返回距离，数组Y被恢复为输入状态。函数并未修改数组X。
: y3 v, i) f7 r; m' r# l' ]& X  W
首先考察“小问题”，即少于四个点的点集。因为分割过程不会产生少于两点的数组，因此只需要处理两点和三点的情形。对于这两种情形，可以尝试所有的可能性。当点数超过三个时，通过计算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 ]。

+ E. F' I. G, s2 q为实现图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 ]最多有七个。

; `1 {! A) l& V8 \2. 复杂性分析
! Q3 L$ v* L: _4 `
* `+ _- N/ z3 _- c7 N令t (n) 代表处理n 个点时，函数close 所需要的时间。当n＜4时，t (n) 等于某个常数d。当n≥4时，需花费(n) 时间来完成以下工作：将点集划分为两个部分，两次递归调用后重构Y，淘汰距分割线很远的点，寻找更好的第三类点对。两次递归调用需分别耗时t (「n /2ù」和t (?n /2?).
3 R; s' z8 C; t% j! C5 S# Y
' }- {  Z0 i' A: q7 W7 i这个递归式与归并排序的递归式完全一样，其结果为t (n) = (nl o gn)。另外，函数c l o s e s t还需耗时(nl o gn)来完成如下额外工作：对X进行排序，创建Y和Z，对Y进行排序。因此分而治之最近点对求解算法的时间复杂性为(nl o gn)。</P>