最近对问题和凸包问题

最近对问题

最近点对问题要求在一个包含n个点的集合中，找出距离最近的两个点

蛮力法

很显然，蛮力法应该是这样，分别计算每一对点之间的距离，然后找出距离最小的那一对。时间复杂度为O(n²)

double BruteForceCloestPoints(Node a[])
{
	int mindist = inf;
	for(int i = 0;i<n-1;++i){
		for(int j = i+1;j<n;++j){
			double curdist = sqrt(pow(a[i].x-a[j].x,2)+pow(a[i].y-a[j].y,2));
			mindist = min(mindist,curdist);
		}
	}
	return mindist;
}

分治法

我们将点按照x轴排序后，就会发现，最左边的那些结点和最右边的那些结点显然不是最短距离，所以最短距离肯定存在于某个区间上面，且子问题没有交集，所以我们可以使用分治法，将区间一分为二，然后求解子问题。
当然有可能存在一个点在中线左边，一个点在中线右边，这时我们求完子问题后得到最小距离还有再求位于两端的距离小于d的点的最小距离与子问题得出的最小距离比较。当n<=3的时候，返回由蛮力算法求出的最小距离

typedef struct 
{
	int x;
	int y;
}point;

double distance(const point& a, const point& b)
{
	return sqrt(pow(a.x - b.x, 2) + pow(a.y - b.y, 2));
}

double merge(point a[], int low, int high,double d)
{
	double mindist = d;
	vector<point> left, right;
	int m = a[(low + high) / 2].x;
	for (int i = low; i <= high; ++i){
		if (a[i].x - m < 0 && a[i].x - m > -d)
			left.push_back(a[i]);
		if (a[i].x - m > 0 && a[i].x - m < d)
			right.push_back(a[i]);
	}
	for (auto& p : left) {
		for (auto& q : right) {
			mindist = min(mindist, distance(p,q));
		}
	}
	return mindist;
}

double EfficientClosestPair(point a[], int low, int high)
{
	if (high - low < 2)
		return distance(a[low],a[low + 1]);
	if (high - low < 3)
		return min(distance(a[low],a[low+1]), min(distance(a[low], a[low + 2]), distance(a[low+1], a[low + 2])));
	int mid = (low + high) / 2;
	double d1 = EfficientClosestPair(a, low, mid);
	double d2 = EfficientClosestPair(a, mid + 1, high);
	double d = min(d1, d2);
	return merge(a, low, high,d);
}

凸包问题

凸包问题被视为计算几何中最重要的问题之一，定义如下：

在平面或者高维空间的一个给定点集合中寻找凸包

该问题之所以如此突出，是由于许多各种各样的应用要么本身就是凸包问题，要么其中一部分需要按照凸包问题来解决。大多数此类应用都是基于这样一个原理：凸包能方便地提供目标形状或给定数据集的一个近似。

凸集合的概念：

对于平面上的一个点集合（有限的或无限的），如果以集合中任意两点p和q为端点的线段上所有点都属于该集合，我们说这个集合是凸的。

下面这个凸包的正式定义可以应用于任意集合，包括那些正好位于一条直线上的点的集合。

一个点集合S的凸包是包含S的最小凸集合（“最小”意指，S的凸包一定是所有包含S的凸集合的子集）

如果S是凸的，它的凸包很明显是它本身。如果S是两个点组成的集合，它的凸包是连接着两个点的线段。如果S是由三个不同线的点组成的集合，它的凸包是以这三个点为顶点的三角形；如果三点同线，凸包是以距离最远的两个点为端点的线段。对于更大集合的凸包的例子，请见下图

凸包问题是为一个有n个点的集合构造凸包问题。为了解决该问题，需要找出某个点，它们将作为这个集合的凸多边形的顶点。

贪心法

对于一个n个点集合中的两个点p_i和p_j，当且仅当该集合中的其他点都位于穿过这两点的直线的同一边时，它们的连线是该集合凸包边界的一部分。对每一对点都做一遍检验之后，满足条件的线段构成了该凸包的边界

为了实现这个算法，需要用到一些解析集合的基本知识。首先，在坐标平面上穿过两个点(x₁,y₁)，(x₂,y₂)的直线是由下列方程定义的：
ax+by = c
其中 a = y₂ - y₁,b = x₁-x₂ ,c = x₁y₂ - y₁x₂。
其次，这样一根直线把平面分为两个半平面：其中一个半平面中的点都满足ax+by>c，而另一个半平面中的点都满足ax+by<c。因此，为了检验某些点是否位于这条直线的同一边，只需要把每个点代入ax+by-c,检验这个表达式的符号是否相同。

typedef struct
{
	int x;
	int y;
}point;
typedef struct
{	
	point one;
	point two;
}Edge;

void BruteForceClosure(const vector<point>& invec,vector<point>& outvec)
{
	int a, b, c;
	int last;
	bool first = true;
	bool can = true;
	for (int i = 0; i < invec.size() - 1; ++i) {
		for (int j = i + 1; j < invec.size(); ++j) {
			a = invec[j].y - invec[i].y;
			b = invec[i].x - invec[j].x;
			c = invec[i].x * invec[j].y - invec[i].y * invec[j].x;
			for (int k = 0; k < invec.size(); ++k) {
				if (k != i && k != j) {
					int ck = a * invec[k].x + b * invec[k].y - c;
					if (first) {
						first = false;
						last = ck;
					}
					else {
						if (last * ck < 0) {
							can = false;
							break;
						}
					}
				}
			}
			if (can) {
				Edge e;
				e.one = invec[i];
				e.two = invec[j];
				outvec.emplace_back(e);
			}
		}
	}
}

分治法（快包算法）

假设集合S是平面上上n>1个点构成的。我们还假设这些点按照他们的x轴坐标升序排列的，如果x轴坐标相同，则按照y轴坐标升序排列。不难证明这样一个几何学上的明显事实，最左面的点p₁和最右面的点p_n一定是该集合的凸包顶点。这条直线把点分为两个集合：S1是位于直线左侧或在直线上的点构成的集合，S2是位于直线右侧或在直线上的点构成的集合。

S的凸包的边界是由下面两个多角形链条构成的：一条“上”边界和一条“下”边界。“上”边界称为上包，是一系列线段的序列，“下”边界称为下包，是一系列线段的序列

快包算法：

到现在我们还需要具备一个几何知识， 即：
如果我们有三个点A(x1,y1),B(x2,y2), C(x3,y3)

我们可以用以下这个行列式求面积 (也就是我们常说的"叉积")

这个行列式对于平面上任意三角形, 求解面积都很方便, 所得结果是三角形ABC面积的两倍 (可以回忆高中知识~)

而且, 当A->B->C为顺时针顺序时, 该值为正, 反之则为负

struct point
{
	int x;
	int y;
}p[nmax];

int visited[nmax];

int Djudge(point a1, point a2, point a3)
{
	int calculate = a1.x * a2.y + a3.x * a1.y + a2.x * a3.y - a3.x * a2.y - a2.x * a1.y - a1.x * a3.y;
	return calculate;
}

bool cmp(const point& a, const point& b)
{
	if (a.x != b.x)
		return a.x < b.x;
	else
		return a.y < b.y;
}

void DealConvex(int first, int last)
{
	int max = 0, index = -1;
	//求上包
	if (first < last) {
		//注意两端没有必要计算
		for (int i = first + 1; i < last; ++i) {
			int calcu = Djudge(p[first], p[i], p[last]);
			if (calcu == 0) visited[i] = 1;
			if (calcu > nmax) {
				max = calcu;
				index = i;
			}
		}
	}
	//求下包
	else {
		for (int i = first - 1; i > last; --i) {
			int calcu = Djudge(p[first], p[i], p[last]);
			if (calcu == 0) visited[i] = 1;
			if (calcu > nmax) {
				max = calcu;
				index = i;
			}
		}
	}
	if (index = -1) {
		visited[index] = 1;
		DealConvex(first, index);
		DealConvex(index, last);
	}
}

void PrintPoint(int n) {
	int mark[nmax];
	int t = 0;
	point ans[nmax];
	for (int i = 0; i < n; i++)
	{
		if (visited[i] == 1)
		{
			ans[t].x = p[i].x;
			ans[t].y = p[i].y;
			t++;
		}
	}
	//顺时针输出
	mark[0] = mark[t - 1] = 1; //数组mark避免重复检查降低效率
	for (int i = 1; i < t - 1; i++)
	{
		mark[i] = 0;
	}
	cout << ans[0].x << " " << ans[0].y << endl;
	for (int i = 1; i < t - 1; i++)
	{
		int d = Djudge(ans[0], ans[t - 1], ans[i]);
		if (d >= 0)
		{
			cout << ans[i].x << " " << ans[i].y << endl;
			mark[i] = 1;
		}
	}
	cout << ans[t - 1].x << " " << ans[t - 1].y << endl;
	for (int i = 1; i < t; i++)
	{
		if (mark[i] != 1)
		{
			int d = Djudge(ans[0], ans[t - 1], ans[i]);
			if (d < 0)
			{
				cout << ans[i].x << " " << ans[i].y << endl;
			}
		}
	}
}