最近点对问题（算法与数据结构设计）

课题内容和要求

最近点对问题，在二维平面上输入n个点列P。其中任一点pi=（xi，yi），编写程序求出最近的两个点。使用穷举法实现，算法复杂度O(n2)；优化算法，以O(nlog2n)实现这一问题

数据结构说明

Point_X结构体，用于表示点集，n用来存放点集中点的数量，dis_min用来存放求得的最近点对的距离，p数组用于存放所有点的x轴和y轴坐标。

Point结构体，用于存放点的x轴和y轴坐标。

Point_Y结构体，辅助数据结构，用于求分治法中计算跨越左子部分和右子部分存在的解，其中flag用于标记属于左子部分或右子部分，0表示左子部分，1表示右子部分。

算法设计

穷举法求最近点对，利用双重循环嵌套求每个点与其余点的距离，每当有最小的距离出现，记录两个点以及它们的最小距离。

分治法求最近点对，在主程序中对点集X按x轴由小到大排序，以及对点集X的复制点集Y按y轴由小到大排序，在递归函数中不断递归求中点X[mid]并分划为两部分，直至点被分划到少于或等于3个点，直接求它们中的最近点对以及最短距离min。递的部分完成后开始归，从左右子部分当中取最小的min，并在点集Y中直接取属于（X[mid]-min,X[mid]+min）的点存放到点集Y_中，并把属于（X[mid]-min，X[mid]）的左子部分做标记flag=0，把属于（X[mid+1]，X[mid]+min）的右子部分做标记flag=1。接着在点集Y_中找属于左子部分的点，每次找到一个属于左子部分的点，在此基础向上找（因为按y轴排序）四个属于右子部分的点，向下找四个属于右子部分的点，一共8个点，计算左子部分找到点与其对应向上向下的8个点的距离求最短距离。

详细设计

程序包含文件main.c、point.h、point.c

main.c

#include <stdlib.h>
#include <stdio.h>
#include "point.h"
#include <math.h>
int main()
{
    Point_X *X = (Point_X*)malloc(sizeof(Point_X));
    #if 1   // 免输入直接出结果，已给出点集合
    double p[10][2] = {{0,5},{-2,2},{0,1},{2,1},{-5,0},{-0.5,0},{0,0},{0.5,0},{5,0},{0,-5}};
    if(!Init(X,10))
    {
        printf("点集初始化失败\n");
        return 0;
    }
    Point *Y = (Point*)malloc((X->n)*sizeof(Point));
    if(!Y)
    {
        printf("申请空间失败——辅助数据Y\n");
    }
    for(int i = 0;i<(X->n);i++)
    {
        X->p[i].x = p[i][0];
        X->p[i].y = p[i][1];
        Y[i].x = p[i][0];
        Y[i].y = p[i][1];
    }
    #elif 0   // 需要用户输入
    int n = 0;  // 存放点集的数量，用来初始化点集
    printf("请输入点集中点的数量：\n");
    scanf("%d",&n);
    if(!Init(X,n))
    {
        printf("点集初始化失败\n");
        return 0;
    }
    Point *Y = (Point*)malloc((X->n)*sizeof(Point));
    if(!Y)
    {
        printf("申请空间失败——辅助数据Y\n");
    }
    for(int i = 0;i<(X->n);i++)
    {
        printf("请输入第%d个点的x轴：",i+1);
        scanf("%lf",&(X->p[i].x));
        printf("请输入第%d个点的y轴：",i+1);
        scanf("%lf",&(X->p[i].y));
        Y[i].x = X->p[i].x;
        Y[i].y = X->p[i].y;
    }
    getchar();  // 将上一个输入的回车接收
    #endif
    // 点对1，用于记录穷举法得出的最近点对
    Point pair_pp1;
    Point pair_pp2;
    // 点对2，用于记录分治法求得出得最近点对
    Point pair_p1;
    Point pair_p2;
    double min = INFINITY;
    for(int i = 0;i<X->n-1;i++)    // 穷举法求最近点对
    {
        for(int j=i+1 ;j<X->n;j++)
        {
            double dis = distance(X->p[i],X->p[j]);
            if(dis<min) // 当dis小于min值更新最近点对的点以及两点之间的距离
            {
                pair_pp1 = X->p[i];
                pair_pp2 = X->p[j];
                min = dis;
            }
        }
    }
    printf("以下为穷举法求得的最近点对：\n");
    printf("点集中最近的两个点，点1：(%.2lf,%.2lf),点2：(%.2lf,%.2lf),两点距离：%.2lf\n",pair_pp1.x,pair_pp1.y,pair_pp2.x,pair_pp2.y,min);
    printf("---------------------------\n");
    QuickSort_X(X->p,0,X->n-1); // 调用针对x轴排序的快速排序算法对点集X排序
    QuickSort_Y(Y,0,X->n-1);    // 调用针对y轴排序的快速排序算法对点集Y排序
    close_point(X,Y,0,X->n-1,&pair_p1,&pair_p2);    // 调用分治法求最近点对
    printf("以下为分治法求得的最近点对：\n");
    printf("点集中最近的两个点，点1：(%.2lf,%.2lf),点2：(%.2lf,%.2lf),两点距离：%.2lf\n",pair_p1.x,pair_p1.y,pair_p2.x,pair_p2.y,X->dis_min);
    printf("输入回车结束程序。\n");    
    getchar();  // 输入回车结束程序
    free(Y);
    free(X->p);
    free(X);
}
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
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85

point.h

#ifndef __POINT_H__
#define __POINT_H__
typedef struct point
{
    double x;
    double y;
}Point;
typedef struct point_x
{
    Point *p;   // 指向点数组Point的头指针
    double dis_min; // 记录点集中两点最小的距离
    int n;  // 点的数量
}Point_X;
typedef struct point_y  // 辅助数据结构
{
    Point p;
    int flag;   // 表示位，0表示属于左半边，1表示属于右半边
}Point_Y;
int Init(Point_X *X, int n);    // 初始化点集
void Swap(Point *p,int i, int j);   // 交换函数
int Partition_X(Point *p,int left,int right);   // 分化函数，针对x轴
void QuickSort_X(Point *p, int left, int right);    // 快速排序，针对x轴进行排序
int Partition_Y(Point *p,int left,int right);   // 分划函数，针对y轴分划
void QuickSort_Y(Point *p, int left, int right);    // 快速排序，针对y轴排序
double distance(Point p1, Point p2);    // 求点的距离
void close_point(Point_X *X, Point *Y, int left, int right, Point *pair_p1,Point *pair_p2); // 分治法求最近点对
#endif
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

point.c

#include "point.h"
#include <math.h>
#include <stdlib.h>
#include <stdio.h>
/**
 * @brief 初始化点集
 * @param X 指向点集的指针
 * @param n 点的数量
 * @return 0,失败；1，成功
*/
int Init(Point_X *X, int n)
{
    X->p = (Point*)malloc(n*sizeof(Point));
    if(!X->p)
    {
        return 0;
    }
    X->dis_min = INFINITY;
    X->n = n;
    return 1;
}
/**
 * @brief 交换函数
 * @param p 点数组
 * @param i 数组下标i
 * @param j 数组下标j
 * @return void
*/
void Swap(Point *p,int i, int j)
{
    Point tmp;
    if(i==j)
    {
        return;
    }
    tmp = p[i];
    p[i] = p[j];
    p[j] = tmp;
}
/**
 * @brief 分化函数，针对x轴
 * @param p 点数组
 * @param left 左边界
 * @param right 右边界
 * @return j，分划元素下标
*/
int Partition_X(Point *p,int left,int right)
{   
    int i = left;
    int j = right+1;
    do
    {
        do
        {
            i++;
        }while(p[i].x<p[left].x&&i<=right);
        do
        {
            j--;
        } while(p[j].x>p[left].x);
        if(i<j)
        {
            Swap(p,i,j);
        }             
    }while(i<j);
    Swap(p, left, j);
    return j;
}
/**
 * @brief 快速排序，针对x轴进行排序
 * @param p 点数组
 * @param left 左边界
 * @param right 右边界
 * @return void
*/
void QuickSort_X(Point *p, int left, int right)
{
    if(left<right)
    {
        int j = Partition_X(p, left, right);
        QuickSort_X(p, left, j-1);
        QuickSort_X(p, j+1,right);
    }
}
/**
 * @brief 分划函数，针对y轴分划
 * @param p 点数组
 * @param left 左边界
 * @param right 右边界
 * @return j，分划元素下标
*/
int Partition_Y(Point *p,int left,int right)
{   
    int i = left;
    int j = right+1;
    do
    {
        do
        {
            i++;
        }while(p[i].y<p[left].y&&i<=right);
        do
        {
            j--;
        } while(p[j].y>p[left].y);
        if(i<j)
        {
            Swap(p,i,j);
        }             
    }while(i<j);
    Swap(p,left,j);
    return j;
}
/**
 * @brief 快速排序，针对y轴排序
 * @param p 点数组
 * @param left 左边界
 * @param right 右边界
 * @return void
*/
void QuickSort_Y(Point *p, int left, int right)
{
    if(left<right)
    {
        int j = Partition_Y(p, left, right);
        QuickSort_Y(p, left, j-1);
        QuickSort_Y(p, j+1,right);
    }
}
/**
 * @brief 求两点的距离
 * @param p1 点1
 * @param p2 点2
 * @return 返回两点的距离
*/
double distance(Point p1, Point p2)
{
    return sqrt(pow((p1.x-p2.x),2)+pow((p1.y-p2.y),2));
}
/**
 * @brief 分治法求最近点对
 * @param X 对于x轴已经排序好的点集
 * @param Y 对于y轴已经排序好的点集
 * @param left X点集左边界
 * @param right X点集右边界
 * @param pair_1 最近点对的点1
 * @param pair_2 最近点对的点2
 * @return void
*/
void close_point(Point_X *X, Point *Y, int left, int right, Point *pair_p1,Point *pair_p2)
{
    if(right-left<=2)   // 少于或等于三个点直接计算
    {
        if(left==right) // 只有一个点，直接返回
        {
            return;
        }
        for(int i = left;i<right;i++)
        {
            for(int j = i+1;j<right+1;j++)
            {
                double dis = distance(X->p[i],X->p[j]);
                if(X->dis_min>dis)
                {
                    X->dis_min = dis;
                    pair_p1->x = X->p[i].x;
                    pair_p1->y = X->p[i].y;
                    pair_p2->x = X->p[j].x;
                    pair_p2->y = X->p[j].y;
                }
            }
        }
        return;
    }
    int mid = (right+left)/2;   // 分划点
    close_point(X,Y,left,mid,pair_p1,pair_p2);  // 对左边集合进行计算
    close_point(X,Y,mid+1,right,pair_p1,pair_p2);   // 对右边进行计算
    int num_point = 1;  // 用于计算在范围（X[mid]-min，X[mid]+min）内的点的数量
    int mid_ = mid-1;
    while(X->p[mid_].x>(X->p[mid].x-X->dis_min)&&mid_>=0) // 计算在（X[mid]-min，X[mid]）范围的点
    {
        num_point++;
        mid_--;
    }
    mid_ = mid+1;
    while(X->p[mid_].x<X->p[mid].x+X->dis_min&&mid_<=X->n)  // 计算在（X[mid]，X[mid]+min）范围的点
    {
        num_point++;
        mid_++;
    }

    // 在Y点集中提取属于（X[mid]-min,X[mid]+mid）的点，并创建Y_集合存放
    Point_Y *Y_ = (Point_Y*)malloc(num_point*sizeof(Point_Y));
    if(!Y_)
    {
        printf("申请空间失败——Y_\n");
    }
    int j = 0;  // Y_点集下标
    for(int i = 0;i<X->n;i++)
    {
        if(Y[i].x>((X->p[mid].x)-(X->dis_min))&&Y[i].x<=X->p[mid].x)  // 属于（X[mid]-min,X[mid]）存入Y_集合并做标记
        {
            Y_[j].p = Y[i];
            Y_[j].flag = 0;
            j++;
        }
        if(Y[i].x>X->p[mid].x&&Y[i].x<X->p[mid].x+X->dis_min)   // 属于（X[mid],X[mid]+min）存入Y_集合并做标记
        {
            Y_[j].p = Y[i];
            Y_[j].flag = 1;
            j++;
        }
    }
    for(int i = 0;i<num_point;i++) // 计算（X[mid]-min,X[mid]+min）范围内的最小距离
    {
        if(Y_[i].flag==0)   // 取点集合Y_属于（X[mid]-min,X[mid]）的点，并计算每个点与属于（X[mid],X[mid]+min）的最多8个点进行计算
        {
            int four = 0;   // 用于计数
            for(int j = i-1;j>=0;j--)   // 向前找最多4个点进行计算
            {
                if(Y_[j].flag==0)
                {
                    continue;
                }
                four++;
                double dis = distance(Y_[i].p,Y_[j].p);
                if(X->dis_min>dis)
                {
                    X->dis_min = dis;
                    pair_p1->x = X->p[i].x;
                    pair_p1->y = X->p[i].y;
                    pair_p2->x = X->p[j].x;
                    pair_p2->y = X->p[j].y;
                }
                if(four==4) // 最多找4个点，退出该循环
                {
                    four=0;
                    break;
                }
            }
            four=0;
            for(int j = i+1;j<num_point;j++)    // 向后找最多4个点进行计算
            {
                if(Y_[j].flag==0)
                {
                    continue;
                }
                four++;
                double dis = distance(Y_[i].p,Y_[j].p);
                if(X->dis_min>dis)
                {
                    X->dis_min = dis;
                    pair_p1->x = X->p[i].x;
                    pair_p1->y = X->p[i].y;
                    pair_p2->x = X->p[j].x;
                    pair_p2->y = X->p[j].y;
                }
                if(four==4) // 最多找4个点，退出该循环
                {
                    four = 0;
                    break;
                }
            }
four=0;

        }
    }
    free(Y_);
    return;
}
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
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270

测试数据及其结果分析

例子1，程序中自带的例子：
(0,5)(-2,2)(0,1)(2,1)(-5,0)(-0.5,0)(0,0)(0.5,0)(5,0)(0,-5)

运行结果：

例子2，手动输入（程序中默认启用例子1，需要在main.c文件中将条件编译中if后的1改为0，elif后的0改为1）：
(0,0)(1,1)(0.25,0)(1.1,1)(-2,1)(2,2)(3,3)(-1,-1)

运行结果：