距离最近的点对

韩冰

<>给定n 个点（xi，yi）（1≤i≤n），要求找出其中距离最近的两个点。
. E0 c" J3 w  Y6 L2 g
例14-7 假设在一片金属上钻n 个大小一样的洞，如果洞太近，金属可能会断。若知道任意两个洞的最小距离，可估计金属断裂的概率。这种最小距离问题实际上也就是距离最近的点对问题。

( ]* b( V5 D! f/ P: H: |: S) H通过检查所有的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进行递归求解。
9 @9 K5 d* S0 ?2 x# ~! M

if (n较小) {用直接法寻找最近点对
" Q8 c, Z" f9 x$ s$ N9 }/ Q+ n
R e t u r n ; }

( k* c" L' m+ D7 @$ p6 y// n较大
& n( ]' X  j! ?# O4 j4 h
将点集分成大致相等的两个部分A和B

$ D* p8 j. X! G- [# ?# @9 F确定A和B中的最近点对

确定一点在A中、另一点在B中的最近点对

从上面得到的三对点中，找出距离最小的一对点
4 t, h- v! t2 K  G7 k
图14-13 寻找最近的点对
0 V% o: h2 [6 A4 k, l- H- F
5 @5 F' X2 r( v: o
/ \% x5 X1 L, Q' k3 z% D$ p为了确定第三种情况下的最近点对，需要采用一种不同的方法。这种方法取决于点集是如何被划分成A、B的。一个合理的划分方法是从xi（中间值）处划一条垂线，线左边的点属于A，线右边的点属于B。位于垂线上的点可在A和B之间分配，以便满足A、B的大小。

4 T$ P) ?( ~7 m: @8 @例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。

+ J- t) \4 K% M: j, S$ h设是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）。

例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) 是最近的点对。
1 J) f7 {+ l4 s4 v  p/ p. r: ]
为了确定一个距离更小的第三类点，RA 中的每个点最多只需和RB 中的6个点比较，如图1 4 - 1 6所示。
! D6 \4 z) B- E: e; v4 D+ D
1. 选择数据结构
1 p$ b& O( c" Z
, h( P( C. e  s. z( C- w为了实现图1 4 - 1 3的分而治之算法，需要确定什么是“小问题”以及如何表示点。由于集合中少于两点时不存在最近点对，因此必须保证分解过程不会产生少于两点的点集。如果将少于四点的点集做为“小问题”，就可以避免产生少于两点的点集。
6 z$ Q5 ]- n3 R2 z
7 L( ]  d) C& b! K0 y4 J每个点可有三个参数：标号， x 坐标，y 坐标。假设标号为整数，每个点可用P o i n t l类（见程序1 4 - 8）来表示。为了便于按x 坐标对各个点排序，可重载操作符＜=。归并排序程序如1 4 -3所示。
4 j! A* t7 e  v: c
程序14-8 点类

0 v- Z1 J5 h' [, m4 V% hclass Point1 {
4 @" V* M, g7 s3 @% t
- {* g# O' }7 p7 _friend float dist(const Point1&amp;, const Point1&amp;);
) W, u- T4 x8 n% g
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;);

; ]  h. E- M0 {# Tfriend void main();
9 S6 x' c# d* P
p u b l i c :
8 i; V2 f3 T5 r4 g, w
int operator&lt;=(Point1 a) const

# B& M$ V; N; E6 I- R( D% \# q{return (x &lt;= a.x);}

/ ]$ |6 \# @2 l' W6 [% dp r i v a t e :

int ID; // 点的编号
/ R$ e1 }3 r! h; }; i! {8 s: J1 ^
# G6 R- g9 @- y& F+ g8 zfloat x, y; // 点坐标

: q, N3 X) e7 P. c* h  b! \} ;
1 d9 N7 m6 `# g8 a
2 g' P5 K; a7 U" s. rclass Point2 {
% J5 ?. M  A0 I
! o7 V. m6 i/ A' c' X% \* Yfriend float dist(const Point2&amp;, const Point2&amp;);

8 b  h; b" l0 j" U, V4 c# sfriend void close(Point1 *, Point2 *, Point2 *, int, int, Point1&amp;, Point1&amp;, float&amp;);

% Z4 X4 y" H+ D$ k  [: cfriend bool closest(Point1 *, int, Point1&amp;, Point1&amp;, float&amp;);
; @1 h: O+ g1 [$ I& p& W
friend void main();
4 H/ p8 r% N; b/ D. ~
p u b l i c :

int operator&lt;=(Point2 a) const
8 d1 e, e4 x( |9 R3 v$ z
' m! u+ L- v3 y& K/ q  j{return (y &lt;= a.y);}
; N5 S7 Q, ]8 J
p r i v a t e :
8 D5 L$ |/ l' r) G! U, ~
int p; // 数组X中相同点的索引
" M, Z% u/ t% g3 v7 w' ?# k% @
1 _. z' l" g) Z. T6 v' X9 T5 ufloat x, y; // 点坐标

* K7 W/ o& a0 M3 a} ;
7 f& A- U( C; Y  ^3 B
所输入的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类的友元。

程序14-9 计算两点距离

* Z) a& G  h' {4 l. Ztemplate<CLASS T>
* }4 Y  b! ]. s0 V4 s3 E
& ]/ o2 M" |) v+ G& j- }inline float dist(const T&amp; u, const T&amp; v)

2 j# P3 @$ o6 K8 Z; c{ / /计算点u 和v之间的距离

; Z5 i# B2 e  ?float dx = u.x-v. x ;
/ ]+ D- ~5 ?& a3 P% ]; w( i$ ?
8 @, r6 b: r1 Mfloat dy = u.y-v. y ;
1 W% D  w7 P7 {! X3 H4 F- U  a( s
return sqrt(dx * dx + dy * dy);
. d" I# s* D1 c  B; Q+ T, l
, W0 D! G8 e: g, A& j}

8 N) k- j$ `+ _* b/ N$ I+ 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 )，该函数实际求解最近点对。

程序14-10 预处理及调用c l o s e
4 {9 @1 u: v! ?" s( p
( Y( r4 X( ?" n8 p$ dbool closest(Point1 X[], int n, Point1&amp; a, Point1&amp; b, float&amp; d)
6 h( v1 {* w3 A# c1 l3 S0 C+ l* J
- F3 e5 B: c" W, @; r: K* o2 c{// 在n &gt;= 2 个点中寻找最近点对
8 g" `( n9 `4 s7 K4 G: |* r! D- d
// 如果少于2个点，则返回f a l s e
7 b2 G; @) y% t  o
# @' e2 _' K/ z* A; ]/ i# c// 否则，在a 和b中返回距离最近的两个点
+ m* Z+ U, S0 P4 P' O; E
if (n &lt; 2) return false;
" `: P* X  K5 D/ g4 T
// 按x坐标排序

) |1 P! ~& p7 O% `5 s8 h, OM e r g e S o r t ( X , n ) ;

// 创建一个按y坐标排序的点数组
1 x! b6 k' x% p; i1 D! H, L7 k
0 j% I: l3 {) }/ D! @Point2 *Y = new Point2 [n];

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

// 将点i 从X 复制到Y

Y.p = i;
$ ^, g9 x" g, B6 R2 b6 A
Y.x = X.x;
0 ^3 p  c4 Q+ s4 ~; S# P/ J4 k7 f* W
, U7 {' p4 w) ^9 `+ g, J, h6 YY.y = X.y;
) q3 @* I0 ~( m! _2 h% R
% X; y6 d6 C4 H$ c}
- N- i" e, Y  \& P* o. a
9 z# W( O- D5 V# @! a0 XM e r g e S o r t ( Y,n); // 按y坐标排序

/ ~! k, l+ j& w/ D// 创建临时数组

9 ~5 ^( g+ ]$ ]! V% [( H% hPoint2 *Z = new Point2 [n];
' k0 {9 ?+ z: q' |- H! N
4 X( V! t, G) k" Y# @// 寻找最近点对
5 l' C( D5 g3 F  K
c l o s e ( X , Y, Z , 0 , n - 1 , a , b , d ) ;
- M+ }2 j' _9 m5 ~$ ]
) ]9 z7 [& _' I& M/ A// 删除数组并返回

delete [] Y;
; N, M) r. u" V2 \( P
7 L; V8 F& S3 Z0 h( ^/ i3 H/ Adelete [] Z;

5 Y" \- Y" o& ^2 M- l2 g& f$ t* vreturn true;
  M# u4 k2 w$ W
7 Q. D3 r* W# O$ \5 m}

程序1 4 - 11 计算最近点对
% c) ~5 U. o3 S8 U
3 x3 v! @( ?! A7 d4 Uvoid close(Point1 X[], Point2 Y[], Point2 Z[], int l, int r, Point1&amp; a, Point1&amp; b, float&amp; d)
6 g) t( Q, i& _$ m
: R3 V. H( N: X1 C& j: I6 k{//X[l:r] 按x坐标排序

7 z5 P' h+ x' x7 t//Y[l:r] 按y坐标排序
( n6 `4 y8 ^2 A- }  R5 l
if (r-l == 1) {// 两个点

6 U& L; v6 q; S8 k* @2 na = X[l];
" a" b  e5 X- l! r6 O
b = X[r];

+ F. f0 Z9 m% pd = dist(X[l], X[r]);

) ^/ f5 N' v1 Cr e t u r n ; }
% W, H7 w4 P4 V0 P$ `
2 R6 X- e( j' {, D; }! {if (r-l == 2) {// 三个点
/ f* V% M: s" J9 l  k+ _
; x* s0 Z6 y. U* N  Q* `// 计算所有点对之间的距离

float d1 = dist(X[l], X[l+1]);
% e/ `& k  c% ^' Z, {8 R. h
float d2 = dist(X[l+1], X[r]);
$ W! E! Y6 d: O. j
float d3 = dist(X[l], X[r]);

// 寻找最近点对
9 c( Y' T8 Q6 R
$ e" t# o4 E' Y0 `8 Mif (d1 &lt;= d2 &amp;&amp; d1 &lt;= d3) {

; E- I: T! W! i5 e+ K) H+ v1 ?/ }$ p! }a = X[l];

b = X[l+1];

d = d1;
& l1 M% \: B) ?* K
r e t u r n ; }

if (d2 &lt;= d3) {a = X[l+1];
$ ?" c8 d3 @5 U" ]+ V2 A
* E8 [  J7 H1 H2 |5 M" rb = X[r];

d = d2;}
' ~. b' B8 I! u9 f# t1 Y
else {a = X[l];

! A1 ?3 {- P3 a7 T; s8 Kb = X[r];
9 u! C/ D: x; ^. x
3 @! a# w9 V" z! B- Y$ |d = d3;}
  f# ?1 U4 N; ?' [( C- {  \1 A' p$ o
r e t u r n ; }

0 {  Y# F: {, W3 ~* C, ^2 |/ /多于三个点，划分为两部分
8 @( G0 p* d# p
( ~; p& J, z& o2 z/ Wint m = (l+r)/2; // X[l:m] 在A中，余下的在B中

5 F" ]3 M9 E3 {* |// 在Z[l:m] 和Z [ m + 1 : r ]中创建按y排序的表
+ N2 o  W+ J8 n& Z' l, @
int f = l, // Z[l:m]的游标

g = m+1; // Z[m+1:r]的游标
8 x/ x: p; e# H6 c; T
4 v/ X2 j* S; L3 k3 e* [for (int i = l; i &lt;= r; i++)

if (Y.p &gt; m) Z[g++] = Y;
$ Z, M2 G% V) n4 @/ U4 e7 A7 ^1 C
else Z[f++] = Y;

6 \  f& K3 ~2 p: n8 P# u0 f// 对以上两个部分进行求解

c l o s e ( X , Z , Y, l , m , a , b , d ) ;
8 [1 S% o7 g7 }, W4 E9 r
float dr;
/ u7 @7 d7 x" v8 Z7 }/ [
% s- _/ R: e6 j: W/ rPoint1 ar, br;

, I3 p" u2 n0 t( Q- Qc l o s e ( X , Z , Y, m + 1 , r, a r, b r, d r ) ;
- E& c$ I$ I5 b# c$ q; e" ~) z
// (a,b) 是两者中较近的点对
. E: n/ v7 w$ S5 ~5 F
; q% ^$ e1 l3 k9 Vif (dr &lt; d) {a = ar;

% L) d( ]8 r% ?% ?# X* Y5 E7 Ab = br;

$ z  `  }; q* f/ g1 a2 N( A1 Vd = dr;}
4 V6 M  M9 L6 E% z
M e r g e ( Z , Y,l,m,r);// 重构Y

/ /距离小于d的点放入Z
9 }1 S* [; ~5 @
$ ]8 r, H0 n, J1 Iint k = l; // Z的游标

$ n- ?! m$ u$ b" i. Lfor (i = l; i &lt;= r; i++)

if (fabs(Y[m].x - Y.x) &lt; d) Z[k++] = Y;
% _6 E8 N" ^  S( O
// 通过检查Z [ l : k - 1 ]中的所有点对，寻找较近的点对

- ^2 x" j3 a- |( \9 dfor (i = l; i &lt; k; i++){

; K, \: \( F0 W3 ~7 Hfor (int j = i+1; j &lt; k &amp;&amp; Z[j].y - Z.y &lt; d;
8 D: \6 [. k% b5 x- P! T
j + + ) {
/ `" ]2 j) P/ x. U
' L7 e4 h) @- K/ ]% mfloat dp = dist(Z, Z[j]);
5 P4 ~2 o) {+ y9 t
2 G0 L5 n. G. pif (dp &lt; d) {// 较近的点对

d = dp;

, l. W, u% W# a: T. |a = X[Z.p];
  \. m9 d' Y' y9 Z
b = X[Z[j].p];}

' Z2 C9 E8 I  R+ {1 T}

" M4 i0 g4 h8 {1 b$ S! Z}
+ q6 E# P0 s5 y' K
7 t, b, ?1 j# ]; D2 X}

函数c l o s e（见程序1 4 - 11）用来确定X[1:r] 中的最近点对。假定这些点按x 坐标排序。在Y [ 1 : r ]中对这些点按y 坐标排序。Z[ 1 : r ]用来存放中间结果。找到最近点对以后，将在a, b中返回最近点对，在d 中返回距离，数组Y被恢复为输入状态。函数并未修改数组X。
$ m. a( ?+ v" N3 B5 e5 h
首先考察“小问题”，即少于四个点的点集。因为分割过程不会产生少于两点的数组，因此只需要处理两点和三点的情形。对于这两种情形，可以尝试所有的可能性。当点数超过三个时，通过计算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+ S5 U/ w8 M) L& S; J
6 J  W6 A: m) e7 x2 V8 L; X# g  a为实现图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

8 E1 q! k+ X) R+ J$ O+ n还是在RB中）与Z[j] 配对来实现，其中i＜j 且Z [ j ] . y - Z [ i ] . y＜ 。对每一个Z [ i ]，在2 × 区域内所检查的点如图1 4 - 1 7所示。由于在每个2 × 子区域内的点至少相距。因此每一个子区域中的点数不会超过四个，所以与Z [ i ]配对的点Z [ j ]最多有七个。
: V: J" \. ?4 M! g) P! `! n
2. 复杂性分析
. V3 T3 z* s' a, o
& f5 b5 L. u3 r$ a4 q8 I- z! ~令t (n) 代表处理n 个点时，函数close 所需要的时间。当n＜4时，t (n) 等于某个常数d。当n≥4时，需花费(n) 时间来完成以下工作：将点集划分为两个部分，两次递归调用后重构Y，淘汰距分割线很远的点，寻找更好的第三类点对。两次递归调用需分别耗时t (「n /2ù」和t (?n /2?).
! K) S% u0 x& \( d0 R: n
这个递归式与归并排序的递归式完全一样，其结果为t (n) = (nl o gn)。另外，函数c l o s e s t还需耗时(nl o gn)来完成如下额外工作：对X进行排序，创建Y和Z，对Y进行排序。因此分而治之最近点对求解算法的时间复杂性为(nl o gn)。</P>