距离最近的点对

韩冰

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

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

, j3 h: i4 l! 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进行递归求解。
; _. f' R3 |% T# U9 k7 v$ \
4 {$ ?: L/ b; H/ o! x( B- F" r0 ?- n: s
if (n较小) {用直接法寻找最近点对

R e t u r n ; }
. p1 p# w& [* R; n- ~0 t; R
& N7 ^- K" }# e" s; s1 Z; m// n较大

0 ]  [+ w( L' U6 l将点集分成大致相等的两个部分A和B
0 f& C5 s! |8 Z8 F5 e
! o" W% ?0 p( ]9 B确定A和B中的最近点对

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

0 z$ H* \: W; L9 u% x从上面得到的三对点中，找出距离最小的一对点

  X  a: P* s+ n# o" ?2 h/ g图14-13 寻找最近的点对
! T6 {5 v# d9 p: f: t5 v
  b4 A/ p' v3 g) X6 k' u
9 x4 r9 C, x# B) E3 z  f9 d# F为了确定第三种情况下的最近点对，需要采用一种不同的方法。这种方法取决于点集是如何被划分成A、B的。一个合理的划分方法是从xi（中间值）处划一条垂线，线左边的点属于A，线右边的点属于B。位于垂线上的点可在A和B之间分配，以便满足A、B的大小。
( ?' y# R: O% f/ e! _! [$ N
例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。
3 X& u# }4 N, Z, c  U! i( J1 Q7 X9 O8 c
, ~2 @' f! X; k% n( w设是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）。
* q5 R! @' _8 p4 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) 是最近的点对。

为了确定一个距离更小的第三类点，RA 中的每个点最多只需和RB 中的6个点比较，如图1 4 - 1 6所示。
, O  }4 W6 }( V0 O4 m% c0 J4 k. ?
1. 选择数据结构
# C8 S+ p5 H% v! [3 U4 m8 `& s* H! o
为了实现图1 4 - 1 3的分而治之算法，需要确定什么是“小问题”以及如何表示点。由于集合中少于两点时不存在最近点对，因此必须保证分解过程不会产生少于两点的点集。如果将少于四点的点集做为“小问题”，就可以避免产生少于两点的点集。

每个点可有三个参数：标号， x 坐标，y 坐标。假设标号为整数，每个点可用P o i n t l类（见程序1 4 - 8）来表示。为了便于按x 坐标对各个点排序，可重载操作符＜=。归并排序程序如1 4 -3所示。
' W8 c% k# o) A/ a* i
" d5 ~" W) ]/ ~程序14-8 点类

* P. V" A& j/ m2 H/ [3 o3 P2 Rclass Point1 {

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

' D; [- c% d, Z( O4 ^# vfriend void close(Point1 *, Point2 *, Point2 *, int, int, Point1&amp;, Point1&amp;, float&amp;);

friend bool closest(Point1 *, int, Point1&amp;, Point1&amp;,float&amp;);
$ b, }/ o7 E0 e% Y" n, H
friend void main();
7 u4 w+ h7 Q5 R' C6 `$ {- u4 {
p u b l i c :
  Z0 p0 X) L- d$ l# s& L& s$ @6 t# F
int operator&lt;=(Point1 a) const
  r- R# f7 ~) B) x$ p4 d) H
! A+ X6 m* @4 E{return (x &lt;= a.x);}
  @; I2 W. e& K
  j% \  L8 f) a( ?$ r" }p r i v a t e :
$ ]& N  Q, j1 y$ g( G7 A# a/ U
int ID; // 点的编号
; A% c9 P6 O# u! m/ N9 u
float x, y; // 点坐标
% r% w$ m" J; N6 ]4 ^( K% d) ^
- R9 Q4 c: C: C% u0 q9 M7 h8 m1 Y} ;
" }. F% o8 A) h, Z0 U
. B9 g: V; p- @# ]0 Hclass Point2 {
/ `7 L/ `0 P, q
friend float dist(const Point2&amp;, const Point2&amp;);
% p3 |7 ^" D8 F* y' w/ U
" c% P. P4 s. L1 b2 G2 f+ Ifriend void close(Point1 *, Point2 *, Point2 *, int, int, Point1&amp;, Point1&amp;, float&amp;);
; X0 \, M. h. T/ X  H
; V& C2 J% f6 i, D4 r, z' T; d9 Mfriend bool closest(Point1 *, int, Point1&amp;, Point1&amp;, float&amp;);

friend void main();

p u b l i c :
/ f% ~$ m* |7 D9 x" x' a8 K* N8 D
int operator&lt;=(Point2 a) const
' ~, L, U/ u5 a; l4 c
{return (y &lt;= a.y);}
+ E+ P& u8 A7 W1 C- v
p r i v a t e :
7 l5 |- r  k$ z$ b3 J
int p; // 数组X中相同点的索引
4 C: t! k3 f, X- u6 u0 q* w
* p5 T2 L9 A9 k( O! B5 cfloat x, y; // 点坐标

5 D5 n  I$ s- L* F9 m} ;

所输入的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 计算两点距离

template<CLASS T>
( Y3 s7 S4 F! v- a9 H5 N5 _
  p% q/ e- y$ ^inline float dist(const T&amp; u, const T&amp; v)
: V6 e1 I5 x) I% f$ g- |1 Z1 |
{ / /计算点u 和v之间的距离
. Q% r: e% F2 h
& y; c) C4 R$ d; f* k1 r8 Kfloat dx = u.x-v. x ;
! L7 C" V/ V4 c& r2 `8 u/ W
float dy = u.y-v. y ;
6 U- m- x5 K( J8 i. f* p
3 V; E' i" q& o+ Ereturn sqrt(dx * dx + dy * dy);

}

$ P2 S5 L, O6 V5 c3 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 )，该函数实际求解最近点对。

程序14-10 预处理及调用c l o s e

bool closest(Point1 X[], int n, Point1&amp; a, Point1&amp; b, float&amp; d)
7 _- P7 t! ^, f2 P1 b+ P
{// 在n &gt;= 2 个点中寻找最近点对

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

& b$ D" p7 |9 b7 c& E4 \$ Q2 @// 否则，在a 和b中返回距离最近的两个点
; d% j8 k: g& z. p: g: t, V9 D
( S2 d& t5 ^1 C! S- ?* `' m7 Tif (n &lt; 2) return false;
/ i% o3 l$ Z! _' _
// 按x坐标排序

6 z% @3 C9 k5 s  }4 K. iM e r g e S o r t ( X , n ) ;
5 P9 Z3 S+ p. ~! g. ~$ g( M4 P, B
" [. J3 C7 j* {// 创建一个按y坐标排序的点数组

# Y" v$ T" A4 z6 [  N5 xPoint2 *Y = new Point2 [n];
$ S: u! F9 q: s- B/ C
# L' Y9 R# ]$ afor (int i = 0; i &lt; n; i++) {
7 \' ?$ ~' V2 x$ L( H4 R6 N. o( h
// 将点i 从X 复制到Y
- [* x# K4 o: j5 I
Y.p = i;
9 u& M5 m- M; t4 d) u9 Z' E
Y.x = X.x;

Y.y = X.y;
  R% b: A5 z6 K, N+ u9 r
}

M e r g e S o r t ( Y,n); // 按y坐标排序

// 创建临时数组
- d, K7 {7 f! ^- w3 P: w: Z
6 \! z$ z2 q. m& j# e0 }$ cPoint2 *Z = new Point2 [n];
" R8 Z2 A$ [0 _: R
  s, M6 K; p# `, \. v& y- m// 寻找最近点对

c l o s e ( X , Y, Z , 0 , n - 1 , a , b , d ) ;
% K- \5 J$ v5 x8 ]/ U: w2 x1 H- x) m
& b& b6 U: c# `# H, e// 删除数组并返回
! f7 I: |" Z. r1 y6 n; }& S& o
delete [] Y;

6 J9 ]- C, ]( G% j+ Z. |delete [] Z;
, l' g0 z1 S6 H/ e% P: \3 J
return true;

8 C& W( M0 K4 r1 G: p}
% r' P- f) z4 o5 h  }1 G
5 f" z& J5 t+ ^/ k# H7 B" w" R# c程序1 4 - 11 计算最近点对
$ }( u; L/ b* p7 n( N% E
* r9 h* K6 N" x" Zvoid close(Point1 X[], Point2 Y[], Point2 Z[], int l, int r, Point1&amp; a, Point1&amp; b, float&amp; d)

/ d1 c9 U( b2 C{//X[l:r] 按x坐标排序

! u' E; p* ?9 v; q7 f//Y[l:r] 按y坐标排序
" ^0 @0 T( f( H8 w
% ?: R" C3 v- M) j* t  ~% `if (r-l == 1) {// 两个点
2 Q" K  n0 x2 @
a = X[l];
: u; G$ a$ w, F+ g+ m
; C6 t4 M+ n, ~6 }  T( [, ~b = X[r];

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

r e t u r n ; }

( B8 _# q* j7 ]5 F7 xif (r-l == 2) {// 三个点

2 n( o/ H; g3 m// 计算所有点对之间的距离
$ l7 Z+ P: Q& t  w8 D3 [# b; _4 N
float d1 = dist(X[l], X[l+1]);
% {; l3 I, B; V$ m8 K
float d2 = dist(X[l+1], X[r]);
& r; s! R6 G% ]) D* C" n
, I, d% ]: v6 D% I1 }- X( _, Mfloat d3 = dist(X[l], X[r]);

) O; `4 q6 S) e" g* ^3 J6 \' m0 M// 寻找最近点对
, U% H6 _+ N' Q
- B# Y1 p! Y+ B- l* O& l/ g; nif (d1 &lt;= d2 &amp;&amp; d1 &lt;= d3) {

a = X[l];
4 ]/ u7 ?! e& G! L) K4 l
; N, @' N" s, g1 |  h% }b = X[l+1];

d = d1;
7 `3 O7 E% O5 y# O5 e
r e t u r n ; }
+ }% {* {3 L9 x
5 w: x$ H& q' O  l7 I0 gif (d2 &lt;= d3) {a = X[l+1];
( `, @) Y6 F) u8 \! k" F
b = X[r];
$ q+ [( N8 {# u8 h$ P
d = d2;}
5 k+ `  a! [8 n( v# f
else {a = X[l];

b = X[r];
0 ?) G) _& q- T$ x# d4 m* s
: E7 o" G! k( F5 w0 l( M$ k& od = d3;}

) |3 @+ m- r- s9 K. ^" Q! Tr e t u r n ; }

4 Q, ~) ]; A1 j, x% T+ h9 |: W5 [/ /多于三个点，划分为两部分

int m = (l+r)/2; // X[l:m] 在A中，余下的在B中
8 b1 P5 }: K, l/ Z2 I% H- M" u
// 在Z[l:m] 和Z [ m + 1 : r ]中创建按y排序的表

int f = l, // Z[l:m]的游标
6 S0 C# I  F, f& V% f2 q  x
% S* K3 K; m+ ?3 Z( K) c# G7 Q9 ?g = m+1; // Z[m+1:r]的游标
% [& X  ^: P: a; `
for (int i = l; i &lt;= r; i++)
8 `4 z# d6 `  ^6 l5 v- u
, M) g& Q. e$ F" d9 c  Q  M" dif (Y.p &gt; m) Z[g++] = Y;
! D2 A( \! T9 O& y' I
! n" b' {  P! h6 R4 T$ Yelse Z[f++] = Y;
, P6 j( J5 P& C" A% G: B) M# O
// 对以上两个部分进行求解

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

float dr;
% K) \5 X* Q$ H% c% _
  l. Z) a* l( S' d, vPoint1 ar, br;
" m: x) r! n" y8 r( l
c l o s e ( X , Z , Y, m + 1 , r, a r, b r, d r ) ;
% c& M6 L8 C) Y8 l# u3 f
6 k; `5 j0 a# @) [% h# m6 A// (a,b) 是两者中较近的点对
+ g7 M9 a2 r& z
- Q' q" X" W+ Jif (dr &lt; d) {a = ar;

b = br;
6 R2 L* g8 D. W& q. i7 Q+ \" _7 j
$ q1 P+ H6 c- ?3 U* Xd = dr;}
& d' b7 u. \; M7 U4 T8 J4 b! B4 W; P8 ]
M e r g e ( Z , Y,l,m,r);// 重构Y

9 E$ h9 S4 ^" E" \6 d* f/ /距离小于d的点放入Z
8 Z% W, S- s) j0 M
2 {6 {6 j4 ~6 Fint k = l; // Z的游标
5 b3 X! W2 X5 |  {9 |: W! r! a* F3 {
for (i = l; i &lt;= r; i++)
7 e7 q5 i$ w1 g$ k' S2 L
2 [9 n& u0 {6 nif (fabs(Y[m].x - Y.x) &lt; d) Z[k++] = Y;
2 l$ f' W2 C7 F, j( O
( v" ^- }' ~. i# V# ^// 通过检查Z [ l : k - 1 ]中的所有点对，寻找较近的点对
: n0 O' U5 \1 u8 f2 W$ `
for (i = l; i &lt; k; i++){
; z7 }" [) K9 m; F
for (int j = i+1; j &lt; k &amp;&amp; Z[j].y - Z.y &lt; d;

9 o' R" P; ]: }- g& y# m4 ~6 A% n9 oj + + ) {
+ a5 L4 W7 }9 {
float dp = dist(Z, Z[j]);
+ d+ r' K) B8 U" u/ X
if (dp &lt; d) {// 较近的点对
$ D9 v/ o) y6 \" {2 Q1 M
d = dp;
% D- N% c5 _" ^( R
a = X[Z.p];

: t+ m# b5 X+ F, h! Ob = X[Z[j].p];}
. K; _( T6 R. ^& u. `# I+ b& O) O  e
}
3 B! M+ [' |2 q6 D1 d. M
}
8 l, h3 X- y4 R1 E3 d# y
" {1 h' d8 O- Z* O4 B}
/ r! ~# R( J: T( S  b2 d5 y
- Y) i/ H6 B+ Q; u函数c l o s e（见程序1 4 - 11）用来确定X[1:r] 中的最近点对。假定这些点按x 坐标排序。在Y [ 1 : r ]中对这些点按y 坐标排序。Z[ 1 : r ]用来存放中间结果。找到最近点对以后，将在a, b中返回最近点对，在d 中返回距离，数组Y被恢复为输入状态。函数并未修改数组X。
' F( I& ?& P5 n: O
# Z2 ~: e2 R5 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 ]。

# d; ^& Q& `; R; c; C& 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 ]最多有七个。
: b/ A& Y( j; O$ m2 S$ C
2. 复杂性分析

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

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