平面最近点对

最近点对(Closest Pair of Points)问题是计算几何的问题:在平面上给出 n 个点,求出距离最短的那对点

该问题可以用 分治法 来解决,优化后时间复杂度可以为O(nlogn)

下面是算法过程:

1.将所有的点分别通过x坐标排序存入Px,通过y坐标排序存入Py

2.取Px中间的点 m,以m为中间的点将Px分成两半

3.递归求解这两半的最短距离解假设为 d

4.创建数组 strip[] 存入所有点中到 m 的距离至少为 d 的点

5.找出 strip[] 中最小的点 d1

6.返回 d1 d 中更小的那个。

从简单的说,就是分成两半,分别求出最小值,难就难在要处理中间左右两半交叉的那部分。

下面是例子

#include <stdio.h>
#include <string.h>
#include <math.h>
#include <stdint.h>
#include <float.h>

#define dl_min(x,y) ((x) < (y) ? (x) : (y));

typedef struct {
    int x,y;
} point_t;

//求两点的距离,即求垂直三角形的最长边
float dist(point_t *p1, point_t *p2) 
{
    int xb = (p1->x - p2->x);
    int yb = (p1->y - p2->y);
    return sqrt( xb*xb + yb*yb );
} 

float brute_force(point_t *p_x, int n)
{
    int i,j;
    float dt;
    float min = FLT_MAX; 
    for (i = 0; i < n; ++i) {
        for (j = i+1; j < n; ++j){
            dt = dist(&p_x[i], &p_x[j]);
            if (dt < min)
                min = dt; 
        }
    }
    return min;
}

int compare_x(void *p1, void *p2)
{
    point_t     *pp1, *pp2;

    pp1 = p1;
    pp2 = p2;

    if (pp1->x < pp2->x)
        return -1;
    else if (pp1->x > pp2->x)
        return 1;

    return 0;
}

int compare_y(void *p1, void *p2)
{
    point_t     *pp1, *pp2;

    pp1 = p1;
    pp2 = p2;

    if (pp1->y < pp2->y)
        return -1;
    else if (pp1->y > pp2->y)
        return 1;

    return 0;
}

float strip_closest(point_t *strip, int size, float d) 
{ 
    int i,j;
    float min = d;  // Initialize the minimum distance as d 
    float tmp;
  
    // 一个个循环比较,这个已经被证明最多循环6次
    for (i = 0; i < size; ++i) 
        for (j = i+1; j < size && (strip[j].y - strip[i].y) < min; ++j) {
            tmp = dist(&strip[i],&strip[j]);
            if (tmp < min) 
                min = tmp;
        }
  
    return min; 
}

void dump_p(point_t *p, int n)
{
    int i;

    for (i = 0; i < n; i++) {
        printf("(%d,%d) ", p[i].x, p[i].y);
    }

    printf("\n");
}

float do_closest_p(point_t *p_x, point_t *p_y, int n)
{
    int i;

    //3以内的点直接遍历比较
    if (n <= 3)
        return brute_force(p_x, n);

    int mid = n / 2;
    point_t *mid_p = &p_x[mid];

    point_t pyl[mid+1];
    point_t pyr[n-mid-1];

    //这里有个技巧,遍历按照y轴排序好的p_y,那么得到的两半数据都是按照y排序好的
    int li = 0,ri = 0;
    for (i = 0; i < n; i++) {
        //控制左边最多存储 mid+1 个点
        if (p_y[i].x <= mid_p->x) {
            if (li == mid)
                break;
            pyl[li++] = p_y[i];
        } else
            pyr[ri++] = p_y[i];
    }

    //剩余的存入右边
    for (i = 0; i < n; i++)
        pyr[ri++] = p_y[i];


    //递归求取左右边最小距离
    float dl = do_closest_p(p_x, pyl, mid);
    float dr = do_closest_p(p_x+mid, pyr, n - mid);
    float d = dl_min(dl, dr);

    //处理中间的情况
    point_t strip[n];
    int     j = 0;
    for (i = 0; i < n; i++)
        if (abs(p_y[i].x - mid_p->x) < d)
            strip[j++] = p_y[i];


    return dl_min(d, strip_closest(strip, j, d));
}

float closest_p(point_t *p, int n)
{
    point_t p_x[n];
    point_t p_y[n];

    memcpy(p_x, p, sizeof(point_t) * n);
    memcpy(p_y, p, sizeof(point_t) * n);

    qsort(p_x, n, sizeof(point_t), compare_x);
    qsort(p_y, n, sizeof(point_t), compare_y);

    dump_p(p_x, n);
    dump_p(p_y, n);

    return do_closest_p(p_x, p_y, n);
}

int main(int argc, char **argv)
{
    point_t p[] = {{1,3},{2,1},{3,4},{4,7},{6,2},{7,5}};
    int n = sizeof(p) / sizeof(point_t);

    printf("distance is %f\n", closest_p(p, n));

    return 0;
}
[root@izbp1irxwqt7ei21awv6wvz ccc]# gcc main.c -lm
[root@izbp1irxwqt7ei21awv6wvz ccc]# ./a.out 
(1,3) (2,1) (3,4) (4,7) (6,2) (7,5) 
(2,1) (6,2) (1,3) (3,4) (7,5) (4,7) 
distance is 2.236068

时间复杂度为 O(nlogn)


参考:

wiki-Closest_pair_of_points_problem

geeksforgeeks-closest-pair-of-points-onlogn-implementation/

上一篇: 分治法
下一篇: 图算法
作者邮箱: 203328517@qq.com