最近点对问题(一维、二维)

（最近点对解可能不唯一，只求出一对）

1、一维空间最近点对问题

1.1 蛮力法

• 用O(nlogn)时间对他们进行排序
• 走过排序好的表，计算每一个点到跟在它后面的点的距离，容易求出这些距离的最小值，时间复杂度为O(n-1)
• 故T(n)=O(nlogn)

1.2 分治法

• 将所有点从小到大排序
• 取排序好的点的中间点m，以它为分界线将点集合S划分为两个子集S1、S2（如下图所示）。
• 递归地在S1、S2上找出其最接近点对 {p1， p2} 和 {q1， q2},并设d=min{|p1-p2|, |q1-q2|},S中最接近点对或者是 {p1， p2} ，或者是 {q1， q2}，或者是某个{p3，q3}，其中p3在S1集合里，q3在S2集合里
  
• 如果S最接近点对是{p3， q3 }，即 |p3-q3|<d, 则p3和q3两者与m的距离不超过d，即p3-m<d，q3-m<d，即p3在 (m-d, m]区间内，q3在(m, m+d]区间内
• 由于在S1中，每个长度为d的半闭区间之多包含一个点（否则必有两点距离小于d），并且m是S1和S2的分割点，因此 (m-d, m]中至多包含 S中的一个点。由图可以看出，如果(m-d, m]中有S中的点，则此点就是S1中最大点。同理，如果(m, m+d]中有S中的点，则此点就是S2中最小点
• 因此，可以用线性时间就能找到区间(m, m+d]和(m, m+d]中所有点，即p3和q3。从而用线性时间就可以将S1的解和S2的解合并成S的解

2、 二维空间最接近点对问题

2.1 蛮力法

def ClosePoints(Points_X, leftindex, rightindex):
    mindist = np.inf
    close_points = []
    for i in range(leftindex, rightindex+1):
        for j in range(i+1, rightindex+1):
            distance = Distance(Points_X[i], Points_X[j])
            if distance < mindist:
                # 每次保存最近两个点的信息之前先清空列表
                close_points.clear()
                mindist = distance
                # 保存最近两个点的信息
                close_points.append(Points_X[i])
                close_points.append(Points_X[j])
    return mindist, close_points

2.2 分治法

2.2.1 求解思路

对于给定的点集a[0…n-1]，采用分治法求最近点对距离的步骤如下：

• 对a中的所有点按x坐标从小到大排列，将a中点复制到b中，对b中的所有点按y坐标从小到大排序。设a中最近点对距离d
• 如果a中点数少于4，则采用蛮力法直接计算各点的最近距离d
• 求出a中间位置的点a[midindex]，以此位置画一条中轴线l（对应的x坐标为a[midindex].x），将a所有点分割为点数大致相同的两个子集：左部分包含a[0…midindex]的点，右部分包含a[midindex+1…n-1]的点。同样将b中的点相应分为两部分leftb和rightb，左部分称为S1（含a[0…midindex]和leftb），右部分为S2（含a[midindex+1…n-1]和rightb），如下图所示。
• 递归调用求出S1中点集的最近点对的距离为d1，递归调用求出S2中点集的最近点对的距离为d2，求出当前最近点对的距离为d=MIN(d1，d2)

【注意】

• 显然S1和S2中任意点对之间的距离小于或等于d，但S1、S2交界的垂直带形区（由所有与中轴线的x坐标值相差不超过d的点构成）中的点对之间的距离可能小于d。将b中所有落在垂直带形区的带你复制到b1中，对于b1中的任一点p，仅需考虑紧随p后的7个点，计算出从p到这7个点的距离，并和d进行比较，将最小的距离存放在d中，最后求得的d即为a中所有点的最近点对问题。
• 对于b1中的点p，为什么只需要考虑紧随p后的7个点呢？如果pl属于Pl，pr属于Pr，且pl和pr的距离小于d，则他们必定位于以 l 为中轴线的d*2d的矩形区域内，如下图所示，该区域最多只有8个点【好像是6个点？？？】（左右阴影正方形中最多有4个点，否则他们距离小于d，与Pl、Pr中所有点的最小距离大于或等于d矛盾）。所以为了求Pl和Pr中点之间的最小距离。只需考虑每个点p之后的7个点就可以了。

2.2.2 例子

以此图为例：

• 对a[0…9]按x坐标从小到大排序结果为a9(1，6)，a1(3，7)，a3(3，4)，a8(3，0)，a0(4，10)，a4(5，6)，a5(5，4)，a6(6，3)，a7(8，1)，a2(9，7)
• 对b[0…9]按y坐标从小到大排序结果为a8(3，0)，a7(8，1)，a6(6，3)，a3(3，4)，a5(5，4)，a9(1，6)，a4(5，6)，a1(3，7)，a2(9，7)，a0(4，10)
• 取中间位置midindex=（0+9）/ 2 = 4，左部分为{a9，a1，a3，a8，a0}，右部分为{a4，a5，a6，a7，a2}
• 对于整个序列的左部分{a9，a1，a3，a8，a0}，取中间位置midindex=2；它左部分为{a9，a1，a3}，点数少于4，采用蛮力法求出a9和a1的最小距离为d1=2.23607。它右部分为{a8，a0}，点数少于4，采用蛮力法求出a8和a0的最小距离为d2=10.0499。则d=MIN(d1，d2)=2.23607。中间部分为{a9，a1，a3，a8}，通过对应的b1点集求出其中最近点对距离，并和d进行比较，最终结果为d=2.23607，最近点对为a9、a1
• 对于整个序列的右部分{a4，a5，a6，a7，a2}，取中间位置midindex=7；它左部分为{a4，a5，a6}，点数少于4，采用蛮力法求出a5和a6的最小距离为d1=1.41421。它右部分为{a7，a2}，点数少于4，采用蛮力法求出a7和a2的最小距离为d2=6.08276。则d=MIN(d1，d2)=1.41421。中间部分为空。最终结果为d=1.41421，最近点对为a5、a6
• 整个序列的中间部分为{a3，a0，a5，a2}，通过对应的b1点集求出其中最近点对距离，并和d进行比较，最终结果为d=1.41421，最近点对为a5、a6

2.2.3 代码

C++：
close_points.app:

//求解最近点对问题的算法
#include "Fundament.cpp"
#define min(x,y) ((x)<(y)?(x):(y))
//以下为求最近点对的两个算法和主函数
double ClosestPoints(vector<Point> a,int leftindex,int rightindex)
//蛮力法求a[leftindex..rightindex]中的最近点对距离
{
	int i,j;
	double d,mindist=INF;
    for (i=leftindex;i<=rightindex;i++)
        for (j=i+1;j<=rightindex;j++)
        {
			d=Distance(a[i],a[j]);
			if (d<mindist)
				mindist=d;
		}
	return mindist;
}
void display(vector<Point> a,int leftindex,int rightindex)	//输出a中的部分元素，用于调试
{
	int i;
	for (i=leftindex;i<=rightindex;i++)
		a[i].disp();
	printf("\n");
}
//分治递归求点对距离
bool pointxcmp(Point &p1,Point &p2)		//用于点按x坐标递增排序
{
	return p1.x<p2.x;
}
bool pointycmp(Point &p1,Point &p2)		//用于点按y坐标递增排序
{
	return p1.y<p2.y;
}
double ClosestPoints11(vector<Point> &a,vector<Point> b,int leftindex,int rightindex)
//递归求a[leftindex..rightindex]中最近点对的距离
{
	vector<Point> leftb,rightb,b1;
	int i,j,midindex;
	double d1,d2,d3=INF,d;
	if ((rightindex-leftindex+1)<=3)		//少于4个点,直接用蛮力法求解
	{
		d=ClosestPoints(a,leftindex,rightindex);
		return d;
	}
    midindex=(leftindex+rightindex)/2;	//求中间位置
	for (i=0;i<b.size();i++)			//将b中点集分为左右两部分
		if (b[i].x<a[midindex].x)
			leftb.push_back(b[i]);
		else
			rightb.push_back(b[i]);
    d1=ClosestPoints11(a,leftb,leftindex,midindex);
	d2=ClosestPoints11(a,rightb,midindex+1,rightindex);
	d=min(d1,d2);	//当前最小距离d=MIN(d1,d2)
    //求中间部分点对的最小距离
	for (i=0;i<b.size();i++)	//将b中间宽度为2*d的带状区域内的子集复制到b1中
		if (fabs(b[i].x-a[midindex].x)<=d)
			b1.push_back(b[i]);
	double tmpd3;
	for (i=0;i<b1.size();i++)	//求b1中最近点对
		for (j=i+1;j<b1.size();j++)
		{
			if (b1[j].y-b1[i].y>=d) break;
			tmpd3=Distance(b1[i],b1[j]);
			if (tmpd3<d3)
				d3=tmpd3;
		}
	d=min(d,d3);
	return d;
}
double ClosestPoints1(vector<Point> &a,int leftindex,int rightindex)
//求a[leftindex..rightindex]中最近点对的距离
{
	int i;
	vector<Point> b;
	vector<Point>::iterator it;
	printf("排序前:\n");
	for (it=a.begin();it!=a.end();it++)
		(*it).disp();
	printf("\n");
    sort(a.begin(),a.end(),pointxcmp);	//按x坐标从小到大排序
	printf("按x坐标排序后:\n");
	for (it=a.begin();it!=a.end();it++)
		(*it).disp();
	printf("\n");

	for (i=0;i<a.size();i++)			//将a中点集复制到b中
		b.push_back(a[i]);

    sort(b.begin(),b.end(),pointycmp);	//按y坐标从小到大排序
	printf("按y坐标排序后:\n");
	for (it=b.begin();it!=b.end();it++)
		(*it).disp();
	printf("\n");
	return ClosestPoints11(a,b,0,a.size()-1);
}
//求最近点对的主函数
int main()
{
	vector<Point> a;
	a.push_back(Point(4,10));
	a.push_back(Point(3,7));
	a.push_back(Point(9,7));
	a.push_back(Point(3,4));
	a.push_back(Point(5,6));
	a.push_back(Point(5,4));
	a.push_back(Point(6,3));
	a.push_back(Point(8,1));
	a.push_back(Point(3,0));
	a.push_back(Point(1,6));
	double mindist=ClosestPoints1(a,0,a.size()-1);
	printf("最近距离=%g\n",mindist);
}

Fundament.cpp：

#include <algorithm>
#include <vector>
using namespace std;
#include <stdio.h>
#include <math.h>
#define INF 32767			
#define MAXN 20

int sign(double x)			//返回数值x的符号
{
	if (x>0) return 1;
	else if (x<0) return -1;
	else return 0;
}
double sqr(double x)		//求x的平方
{
	return x*x;
}

class Point					//点类
{
public:
	double x;				//行坐标
	double y; 				//列坐标
	Point() {}						//默认构造函数
	Point(double x1,double y1)		//重载构造函数
	{	x=x1;
		y=y1;
	}
	void disp()
	{	printf("(%g,%g) ",x,y); }
	friend double Distance(Point p1,Point p2);			//求两个点的距离
};

double Distance(Point p1,Point p2)		//两个点之间的距离
{
	return sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y));
}

Python:

from operator import itemgetter
import numpy as np

# 求解两点距离
def Distance(a_point, b_point):
    return np.sqrt((a_point[0]-b_point[0])**2 + (a_point[1]-b_point[1])**2)

# 蛮力法求解最近两个点的距离
def ClosePoints(Points_X, leftindex, rightindex):
    mindist = np.inf
    close_points = []
    for i in range(leftindex, rightindex+1):
        for j in range(i+1, rightindex+1):
            distance = Distance(Points_X[i], Points_X[j])
            if distance < mindist:
                # 每次保存最近两个点的信息之前先清空列表
                close_points.clear()
                mindist = distance
                # 保存最近两个点的信息
                close_points.append(Points_X[i])
                close_points.append(Points_X[j])
    return mindist, close_points


def ClosePoints11(Points_X, Point_Y, leftindex, rightindex):
    left_Points_Y = []
    right_Point_Y = []
    Point_Y_temp = []

    if (rightindex-leftindex+1) <= 3:   # 如果少于4个点，直接用蛮力法求解
        distance, close_points = ClosePoints(Points_X, leftindex, rightindex)
        return distance, close_points

    midindex = int((leftindex+rightindex)/2)    # 求中间位置
    for i in range(len(Point_Y)):       # 将Points_y分为左右部分
        if Point_Y[i][0] < Points_X[midindex][0]:
            left_Points_Y.append(Point_Y[i])
        else:
            right_Point_Y.append(Point_Y[i])
    distance1, close_points1 = ClosePoints11(Points_X, left_Points_Y, leftindex, midindex)
    distance2, close_points2 = ClosePoints11(Points_X, right_Point_Y, midindex+1, rightindex)
    if distance1 < distance2:
        min_distance = distance1
        close_points = close_points1
    else:
        min_distance = distance2
        close_points = close_points2

    # 求中间部分最小距离
    for i in range(len(Point_Y)):   # b中间宽度为2*d的带状区域内的子集复制到Point_Y_temp中
        if np.fabs(Point_Y[i][0]-Points_X[midindex][0]) <= min_distance:
            Point_Y_temp.append(Point_Y[i])

    # 求Point_Y_temp中最近点
    distance3 = np.inf
    close_points3 = []
    for i in range(len(Point_Y_temp)):
        for j in range(i+1, len(Point_Y_temp)):
            if (Point_Y_temp[j][1]-Point_Y_temp[i][1]) >= min_distance:
                break
            temp_dis = Distance(Point_Y_temp[i], Point_Y_temp[j])
            if temp_dis < distance3:
                distance3 = temp_dis
                close_points3.append(Point_Y_temp[i])
                close_points3.append(Point_Y_temp[j])
    if min_distance > distance3:
        d = distance3
        close_points = close_points3
    else:
        d = min_distance
    return d, close_points


def ClosePoints1(Point, leftindex, rightindex):
    print('排序前：', Point)
    Point_X = Point
    Point_X = sorted(Point_X)
    print('按X坐标排序后：', Point_X)
    Point_Y = Point
    Point_Y = sorted(Point_Y, key=itemgetter(1))
    print('按Y坐标排序后：', Point_Y)
    return ClosePoints11(Point_X, Point_Y, leftindex, rightindex)

if __name__ == '__main__':
    Point = [(4, 10), (3, 7), (9, 7), (3, 4), (5, 6), (5, 4), (6, 3), (8, 1), (3, 0), (1, 6)]
    mindist, close_points = ClosePoints1(Point, 0, len(Point)-1)
    print(close_points)
    print(mindist)