距离最近的点对

韩冰

<>给定n 个点（xi，yi）（1≤i≤n），要求找出其中距离最近的两个点。
2 E* d% F% ^- C6 j1 l+ E! f# }7 I
+ U0 P# q; O5 U. U1 f例14-7 假设在一片金属上钻n 个大小一样的洞，如果洞太近，金属可能会断。若知道任意两个洞的最小距离，可估计金属断裂的概率。这种最小距离问题实际上也就是距离最近的点对问题。
! K7 X* Q5 ~: f- T7 _0 ^( Z# T/ ~% u
通过检查所有的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进行递归求解。
  w, |: q0 g/ h9 w1 x, ], H2 t; N5 Q
9 b4 X- X0 X! _/ Q
if (n较小) {用直接法寻找最近点对
+ |9 E$ L! a4 K* K. t* T/ _; r% w
R e t u r n ; }

// n较大

/ a* ?3 g% H1 |1 P& |! p0 r- |将点集分成大致相等的两个部分A和B

确定A和B中的最近点对
) ^5 o' p. O  b% @* S
2 g$ p8 t& Y; ?4 W2 e5 g确定一点在A中、另一点在B中的最近点对

8 S% t/ |( n  _% R' a9 e3 P从上面得到的三对点中，找出距离最小的一对点

6 M$ w7 F  S: ^5 e% z图14-13 寻找最近的点对


为了确定第三种情况下的最近点对，需要采用一种不同的方法。这种方法取决于点集是如何被划分成A、B的。一个合理的划分方法是从xi（中间值）处划一条垂线，线左边的点属于A，线右边的点属于B。位于垂线上的点可在A和B之间分配，以便满足A、B的大小。
; ?# ~  I0 P% J* r+ Y$ w
例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。

1 M% w& b$ y* B$ T4 {  G. E# v设是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）。
4 s! Q% v& j/ ?
% v- z, O8 ]+ a  p6 B例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) 是最近的点对。
0 l/ i1 B! p4 _1 J
* ~" s  l, i8 d  [# k. n为了确定一个距离更小的第三类点，RA 中的每个点最多只需和RB 中的6个点比较，如图1 4 - 1 6所示。

1. 选择数据结构

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

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

程序14-8 点类

class Point1 {
0 z# `+ e7 U, H$ U6 K4 F9 Y; @
1 n2 e3 d# C! |3 _4 h5 tfriend float dist(const Point1&amp;, const Point1&amp;);
6 ?. K+ s6 |  I0 G% `
friend void close(Point1 *, Point2 *, Point2 *, int, int, Point1&amp;, Point1&amp;, float&amp;);
+ t3 A/ k+ k9 u( u% `1 Q
& O8 Z7 B: U4 ofriend bool closest(Point1 *, int, Point1&amp;, Point1&amp;,float&amp;);

, M3 T  o2 m+ x5 d+ X- ]friend void main();

& I8 {6 J5 j6 Xp u b l i c :

% w1 O9 E) E" @6 k, Yint operator&lt;=(Point1 a) const
/ l- W$ u2 o6 Y9 L0 b4 z
7 [% l4 _  [' O$ t0 i% `2 |+ ^' o/ j{return (x &lt;= a.x);}

p r i v a t e :
$ [; ?. ~5 v- V9 I9 ^$ ?1 t
int ID; // 点的编号

) C) M, Y2 {/ i& A# H- P1 wfloat x, y; // 点坐标
' |/ e) W2 `& C: B" D
8 Q  d5 W: X: R4 }2 m' W$ K/ b7 K} ;
, ~, `& Q/ z, b) A% r' g
) d2 }  T2 a, E+ t5 Rclass Point2 {

friend float dist(const Point2&amp;, const Point2&amp;);
) k0 p9 j5 S# ]# n% D) z
$ v+ V1 \8 D' e  i  E! Rfriend void close(Point1 *, Point2 *, Point2 *, int, int, Point1&amp;, Point1&amp;, float&amp;);

( y8 Y% v5 k5 N1 s* Vfriend bool closest(Point1 *, int, Point1&amp;, Point1&amp;, float&amp;);

2 O: B: g* P; o% s7 bfriend void main();

! U! f0 V. q1 a4 Lp u b l i c :

% d$ ?! g# W! g+ hint operator&lt;=(Point2 a) const

) h4 F' i* K9 H{return (y &lt;= a.y);}
  y1 W+ W9 @$ j  _$ @
p r i v a t e :

- _1 w1 I8 z' C) `6 T1 Gint p; // 数组X中相同点的索引
$ q! R7 S) I1 `0 [- f* f& M  E
float x, y; // 点坐标

} ;

% c1 f  W  @: w: j6 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中的对应点。
. |# Y' k' E! Y7 R" I7 ?& l
7 }; k1 [& Y. a确定了必要的数据结构之后，再来看看所要产生的代码。首先定义一个模板函数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类的友元。
. K' F! E+ }& u7 L* H* R4 r
0 i0 `6 t' t8 f5 v' X$ b程序14-9 计算两点距离
, f6 ~* x+ a$ |8 {* B5 V" ~* R& s
template<CLASS T>

inline float dist(const T&amp; u, const T&amp; v)
  E) r* M' s1 t1 D& n. B
) g2 G: l$ D+ z0 d3 Y* [1 S{ / /计算点u 和v之间的距离
" X% q# S, g' ~8 ~5 H
6 w8 O, p' n; j9 B$ N9 W3 \+ h( T, zfloat dx = u.x-v. x ;
* d# P! e+ j0 |( o! g$ v) D* A7 l# c
float dy = u.y-v. y ;
& ]: x% f. P1 ~5 {
  K5 }1 ]+ \: m) N2 ?$ h( ereturn sqrt(dx * dx + dy * dy);
: X: A; L' c7 {0 T. ]
}

如果点的数目少于两个，则函数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
3 J! }! I* B9 X* w& H+ ^
bool closest(Point1 X[], int n, Point1&amp; a, Point1&amp; b, float&amp; d)

) M* }: Z# _- O) x{// 在n &gt;= 2 个点中寻找最近点对

% M$ d  f4 T* V+ U( v// 如果少于2个点，则返回f a l s e
8 W- {! `1 |% [: H- F; ?0 R- o
// 否则，在a 和b中返回距离最近的两个点

if (n &lt; 2) return false;

// 按x坐标排序
$ o% a( z- l8 F" H: h0 `0 Q  p" B2 x: N
M e r g e S o r t ( X , n ) ;

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

Point2 *Y = new Point2 [n];

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

4 {$ l6 _/ M0 {: m+ v# c// 将点i 从X 复制到Y

5 |! `* B8 L" A2 Q& _/ K0 cY.p = i;
' _& J: Z. ^# W2 B! U+ K5 V& {
3 j. l+ P- Q/ T1 I9 ^" aY.x = X.x;

Y.y = X.y;
  C9 b% v' i+ h
4 {( L& M  C: z5 T% w& y" b' Z5 |}

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

: H! ~, ^. p* C3 Y// 创建临时数组
2 o/ i' U, m+ k
. E% D9 @1 s% B$ `$ lPoint2 *Z = new Point2 [n];

2 u9 Z8 ]$ K  @2 O; N9 j  Q// 寻找最近点对

c l o s e ( X , Y, Z , 0 , n - 1 , a , b , d ) ;
3 u& z  f) r9 z
// 删除数组并返回
: m) u$ M2 f: c- p$ h
& \6 y4 |0 o: w4 edelete [] Y;
3 g  C' ^& C7 M/ m/ M7 p$ v
5 r5 \: B+ \& X3 @delete [] Z;
: C+ o5 a! S, [/ T
! Y1 g  n# l1 Oreturn true;
7 f  a9 Y7 h" `' T$ V
% H! y. w' l8 Q# |  K}

: `9 t  G  B! Z0 b. O程序1 4 - 11 计算最近点对
* K, ]- G% ~, H7 ~3 V' \
void close(Point1 X[], Point2 Y[], Point2 Z[], int l, int r, Point1&amp; a, Point1&amp; b, float&amp; d)
, L; z5 b5 R# T9 L
{//X[l:r] 按x坐标排序

//Y[l:r] 按y坐标排序

if (r-l == 1) {// 两个点
( ^/ w4 F3 G7 Q* M
  K- m9 w1 J% o8 v* F! y1 Y, da = X[l];

b = X[r];
3 c- |$ X4 T- Y9 Q1 w
d = dist(X[l], X[r]);

1 |$ ]9 ~9 c- k) v8 y5 q1 c( br e t u r n ; }
- d, ^/ p+ g0 p. r7 ]' `  J
if (r-l == 2) {// 三个点
9 w, w; N9 l) t; U/ q* ]
// 计算所有点对之间的距离

$ x, L; @7 U7 g: S  _/ x, \float d1 = dist(X[l], X[l+1]);

4 v" `. y' `; N$ s2 A, sfloat d2 = dist(X[l+1], X[r]);

float d3 = dist(X[l], X[r]);
5 K2 u1 k- v3 _3 k) B4 \# d. ^
// 寻找最近点对
  x+ ~* |( g4 O: t3 G
+ _' [: f- h& B- r& K/ Pif (d1 &lt;= d2 &amp;&amp; d1 &lt;= d3) {

a = X[l];
, C: q; c$ `: s! p
  ]7 J. U8 v: l3 f' R0 H. Ob = X[l+1];

# D* Q4 Z* O* q; cd = d1;

r e t u r n ; }

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

b = X[r];
  b0 E7 w# H0 D8 n. P
d = d2;}
( S! M6 [% ]# N
else {a = X[l];

b = X[r];

d = d3;}

. C0 M) L, P0 D5 kr e t u r n ; }

/ /多于三个点，划分为两部分
" f* f1 O5 K+ t/ J. ]2 F
) F6 E7 F3 y! P  I& t' w! e( F* [int m = (l+r)/2; // X[l:m] 在A中，余下的在B中

% t+ M( @  v% a- T9 I9 T0 l// 在Z[l:m] 和Z [ m + 1 : r ]中创建按y排序的表

int f = l, // Z[l:m]的游标
2 ^( g) }, B* I+ M9 S1 l9 {8 }6 ?
g = m+1; // Z[m+1:r]的游标

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

if (Y.p &gt; m) Z[g++] = Y;

) V% O7 G% V; X( V9 @7 \( @else Z[f++] = Y;

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

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

& G" ~. C9 I4 \4 _& Sfloat dr;
9 L9 m% {' V1 V% p; e
Point1 ar, br;

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

0 |3 d; y: b$ R" D) n// (a,b) 是两者中较近的点对
& o' `- b" D5 ^8 c2 O* `
if (dr &lt; d) {a = ar;
5 @7 x8 N( H! {! L; [7 K, o
' X# w2 w" {* `5 v: G  Gb = br;
8 \# U, h$ a! o5 T; L/ C
d = dr;}

M e r g e ( Z , Y,l,m,r);// 重构Y
; C1 w9 F+ @' X0 m  u( x( O6 F
* n$ E9 h" z6 z# i; G' T5 ^0 o/ /距离小于d的点放入Z
: a% o& M. ^9 C# l
3 D  g! W6 x9 Y1 pint k = l; // Z的游标

2 M" g1 j: D* Hfor (i = l; i &lt;= r; i++)
/ p+ c7 E1 o1 I1 J
if (fabs(Y[m].x - Y.x) &lt; d) Z[k++] = Y;
  G$ Q0 C. L1 @, Z
// 通过检查Z [ l : k - 1 ]中的所有点对，寻找较近的点对

; c* D% r8 f7 ifor (i = l; i &lt; k; i++){

for (int j = i+1; j &lt; k &amp;&amp; Z[j].y - Z.y &lt; d;
" @% G& H1 z, {! q
/ E1 b1 @. v0 ]% J: N0 Nj + + ) {

float dp = dist(Z, Z[j]);
2 ^. N+ e) V" ?5 y- I3 `
if (dp &lt; d) {// 较近的点对

d = dp;
% k. ^4 X. s5 i+ D2 ~
a = X[Z.p];

b = X[Z[j].p];}

}

}

}

函数c l o s e（见程序1 4 - 11）用来确定X[1:r] 中的最近点对。假定这些点按x 坐标排序。在Y [ 1 : r ]中对这些点按y 坐标排序。Z[ 1 : r ]用来存放中间结果。找到最近点对以后，将在a, b中返回最近点对，在d 中返回距离，数组Y被恢复为输入状态。函数并未修改数组X。
. V# t6 f0 f8 O' J: O
! D* m, L) A! d& d5 H0 c. d9 p首先考察“小问题”，即少于四个点的点集。因为分割过程不会产生少于两点的数组，因此只需要处理两点和三点的情形。对于这两种情形，可以尝试所有的可能性。当点数超过三个时，通过计算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

, A; Y) q3 P4 B, E6 e, n. ]还是在RB中）与Z[j] 配对来实现，其中i＜j 且Z [ j ] . y - Z [ i ] . y＜ 。对每一个Z [ i ]，在2 × 区域内所检查的点如图1 4 - 1 7所示。由于在每个2 × 子区域内的点至少相距。因此每一个子区域中的点数不会超过四个，所以与Z [ i ]配对的点Z [ j ]最多有七个。
* b+ R. v, B9 [2 r
2. 复杂性分析
7 H2 f% ]9 m! b
4 W4 G& F5 U1 x0 y3 D8 M令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>