距离最近的点对

韩冰

<>给定n 个点（xi，yi）（1≤i≤n），要求找出其中距离最近的两个点。
1 M9 |* O  G+ \/ E- ~1 ~0 B% U  M
例14-7 假设在一片金属上钻n 个大小一样的洞，如果洞太近，金属可能会断。若知道任意两个洞的最小距离，可估计金属断裂的概率。这种最小距离问题实际上也就是距离最近的点对问题。

, L5 j0 l+ z+ W通过检查所有的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 i. p) g3 s: S0 C6 i0 W" |+ Z
if (n较小) {用直接法寻找最近点对

( j2 d6 \' q  U  w: pR e t u r n ; }
/ N" K, _+ C9 Y& j8 e6 M2 ?3 H9 t/ Z
// n较大
- O: H: N# n' g1 k0 }" C
- t  i0 U4 I+ s4 w将点集分成大致相等的两个部分A和B
+ N2 H+ s6 o7 ~. K/ @
确定A和B中的最近点对

确定一点在A中、另一点在B中的最近点对
& z# g. M$ Z) F$ j0 [+ j
' T6 B( T' d0 R& e; }, Z从上面得到的三对点中，找出距离最小的一对点
) d* i( O- _( C! l/ G3 B
0 G4 V- P4 L+ x' B图14-13 寻找最近的点对
$ w) E( D2 O1 K7 D' y8 Z1 M0 B5 ]

为了确定第三种情况下的最近点对，需要采用一种不同的方法。这种方法取决于点集是如何被划分成A、B的。一个合理的划分方法是从xi（中间值）处划一条垂线，线左边的点属于A，线右边的点属于B。位于垂线上的点可在A和B之间分配，以便满足A、B的大小。
/ d( K, u4 Z, O) g+ }3 V9 H, O. Y
例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。
, p# }4 a( P  \: r6 }4 [
设是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）。
) Y8 A; j6 j3 ]: I9 v  m
4 b% q- }4 j2 L例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所示。

5 Q8 n1 ~! Q, n% `! w1. 选择数据结构
$ J7 F2 Q/ ~* R
为了实现图1 4 - 1 3的分而治之算法，需要确定什么是“小问题”以及如何表示点。由于集合中少于两点时不存在最近点对，因此必须保证分解过程不会产生少于两点的点集。如果将少于四点的点集做为“小问题”，就可以避免产生少于两点的点集。

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

程序14-8 点类

( ~/ o6 ?  z  v% R, _7 ?$ U% oclass Point1 {

! z' }. ~+ ~; kfriend float dist(const Point1&amp;, const Point1&amp;);
+ S1 v5 M  H& z6 v
+ v$ j" I" G4 qfriend void close(Point1 *, Point2 *, Point2 *, int, int, Point1&amp;, Point1&amp;, float&amp;);

% g2 Y  B) H+ k5 `friend bool closest(Point1 *, int, Point1&amp;, Point1&amp;,float&amp;);

friend void main();

/ f- Y3 [0 }. N/ X1 Pp u b l i c :

int operator&lt;=(Point1 a) const
" e% J# N/ P  B7 S8 |
{return (x &lt;= a.x);}

% `; M2 `$ a( W! `p r i v a t e :
9 @4 Q* }8 N4 {3 w
6 V0 S4 F$ N  \6 h1 ?int ID; // 点的编号
% a7 s0 y# b6 u. B6 y! t) t
float x, y; // 点坐标
9 U2 x7 L6 U; |9 s" c# N
} ;

class Point2 {

6 u0 b0 {, y7 v0 ?- Efriend float dist(const Point2&amp;, const Point2&amp;);
. L# Y5 b. [  Z5 N
0 [% s1 q% \4 Z; _# }: efriend void close(Point1 *, Point2 *, Point2 *, int, int, Point1&amp;, Point1&amp;, float&amp;);

+ Z% \+ |3 _* d9 t: S: vfriend bool closest(Point1 *, int, Point1&amp;, Point1&amp;, float&amp;);
: x& @/ ^) v& p
: q7 `1 y5 C1 Lfriend void main();

1 q* ^2 L- c8 p" d% Bp u b l i c :

int operator&lt;=(Point2 a) const

8 T2 b$ d9 r0 `( Z6 l{return (y &lt;= a.y);}

" \8 {+ M0 ?1 v8 r- j- H! Hp r i v a t e :
: p; _' ^4 w5 P' e& G
6 R0 G2 |" _( w. o$ L6 E5 H, Bint p; // 数组X中相同点的索引

7 d3 R/ M8 |5 N# E( U7 efloat x, y; // 点坐标

# E0 N' m7 r3 N( f" s* ^# \1 C* S} ;
: U# o+ w# d# P
4 Q4 M) |5 Z- u! f所输入的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类的友元。

9 ]3 M6 x% c# R9 T, w2 y程序14-9 计算两点距离
; m3 R1 D: g1 Z6 ]0 t( B
template<CLASS T>

inline float dist(const T&amp; u, const T&amp; v)
; G" G, I9 f; ?8 k5 x9 M$ O
4 z+ A4 F7 d$ `5 ?( p{ / /计算点u 和v之间的距离

float dx = u.x-v. x ;
9 m8 I+ x9 `. `2 J% n3 C& @* c
float dy = u.y-v. y ;
2 ?7 P) J4 A: a% N+ v
- [, m: w4 |6 p" g6 Q+ Xreturn sqrt(dx * dx + dy * dy);

}
: j& q6 S$ r2 G5 o9 _8 B; Z( D
- a6 T  z  v2 V! ?如果点的数目少于两个，则函数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 )，该函数实际求解最近点对。
3 Q; o8 L/ N3 {) ]" j6 @
+ z3 J) x9 z/ B3 o% K4 R! E7 K程序14-10 预处理及调用c l o s e

, W! }" r& b$ H& |bool closest(Point1 X[], int n, Point1&amp; a, Point1&amp; b, float&amp; d)
/ Z1 y5 p. I+ H7 z% P  g, Z
, F/ y% |! p5 H$ ?& M{// 在n &gt;= 2 个点中寻找最近点对
9 q8 w. S, X2 v
// 如果少于2个点，则返回f a l s e
  A, b+ ~3 f0 w# P
// 否则，在a 和b中返回距离最近的两个点
# J0 c& g) m; V: \% u
. h" H$ E1 T3 {6 vif (n &lt; 2) return false;
! E7 u6 [/ a5 A& y$ [6 L6 D* \- X
// 按x坐标排序
8 ?5 S2 K; a, E' Q: W
M e r g e S o r t ( X , n ) ;

// 创建一个按y坐标排序的点数组

# a& O6 V9 s, f0 X4 M9 i# c+ QPoint2 *Y = new Point2 [n];
5 t6 X  Y! `. A% V
for (int i = 0; i &lt; n; i++) {

8 v- I/ w. R8 [7 D" a// 将点i 从X 复制到Y
( c7 I3 o' H+ H2 n
( c2 [" W: G3 Q: O8 ]. KY.p = i;

Y.x = X.x;

Y.y = X.y;
- ~9 c3 l% _' P+ J
}

M e r g e S o r t ( Y,n); // 按y坐标排序
, y3 u- X5 p/ w/ R
( Z3 t; d" |4 u" e  w// 创建临时数组
2 E1 y7 a7 |* T7 A
5 [& W0 _6 I( W' D( P# |Point2 *Z = new Point2 [n];

// 寻找最近点对

( z( Y* b! N# c4 ?- Nc l o s e ( X , Y, Z , 0 , n - 1 , a , b , d ) ;
  d, L1 ^! B' }$ a
% `1 j# j, \6 N* {) O! L// 删除数组并返回

/ j8 d* ?2 Y0 @+ n0 odelete [] Y;

delete [] Z;
9 v, W8 B: T. C+ Y
return true;
: n; t+ \, _0 v7 C4 s, L
7 l- i; E2 t* Y+ c' F8 L' O}

4 Z" S1 a! P+ [- e- e程序1 4 - 11 计算最近点对

void close(Point1 X[], Point2 Y[], Point2 Z[], int l, int r, Point1&amp; a, Point1&amp; b, float&amp; d)

4 x9 ?; E' A0 y4 h. ?9 e{//X[l:r] 按x坐标排序
! X) h% H' B6 w. d! ?
//Y[l:r] 按y坐标排序

2 T5 K# ]# n5 Lif (r-l == 1) {// 两个点

a = X[l];

b = X[r];

d = dist(X[l], X[r]);
, J, b. o- P% [; @4 J# I
( h. s' I) t* |4 r+ X2 Mr e t u r n ; }
( X/ m- z3 B, p/ |
9 v" ~8 j' q2 l0 r: rif (r-l == 2) {// 三个点

// 计算所有点对之间的距离

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

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

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

// 寻找最近点对

if (d1 &lt;= d2 &amp;&amp; d1 &lt;= d3) {
, R  `! A( l/ K9 |: }' s* i
. B" h. ~- d1 @: Q3 [# Wa = X[l];

$ K7 J9 \% L( f: T7 z' H' w+ ub = X[l+1];
9 C5 s) c4 s1 \, ^0 V
d = d1;
1 l5 U- g8 K9 g) I1 ?
* h% P$ g. H; E* Q5 b' fr e t u r n ; }
! w# \$ Q( o7 ?
/ R# L, S3 c' R! Q# s; ]3 j+ x8 Yif (d2 &lt;= d3) {a = X[l+1];

b = X[r];

( |% }1 j4 l7 J4 }9 j; Fd = d2;}
, U) u3 H( h# x1 y
else {a = X[l];
! ^' c) M- P- J
b = X[r];
: r- o$ c2 x& I3 T# t
d = d3;}

4 h, X4 T* |0 Z3 y/ k! Nr e t u r n ; }
8 g2 h, }, P% b# g+ l  E
& O" E$ ]: t0 D8 }" ], x/ /多于三个点，划分为两部分
6 k- J  M( v& C, g
: \! i! q5 m. H. N- T2 t. oint m = (l+r)/2; // X[l:m] 在A中，余下的在B中
; f' Q% S' D0 K7 w/ t
// 在Z[l:m] 和Z [ m + 1 : r ]中创建按y排序的表
4 h7 \/ X, i; L1 }8 b* ^
: U) w/ z; P5 I6 iint f = l, // Z[l:m]的游标

g = m+1; // Z[m+1:r]的游标
4 [) m8 S6 K  ~' i
for (int i = l; i &lt;= r; i++)
- Z  T; Z1 h8 O! s, ]
2 a8 f+ G1 A, ~$ B- o# Q/ M0 E& Iif (Y.p &gt; m) Z[g++] = Y;

9 L: E+ d) c* c0 t4 @# _else Z[f++] = Y;

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

c l o s e ( X , Z , Y, l , m , a , b , d ) ;

& V& R) \- b( z- H  B9 j- g" Yfloat dr;
# e" M% _0 d2 D
Point1 ar, br;
0 X: W9 @5 o, w
6 z% b6 C  o2 w3 A2 w. gc l o s e ( X , Z , Y, m + 1 , r, a r, b r, d r ) ;

// (a,b) 是两者中较近的点对
0 ?$ P+ h8 w, L7 B& V" n: U( F8 S% }
) o0 F6 W- o" O8 Y, o+ n1 G. Xif (dr &lt; d) {a = ar;
" p* |8 s+ m9 K' m
3 \, x! |3 N" J1 {b = br;
6 u" [( g6 Y& L$ T0 {
d = dr;}

M e r g e ( Z , Y,l,m,r);// 重构Y
9 @6 m7 v$ |4 A7 G4 {( ^: G5 Q
* b/ @2 m6 d2 T0 M! Z/ /距离小于d的点放入Z

' W! O$ G3 N7 l) b: h; `$ Yint k = l; // Z的游标
: {, U9 Q! K8 L1 D+ I: H
for (i = l; i &lt;= r; i++)

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

$ D2 P  U* a/ g! n/ x: |// 通过检查Z [ l : k - 1 ]中的所有点对，寻找较近的点对

& Z+ p) ~. q6 h/ E% @* ?for (i = l; i &lt; k; i++){
: o3 W- \: O; C. b3 s8 ^/ j: P
; l6 f6 S% _, f4 U' o1 ffor (int j = i+1; j &lt; k &amp;&amp; Z[j].y - Z.y &lt; d;
5 F$ T: S9 Q5 l- ?- u
j + + ) {
- F" |7 Z# o! v! H# O$ G
float dp = dist(Z, Z[j]);

if (dp &lt; d) {// 较近的点对
+ e/ K2 k3 C2 U& ]: l9 g7 L- C. H% P
; @' ]% W. P0 n9 Jd = dp;
+ ]& A2 `! X4 \& W1 c9 ]6 g0 N+ o
% _- X5 |2 G" S6 }! u# M! Z6 [a = X[Z.p];
; o6 V" N6 z) |' H3 f
b = X[Z[j].p];}
0 `/ J. z$ V' D& T. N
+ P8 _0 O- [  [: m9 Q}

}

+ r4 F. u6 m6 x( O6 \$ Y& _! r}
& I8 h8 _  D0 f5 _2 {( b8 x" O
. p- \) D  i4 E0 P* O函数c l o s e（见程序1 4 - 11）用来确定X[1:r] 中的最近点对。假定这些点按x 坐标排序。在Y [ 1 : r ]中对这些点按y 坐标排序。Z[ 1 : r ]用来存放中间结果。找到最近点对以后，将在a, b中返回最近点对，在d 中返回距离，数组Y被恢复为输入状态。函数并未修改数组X。
3 T) |! T- A: q6 d
+ i) K; d4 R  K; a$ V首先考察“小问题”，即少于四个点的点集。因为分割过程不会产生少于两点的数组，因此只需要处理两点和三点的情形。对于这两种情形，可以尝试所有的可能性。当点数超过三个时，通过计算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 ]。
0 k0 G" j6 v6 s1 a) M+ b$ s* @
4 s5 m) G4 t0 r+ ]4 O0 l为实现图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 ]最多有七个。

, d" R1 f& T! p+ r. E) @, z2. 复杂性分析
8 d+ u5 }* q# V' s' B
令t (n) 代表处理n 个点时，函数close 所需要的时间。当n＜4时，t (n) 等于某个常数d。当n≥4时，需花费(n) 时间来完成以下工作：将点集划分为两个部分，两次递归调用后重构Y，淘汰距分割线很远的点，寻找更好的第三类点对。两次递归调用需分别耗时t (「n /2ù」和t (?n /2?).

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