2.11 2D平面最近点对问题[closest pair problem]

【本文链接】

http://www.cnblogs.com/hellogiser/p/closest-pair-problem.html

【题目】

给定平面上N个点的坐标,找出距离最近的两个点之间的距离。

【蛮力法】

对于n个点,一共可以组成n(n-1)/2对点对,对这n(n-1)/2对点对逐对进行距离计算,通过循环求得点集中的最近点对。时间复杂度为T(n)=n^2。

C++ Code 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
 
/*
    version: 1.0
    author: hellogiser
    blog: http://www.cnblogs.com/hellogiser
    date: 2014/7/11
*/

struct Point
{
    
double x;
    
double y;
}

double distance(const Point &a, const Point &b) const
{
    
// distance of point a and b
    return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}

double ClosestPairBruteforce(Point P[], int n)
{
    
// O(n^2)
    // input: a list P of n points
    // output: distance of the closest pair of points
    double dmin = DBL_MAX;
    
int i, j;
    
for (i = 0; i < n; i++)
        
for( j = i + 1; j < n; j++)
        {
            d = distance(P[i], P[j]);
            
if (d < dmin)
            {
                dmin = d;
            }
        }
    
return dmin;
}

【分治法】

首先划分集合S为SL和SR,使得SL中的每一个点位于SR中每一个点的左边,并且SL和SR中点数相同。分别在SL和SR中解决最近点对问题,得到d1和d2,分别表示SL和SR中的最近点对的距离。令d=min(d1,d2)。如果S中的最近点对(p,q),p在SL并且q在SR中,那么p和q一定在以L为中心的带状区域内,以L-d和L+d为界,如下图所示:

                       

可以证明,对于[L-d,L]区域中的p点,在[L,L+d]中至多需要计算与6个点之间的距离。(证明略)

思路如下

Pseudo Code  
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
 
/*
    version: 1.0
    author: hellogiser
    blog: http://www.cnblogs.com/hellogiser
    date: 2014/7/11
*/

ClosestPair(S)
    
if S.size=1 return DBL_MAX
    
if S.size2=2 return distance(S[0],S[1])
    
//otherwise,do the following
    let L = median(S)
    divide S into SL 
and SR at L
    d1 = CloestPair(SL)
    d2 = CloestPair(SR)
    d12 = CrossPair(SL,SR)
    
return min(d1,d2,d12)

 时间复杂度为T(n)=2T(n/2)+(n/2)*6,可以得到时间复杂度为O(nlgn)。

具体步骤如下:

Step 0  Sort the points by x (list one) and then by y (list two).
 
Step 1 Divide the points given into two subsets S1 and S2 by a vertical line x = m so that half the points lie to the left and half the points lie to the right.
(Note: set m = (x[N/2]+x[N/2+1])/2 so that no points lie on the split line.)
 
Step 2  Find recursively the closest pairs for the left and right subsets.
 
Step 3   Set d = min{d1, d2}
        We can limit our attention to the points in the symmetric vertical strip of width 2d as possible closest pair. Let C1 and C2 be the subsets of points in the left subset S1 and of the right subset S2, respectively, that lie in this vertical strip. The points in C1 and C2 are stored in increasing order of their y coordinates, taken from the second list.
 
Step 4   For every point P(x,y) in C1, we inspect points in C2 that may be closer to P than dThere can be no more than 6 such points (because dd2)!
 伪代码如下:
Pseudo Code 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
 
/*
    version: 1.0
    author: hellogiser
    blog: http://www.cnblogs.com/hellogiser
    date: 2014/7/11
*/

GetCloestPair(pts, n)
    copy pts[
0..n-1] to ptsByX[0..n-1],ptsByY[0..n-1]
    qsort(ptsByX,cmpX)
    qsort(ptsByY,cmpY)
    ClosestPair(ptsByX, ptsByY, n)

ClosestPair(ptsByX, ptsByY, n)
      
// Base cases
    if (n = 1return INT_MAX
    
if (n = 2return distance(ptsByX[0], ptsByX[1])
    
// Divide S into SL SR by line x = xm
    mid = n/2 -1 
    copy ptsByX[
0 . . . mid] into new array XL in x order
    copy ptsByX[mid+
1 . . . n−1] into new array XR
    copy ptsByY[
0 . . . mid] into new array YL in y order
    copy ptsByY[mid+
1 . . . n−1] into new array YR
     
// XL and YL refer to same points, as do XR,YR.
    // Conquer
    d1 = ClosestPair(XL, YL, floor(n/2))
    d2 = ClosestPair(XR, YR, ceil(n/
2))
    
// Combine sub solutions to final solution
    d12 = CrossPair(ptsByX,XL,XR,n,d1,d2);
    
return min(d1,d2,d12)
其中最为重要的是CrossPair步骤。
 Pseudo Code  
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
 
CrossPair(ptsByX,XL,XR,n,d1,d2)
    mid = n/2 -1
    d =  min(d1, d2)
    xm = (ptsByX[mid]+ptsByX[mid+
1])/2
    
//C1: Select points in XL where x>xm-d
    i = mid
    
while(i>=0&&XL[i].x>xm-d)
            add XL[i] to C1
            i = i-
1
    
//C1=XL[i+1..mid]
    //C2: Select points in XR where x<xm+d
    j = mid+1
    
while(j<=n-1&&XR[j].x<xm+d)
            add XR[j] to C2
            j = j+
1
    
//C2=XL[mid+1..j-1]
    // For given Point P in C1, there are at most 6 points in C2 within distance of d
    minDist = DBL_MAX
    
for(i=0;i<C1.length;i++)
        p = C1[i]
        
for(j=0;j<C2.length;j++)
            q = C2[j]
            
// Make sure Q within d*2d rectangel of P(at most 6 Q)
                if(p.y-d<q.y<p.y+d)
                            dist = distance(p,q)
                            
if(minDist>dist) 
                                    minDist = dist
    
return minDist
可以通过left和right下标来表示C1和C2,这样可以进一步优化为
Pseudo Code 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
 
CrossPair(ptsByX,XL,XR,n,d1,d2)
    mid = n/2 -1
    d = min(d1, d2)
    xm = (ptsByX[mid]+ptsByX[mid+
1])/2
    
//C1: Select points in XL where x>xm-d
    i = mid
    
while(i>=0&&XL[i].x>xm-d)
            i = i-
1
    left = i+
1
    
//C1=XL[left..mid]
    //C2: Select points in XR where x<xm+d
    j = mid+1
    
while(j<=n-1&&XR[j].x<xm+d)
            j = j+
1
    right = j-
1
    
//C2=XL[mid+1..right]
    // For given Point P in C1, there are at most 6 points in C2 within distance of d
    minDist = DBL_MAX
    
for(i=left;i<=mid;i++)
        p = XL[i]
        
for(j=mid+1;j<=right;j++)
            q = XR[j]
            
// Make sure Q within d*2d rectangel of P(at most 6 Q)
                if(p.y-d<q.y<p.y+d)
                            dist = distance(p,q)
                            
if(minDist>dist) 
                                    minDist = dist
    
return minDist
【参考】