距离最近的点对

韩冰

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

例14-7 假设在一片金属上钻n 个大小一样的洞，如果洞太近，金属可能会断。若知道任意两个洞的最小距离，可估计金属断裂的概率。这种最小距离问题实际上也就是距离最近的点对问题。
, {2 [/ Q0 J% h1 C: i2 g
通过检查所有的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进行递归求解。


; V$ E) `, o8 ?4 X% yif (n较小) {用直接法寻找最近点对

R e t u r n ; }

3 m' s% [% R3 F. X: e// n较大

将点集分成大致相等的两个部分A和B
( v; @% \  l  n6 D2 K
确定A和B中的最近点对
0 @" i5 T! K; w# R8 w/ P: c3 C7 G
确定一点在A中、另一点在B中的最近点对
% z: z, A* a: W/ B
( e  W! j. Q5 y9 [' F9 l6 ]从上面得到的三对点中，找出距离最小的一对点

; G7 F+ B- |* F- _" @6 R5 N. }, w图14-13 寻找最近的点对


1 l+ M; p+ S( q. M. ^. j7 H: G# f; s为了确定第三种情况下的最近点对，需要采用一种不同的方法。这种方法取决于点集是如何被划分成A、B的。一个合理的划分方法是从xi（中间值）处划一条垂线，线左边的点属于A，线右边的点属于B。位于垂线上的点可在A和B之间分配，以便满足A、B的大小。
9 C6 y' T8 `3 J/ [
0 C1 C& \) d# s8 D! L# {9 j例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。
" n3 ?: x& b, M5 `* ]
设是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) 是最近的点对。

7 x2 V5 U3 M$ B1 l; z为了确定一个距离更小的第三类点，RA 中的每个点最多只需和RB 中的6个点比较，如图1 4 - 1 6所示。
) V8 v9 `% ~, v4 O1 {( j7 s
8 e, ~) w( I8 F2 H1. 选择数据结构

" u4 i) J7 c" W# s为了实现图1 4 - 1 3的分而治之算法，需要确定什么是“小问题”以及如何表示点。由于集合中少于两点时不存在最近点对，因此必须保证分解过程不会产生少于两点的点集。如果将少于四点的点集做为“小问题”，就可以避免产生少于两点的点集。

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

0 ?  E; z' P- \# Z' U1 r0 y程序14-8 点类
$ B$ u- v, y  N& \
class Point1 {
4 I" m1 Y$ @7 w: s
  @" }& x4 O3 e5 d4 T$ V+ J/ \friend float dist(const Point1&amp;, const Point1&amp;);

friend void close(Point1 *, Point2 *, Point2 *, int, int, Point1&amp;, Point1&amp;, float&amp;);

) D- p" N7 m1 k# zfriend bool closest(Point1 *, int, Point1&amp;, Point1&amp;,float&amp;);
7 i& n1 r1 w1 l
friend void main();

p u b l i c :
+ m* u* C: n& X8 i5 U! x/ E" }
int operator&lt;=(Point1 a) const

2 ~7 D: K& F5 L: x- p{return (x &lt;= a.x);}

p r i v a t e :

int ID; // 点的编号

float x, y; // 点坐标
, y7 z& h7 l4 q7 O0 P. X" O
# w0 v' Z0 J# W} ;

class Point2 {

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

. \. k' c, D( \friend void close(Point1 *, Point2 *, Point2 *, int, int, Point1&amp;, Point1&amp;, float&amp;);
5 G- W) M" W; j% i+ d# x% M$ w
, M4 Z2 {- E5 ?( r, O5 jfriend bool closest(Point1 *, int, Point1&amp;, Point1&amp;, float&amp;);
6 I* |7 E% \% A8 _- F: v8 a
friend void main();

8 k0 C( b- g2 ~p u b l i c :
7 u$ _) E! R7 k: S9 Y6 N- A
% Q8 u/ h0 T% h4 Q5 ?* ~3 xint operator&lt;=(Point2 a) const
0 O5 S, l4 ^- o5 m( J. I4 R
+ j: U- Y1 f! K4 T5 B9 ^{return (y &lt;= a.y);}

p r i v a t e :
5 r, e/ j- {( O$ p
* \: s3 m3 c: e* O# Jint p; // 数组X中相同点的索引

' E. ^/ v. n  Vfloat x, y; // 点坐标
) L7 \7 P% q1 J- K* y6 S* ^
} ;
8 i# m6 M5 |9 J8 \  k
  x# ]8 U; F! B4 g/ F: k; r所输入的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 l% O" O9 n  p- P- W
程序14-9 计算两点距离
! L2 _4 C1 p( y5 i9 F4 k
template<CLASS T>
) q, B  U8 a: g( a' q7 X2 y
inline float dist(const T&amp; u, const T&amp; v)

{ / /计算点u 和v之间的距离
2 L' h& T7 A& z7 a
float dx = u.x-v. x ;

float dy = u.y-v. y ;

return sqrt(dx * dx + dy * dy);
; [; ?3 i$ y2 M3 y( o0 h8 h! e
! e2 E9 M" ?1 O# l- \) @+ S}
! V; A; [# r) o9 `) i
7 p3 H% w/ A- A2 o如果点的数目少于两个，则函数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 )，该函数实际求解最近点对。
9 {+ P& N+ s, m) Z! ^* [. w
, I  [! P4 c# [8 ?2 s# a" u程序14-10 预处理及调用c l o s e
5 L$ k% M* ^" i
bool closest(Point1 X[], int n, Point1&amp; a, Point1&amp; b, float&amp; d)

$ u. L: y) c+ `6 H9 |& c% }0 V( k* t{// 在n &gt;= 2 个点中寻找最近点对

// 如果少于2个点，则返回f a l s e

7 y! i1 ]9 k& ~- S8 h& ]3 n// 否则，在a 和b中返回距离最近的两个点
1 @! d7 U$ n/ q1 _
if (n &lt; 2) return false;
8 x' {7 y1 J4 s' J! F- r6 [  N$ D
// 按x坐标排序

M e r g e S o r t ( X , n ) ;

# {! t' w$ I" s. Q// 创建一个按y坐标排序的点数组
  N6 ^% c. B4 ^  L( S% t, ~
+ r% U8 U4 y" }1 g4 f; k$ `6 _/ oPoint2 *Y = new Point2 [n];
& d/ k$ Y+ B( [% o/ T. h. q8 \! J
for (int i = 0; i &lt; n; i++) {
4 I; O8 e  r4 w0 C; \" ]
3 h4 \1 D, g6 F8 u// 将点i 从X 复制到Y
' G9 B, i9 _: l' H
% s( E  `& r( z# K2 N8 I5 rY.p = i;

Y.x = X.x;
% ^* K0 [- N% y" w& O* x
0 ~: F, k0 _8 Q" [; S% J, a7 m5 n! z6 HY.y = X.y;

9 ?8 p6 j8 U. A/ n  J7 F+ g3 t3 }}

0 ]+ g* F9 y7 n- A2 V& iM e r g e S o r t ( Y,n); // 按y坐标排序
" N* r' W' Q3 s: [3 H2 x
// 创建临时数组
: U6 G( a3 J2 z) K
' b, k6 L5 ^8 f* I1 OPoint2 *Z = new Point2 [n];

5 M" ^7 F% d. o- X8 T4 y: x7 ?" o- Z// 寻找最近点对

c l o s e ( X , Y, Z , 0 , n - 1 , a , b , d ) ;

// 删除数组并返回
# ?0 I' B9 M; D4 b/ I
  m$ S' l; B: v: Mdelete [] Y;
( L- m" c$ J. Y4 H' x
delete [] Z;

return true;

- E6 F/ ^. u1 B- s) h1 R6 r}
& A; M+ R& `4 C1 U; K' H. }
. \, l6 X7 Y' r程序1 4 - 11 计算最近点对

void close(Point1 X[], Point2 Y[], Point2 Z[], int l, int r, Point1&amp; a, Point1&amp; b, float&amp; d)
5 m  w0 V* p7 l, G. {, n) J5 U$ l
{//X[l:r] 按x坐标排序
/ g- d( w3 H" }
//Y[l:r] 按y坐标排序

if (r-l == 1) {// 两个点

( O" ~  @, \/ v; ]  W8 ga = X[l];

, \' Z6 T  o5 m# P4 P0 ?% ?b = X[r];
5 n' B3 [6 E) L+ d0 {
d = dist(X[l], X[r]);

6 H( s7 G% [/ i2 Rr e t u r n ; }
+ L, w: N0 F7 l  V. N9 e9 W! x
+ H5 }1 g- c3 q& w1 g0 M4 Tif (r-l == 2) {// 三个点
% y- d+ z; H7 j
8 H6 _* O' _4 H. J* V// 计算所有点对之间的距离

; {5 H1 Q6 I3 p0 E! Mfloat d1 = dist(X[l], X[l+1]);

3 b$ ]/ h' D8 S: ]$ B; L6 F- R! wfloat d2 = dist(X[l+1], X[r]);
) _- J9 B" \; e( y, o: o) c# o
float d3 = dist(X[l], X[r]);

// 寻找最近点对

4 w, B6 h" l6 O% fif (d1 &lt;= d2 &amp;&amp; d1 &lt;= d3) {
4 z- R  |5 ?5 H- R' N
a = X[l];
; S- \# I; X0 U2 l
b = X[l+1];
/ I/ H, ^$ L# g
d = d1;

2 ^5 P6 y! l! O/ t; l7 Y% w0 Er e t u r n ; }

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

/ A5 d" r  n4 xb = X[r];

% a. w# t& h( z& i  Dd = d2;}

else {a = X[l];
1 ]( w& S+ d% c3 b
b = X[r];

# S+ B4 r6 b5 od = d3;}

$ N! W5 s  [; v7 I3 qr e t u r n ; }

/ /多于三个点，划分为两部分
. m! @4 Z4 {1 d
int m = (l+r)/2; // X[l:m] 在A中，余下的在B中
9 W$ x- o* M( H: H
// 在Z[l:m] 和Z [ m + 1 : r ]中创建按y排序的表

" k  O9 z2 B0 ^. \( X1 n- b! e6 W. q* h) h+ Wint f = l, // Z[l:m]的游标
  s  L+ ~! X0 {' d2 L: x( \
g = m+1; // Z[m+1:r]的游标
5 n0 m/ P2 @+ G1 W1 j5 p4 ^! r
. s  @6 h2 v7 v5 p" J7 R5 s# Pfor (int i = l; i &lt;= r; i++)
! k: V( J+ @2 S# I1 X; K
- c- K4 y' n* u2 n; sif (Y.p &gt; m) Z[g++] = Y;
' G  Y+ Q- Y/ [
* u- h3 ], p! e) L& ?0 v1 nelse Z[f++] = Y;
+ S# [0 ^" [+ e9 X3 [( Y+ R! U
7 b2 v% e0 R+ X9 v6 v% o, N# z" E// 对以上两个部分进行求解
; D% u: S5 d( z2 ~3 Y
c l o s e ( X , Z , Y, l , m , a , b , d ) ;
+ I0 G3 P; K& M/ m9 h1 x! Y! m
float dr;

Point1 ar, br;

2 u2 ^8 X; S" q( Y" |0 \c l o s e ( X , Z , Y, m + 1 , r, a r, b r, d r ) ;
) `6 f9 _5 g9 X9 ^& P2 e
4 q" {3 l2 T3 @5 p) _% J9 g// (a,b) 是两者中较近的点对
; H2 [8 b) h/ }7 |& @
" g" K5 P* U* I3 xif (dr &lt; d) {a = ar;

. i& v9 ]9 `( Y' C$ q$ bb = br;

, h# G/ ]0 g0 I7 z5 Td = dr;}

; ?9 b5 N& w8 xM e r g e ( Z , Y,l,m,r);// 重构Y
" ~* O% }  }9 i4 f' i
& {- s' i+ ]+ Q& q8 ~/ /距离小于d的点放入Z

int k = l; // Z的游标
6 ?, Q1 s- w6 c/ C6 K6 `$ H
# w# c2 U6 K2 l1 q  n6 wfor (i = l; i &lt;= r; i++)
0 F& q, E( J5 |' e/ b; B5 n& F) H
# p: u( a' l' ?2 qif (fabs(Y[m].x - Y.x) &lt; d) Z[k++] = Y;
1 R  i3 q' D6 Z8 m3 G, F
& L3 f9 u3 E/ g5 Z// 通过检查Z [ l : k - 1 ]中的所有点对，寻找较近的点对
5 s" P4 @5 x! i" T/ ~+ v5 B0 G2 x4 j+ E
- J- O+ a0 A* W4 Z% m& V- ]% |for (i = l; i &lt; k; i++){
" i5 o' F4 g* K0 k1 x# ]
for (int j = i+1; j &lt; k &amp;&amp; Z[j].y - Z.y &lt; d;
+ V9 \: `9 T; H& r2 [( i
j + + ) {
5 i% V! P; ^9 l! V# p0 v
float dp = dist(Z, Z[j]);

8 T! X" A1 f" d+ Z+ j8 F# {- Bif (dp &lt; d) {// 较近的点对
' P2 x% a: U3 u+ `
* Z; g. X- ]4 Kd = dp;
" _  y4 z  v1 O/ z* I
a = X[Z.p];
$ f  \4 N4 [* ]/ c; ]! i
b = X[Z[j].p];}
7 u; l, k0 L; {; l3 S
+ x0 @* k8 ]: ]9 x3 T0 Z}

3 }2 t* }3 ?( @}
, `! X  I; f1 y( }
}
0 W; L7 ]$ u% ~$ R* M9 W
, O. J3 l7 {* }- E4 |函数c l o s e（见程序1 4 - 11）用来确定X[1:r] 中的最近点对。假定这些点按x 坐标排序。在Y [ 1 : r ]中对这些点按y 坐标排序。Z[ 1 : r ]用来存放中间结果。找到最近点对以后，将在a, b中返回最近点对，在d 中返回距离，数组Y被恢复为输入状态。函数并未修改数组X。

- K* I6 a& r3 E' W! b" x+ F. M首先考察“小问题”，即少于四个点的点集。因为分割过程不会产生少于两点的数组，因此只需要处理两点和三点的情形。对于这两种情形，可以尝试所有的可能性。当点数超过三个时，通过计算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 ]。
) l) L6 g  N. v7 Q1 W' ~2 Y
: O  m) h# _: K) G5 P9 z为实现图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
2 p2 n) e4 z" R# b- J& h$ d& E
还是在RB中）与Z[j] 配对来实现，其中i＜j 且Z [ j ] . y - Z [ i ] . y＜ 。对每一个Z [ i ]，在2 × 区域内所检查的点如图1 4 - 1 7所示。由于在每个2 × 子区域内的点至少相距。因此每一个子区域中的点数不会超过四个，所以与Z [ i ]配对的点Z [ j ]最多有七个。
, T6 [. e" [. ?7 O) x3 ~) E
2. 复杂性分析
7 |! M. `# w: n1 w$ l9 t
令t (n) 代表处理n 个点时，函数close 所需要的时间。当n＜4时，t (n) 等于某个常数d。当n≥4时，需花费(n) 时间来完成以下工作：将点集划分为两个部分，两次递归调用后重构Y，淘汰距分割线很远的点，寻找更好的第三类点对。两次递归调用需分别耗时t (「n /2ù」和t (?n /2?).

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