距离最近的点对

韩冰

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

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

7 }" R, V- T. V& F* N) P通过检查所有的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进行递归求解。


if (n较小) {用直接法寻找最近点对
) q" {; x9 z- B7 b( R
( v0 Z' N' s9 }0 }3 dR e t u r n ; }

// n较大
( t& p$ a. E0 o8 T* ]) x* R
- d5 b, L* s# T; J将点集分成大致相等的两个部分A和B

: q. x, [: C  q3 z* P" ^确定A和B中的最近点对

/ p2 o! F' [2 z' X% j. O! r& Q确定一点在A中、另一点在B中的最近点对

从上面得到的三对点中，找出距离最小的一对点

& n" C. ]. g& ^  _图14-13 寻找最近的点对
4 w" Q  d# Q7 e. D- C7 s

为了确定第三种情况下的最近点对，需要采用一种不同的方法。这种方法取决于点集是如何被划分成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。
% I% g- {$ u# D* Z+ r; K6 ~! d3 |- F8 i
设是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 ^- T, ~- R. T$ B+ ?& P例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所示。

) `1 u* O( i6 j$ X3 r1. 选择数据结构

- l7 L- g2 B) O# ?7 a( I, I9 `为了实现图1 4 - 1 3的分而治之算法，需要确定什么是“小问题”以及如何表示点。由于集合中少于两点时不存在最近点对，因此必须保证分解过程不会产生少于两点的点集。如果将少于四点的点集做为“小问题”，就可以避免产生少于两点的点集。

每个点可有三个参数：标号， x 坐标，y 坐标。假设标号为整数，每个点可用P o i n t l类（见程序1 4 - 8）来表示。为了便于按x 坐标对各个点排序，可重载操作符＜=。归并排序程序如1 4 -3所示。
- s0 |: o. R6 X  t. w$ a/ P) ?
程序14-8 点类
  @: i6 j5 ]) O+ [$ L7 |
- q, b2 Q" f. L* Z0 P1 Qclass Point1 {

friend float dist(const Point1&amp;, const Point1&amp;);

  K$ v! @  l) w  [, V9 [friend void close(Point1 *, Point2 *, Point2 *, int, int, Point1&amp;, Point1&amp;, float&amp;);
4 a- ~) D3 k  W$ E6 E' T
friend bool closest(Point1 *, int, Point1&amp;, Point1&amp;,float&amp;);
9 |( o9 T3 {6 n* Z
friend void main();

p u b l i c :
9 i9 i8 U  g0 n' s, m$ Q
. ?0 T9 o5 p2 k+ E1 Q' e% rint operator&lt;=(Point1 a) const
0 V( j4 K1 c9 V* U6 R
. C6 c: A. B  _8 i2 Y# L. k{return (x &lt;= a.x);}

p r i v a t e :

" \! l0 B' F) S3 s! b& k' hint ID; // 点的编号
& [0 a; ^8 U2 e
( R- Q/ Z5 I% \" yfloat x, y; // 点坐标

  h4 n. u5 }: K} ;
% m) f' d  j% Z$ ^' b6 L4 ]
class Point2 {
1 C  R  w; ^6 `$ i( X
3 V- o8 t: H4 F5 ]8 kfriend float dist(const Point2&amp;, const Point2&amp;);

friend void close(Point1 *, Point2 *, Point2 *, int, int, Point1&amp;, Point1&amp;, float&amp;);
/ `, X' f! C2 T) l  i, c
friend bool closest(Point1 *, int, Point1&amp;, Point1&amp;, float&amp;);
. Z2 O+ T! {) o1 M6 y) R* E
- d2 G: V+ y; H9 c3 Y" Q+ g" `, rfriend void main();

9 g7 k# Z% C4 k* F3 l4 dp u b l i c :
" N* b# i6 I% ~' `( ?7 {. E0 a4 a8 i
int operator&lt;=(Point2 a) const
. }0 O! m# d' J0 `
{return (y &lt;= a.y);}
* ?* L6 Q4 r, L! |4 I1 R
3 s* k9 V# \5 Wp r i v a t e :
: X2 R( n7 Q0 m
int p; // 数组X中相同点的索引

float x, y; // 点坐标
3 X9 m% E* [: h( m- [) ~/ d
} ;
9 A$ ~) C# N  k* B4 b* ^
+ [2 t/ R& v# q# Y) _- I9 L2 s所输入的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类的友元。
. q" J% g7 B! o4 m8 n& \/ S
( R& e4 z; {( O; O+ [! [8 R程序14-9 计算两点距离

7 t' M* q2 u) p+ B& Ctemplate<CLASS T>
; y1 r$ P4 e8 Q0 z! q
: h, \) A  ?% ~inline float dist(const T&amp; u, const T&amp; v)

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

7 o3 C* F& M/ P  E/ o$ w; tfloat dx = u.x-v. x ;

6 ?% A3 \) q2 |+ Yfloat dy = u.y-v. y ;

return sqrt(dx * dx + dy * dy);
6 m7 G1 E  E+ T7 V; y
}
2 [9 R5 _' d" g
如果点的数目少于两个，则函数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 )，该函数实际求解最近点对。
" N/ s$ z7 d  f: U/ C% U
程序14-10 预处理及调用c l o s e

4 n8 L3 b/ Q. Rbool closest(Point1 X[], int n, Point1&amp; a, Point1&amp; b, float&amp; d)

{// 在n &gt;= 2 个点中寻找最近点对
& g) |# Y( w+ r; R
* [1 d. O: o5 ~& K1 H0 o" g! }// 如果少于2个点，则返回f a l s e

/ f; }) K3 P- U: s/ H// 否则，在a 和b中返回距离最近的两个点
0 i6 ^8 n! u4 x: o5 S! |% u
2 N+ l; |4 s7 |% a2 {/ nif (n &lt; 2) return false;

// 按x坐标排序

* l0 U+ Q7 F2 W4 @; X. qM e r g e S o r t ( X , n ) ;
* P8 S: e0 N3 ?( x
// 创建一个按y坐标排序的点数组

Point2 *Y = new Point2 [n];
( \: O3 R: n$ L% E
% F7 u" c8 k6 d0 N/ B; xfor (int i = 0; i &lt; n; i++) {

1 ?$ R0 s$ u( B9 A" H+ [- Y1 |  x// 将点i 从X 复制到Y
% X1 U+ K) q# `6 n3 {4 Q
  h! L5 s4 v& @: XY.p = i;
# A+ I2 ~& q/ O3 Y9 K0 a; A
4 s9 {  m7 R' r5 }3 K% p  KY.x = X.x;

Y.y = X.y;
" \4 y# @5 |/ p/ r' r4 v" _
' [9 N0 f% V3 M( A4 w3 T}
1 G7 Y: B% j( ]. f+ H! H1 b
M e r g e S o r t ( Y,n); // 按y坐标排序
. c9 k( a0 z- H$ }: @3 A
// 创建临时数组

Point2 *Z = new Point2 [n];

& @' d# k! k& p/ h( r// 寻找最近点对

c l o s e ( X , Y, Z , 0 , n - 1 , a , b , d ) ;
  T/ X: J% p& W( a" Y+ M7 G& A& W
// 删除数组并返回

: D4 @9 z3 x+ i) d% p5 R/ K1 `6 odelete [] Y;

delete [] Z;
! [4 @) e& R  @* o. A
: h( m( z( k2 C5 Treturn true;

* P1 c# X$ ^1 S* x& c1 {}
8 v/ z  {' P) z8 L9 ~3 Y
程序1 4 - 11 计算最近点对
* v  X6 x! M: H3 p/ G' q& y6 a
void close(Point1 X[], Point2 Y[], Point2 Z[], int l, int r, Point1&amp; a, Point1&amp; b, float&amp; d)

# {2 @) `3 p, t# ~{//X[l:r] 按x坐标排序
. V8 ]0 j2 \. S  t! y
//Y[l:r] 按y坐标排序

. {6 O, X; y+ e- ]5 `3 W- Uif (r-l == 1) {// 两个点

& c/ I# k- W; f' P# Ua = X[l];

. \- b: e+ J/ L/ [4 Y5 q1 ^9 Eb = X[r];

d = dist(X[l], X[r]);

r e t u r n ; }
2 r5 b' h: d4 x) g- ~: G- X6 l
: _* _, z9 C: xif (r-l == 2) {// 三个点
8 T* D: p" O  e7 H* o( _
// 计算所有点对之间的距离
% V" y: j% A8 r5 s( g4 T
float d1 = dist(X[l], X[l+1]);
8 F: C' b6 e& ]* g  d
float d2 = dist(X[l+1], X[r]);
& _8 q* g  ]( m7 T7 B
float d3 = dist(X[l], X[r]);
& ]* J& j) [7 n! P" R
# F/ `% v6 a5 C# X  z# z// 寻找最近点对

+ x* V$ u9 i8 b" j) S/ Vif (d1 &lt;= d2 &amp;&amp; d1 &lt;= d3) {
: |( r* o$ X) h9 R
2 J5 Q* d. a, Ra = X[l];

" a% E/ W" X% Y$ A$ U& D6 L% B0 db = X[l+1];
% P* E4 W% A( F7 [- \  g
d = d1;
" c3 ~4 H4 _1 G
4 i4 Q8 y% J1 U5 f, Fr e t u r n ; }

' o0 w* d7 h! mif (d2 &lt;= d3) {a = X[l+1];
: q9 b" ^( m6 w" T0 D8 w9 z. j! s
( z; D9 b# k9 N6 z4 ]& Eb = X[r];
& _+ e/ o% b7 g2 D% P* ^( D' p
9 i- Y2 D7 k* Y$ b4 pd = d2;}

else {a = X[l];

: D5 B' ~: k% `6 t2 V% \) Gb = X[r];
: b7 n- N' F8 d3 k7 N
/ z+ }+ f: v: _2 Id = d3;}
% N' s& M, i; U4 l
2 l# Y8 h" M, ~# m# er e t u r n ; }

; z, q$ Q  a  s, D/ /多于三个点，划分为两部分

$ A, b" B& F& R) I* ~' gint m = (l+r)/2; // X[l:m] 在A中，余下的在B中
( v; R" T$ O, i$ G! @! u. ^0 h1 Q
// 在Z[l:m] 和Z [ m + 1 : r ]中创建按y排序的表

int f = l, // Z[l:m]的游标
% _. r: g* |' [6 [# H6 T4 O) }
( z8 g6 [, O; V2 [+ f- ^9 i* z% Mg = m+1; // Z[m+1:r]的游标

for (int i = l; i &lt;= r; i++)

1 w7 @3 t$ k+ Z/ D; \5 g. Rif (Y.p &gt; m) Z[g++] = Y;
7 Y" Y. `2 ]! z+ D9 {- x
  \( n6 f1 a: {7 k& J3 t, helse Z[f++] = Y;

- V, |9 n5 \6 U4 _3 Z- _0 U' i// 对以上两个部分进行求解
1 o* a$ Z9 u) }* w0 A6 Z
c l o s e ( X , Z , Y, l , m , a , b , d ) ;

float dr;

Point1 ar, br;

* t+ K1 Y. r1 [: sc l o s e ( X , Z , Y, m + 1 , r, a r, b r, d r ) ;

4 g2 X8 Y* B9 T0 x// (a,b) 是两者中较近的点对

if (dr &lt; d) {a = ar;
* b4 P# R  |7 p9 M' o! m
b = br;

d = dr;}

4 @, a, {& i# L/ {8 A0 @* uM e r g e ( Z , Y,l,m,r);// 重构Y
, m- R! n( j- w! z+ T8 h; G
6 K  C+ m' _0 j4 s/ /距离小于d的点放入Z

int k = l; // Z的游标

for (i = l; i &lt;= r; i++)
. o6 N  Z- I! E" W. C
if (fabs(Y[m].x - Y.x) &lt; d) Z[k++] = Y;
% h( x1 l' D. q2 }, B6 @
// 通过检查Z [ l : k - 1 ]中的所有点对，寻找较近的点对

for (i = l; i &lt; k; i++){

for (int j = i+1; j &lt; k &amp;&amp; Z[j].y - Z.y &lt; d;
% u9 I% N% {5 x/ H1 ?
& h1 N  I8 g, E2 J% X( L: d: \! ej + + ) {
4 G) Y. n& x& h" u: z2 A0 @2 j
float dp = dist(Z, Z[j]);

if (dp &lt; d) {// 较近的点对
! U0 \2 q+ [& e% F
+ C. V$ q/ Y: F8 B2 W7 X" V7 pd = dp;

a = X[Z.p];
/ t& a0 h$ g+ z- d" [+ S
+ Z, d) |& e- ?* d) Tb = X[Z[j].p];}

& n) t; r7 I9 R8 t) y9 x}

* y. Y2 H+ s6 R  d( k) m: v}
/ F4 Y. ~0 i* |- D3 X5 R" G
8 b* e( f* `9 A" D5 n, w: g6 \+ u3 y/ x}
! d% ?1 S' D8 }1 S8 @( _: t, A5 d) ?9 |
8 ]( G5 U% f. h函数c l o s e（见程序1 4 - 11）用来确定X[1:r] 中的最近点对。假定这些点按x 坐标排序。在Y [ 1 : r ]中对这些点按y 坐标排序。Z[ 1 : r ]用来存放中间结果。找到最近点对以后，将在a, b中返回最近点对，在d 中返回距离，数组Y被恢复为输入状态。函数并未修改数组X。
) Y! N1 o( z2 D9 @8 V
6 {' y$ c# i. O首先考察“小问题”，即少于四个点的点集。因为分割过程不会产生少于两点的数组，因此只需要处理两点和三点的情形。对于这两种情形，可以尝试所有的可能性。当点数超过三个时，通过计算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
0 n% m4 ^' x0 t: F) O4 l( u; h
8 s" [( T* I# Y# M/ r还是在RB中）与Z[j] 配对来实现，其中i＜j 且Z [ j ] . y - Z [ i ] . y＜ 。对每一个Z [ i ]，在2 × 区域内所检查的点如图1 4 - 1 7所示。由于在每个2 × 子区域内的点至少相距。因此每一个子区域中的点数不会超过四个，所以与Z [ i ]配对的点Z [ j ]最多有七个。

" N" h- k# _: c( u( F& Q' m# v2. 复杂性分析

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

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