计算几何相关问题

1.最近点

// 最近点对.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include <iostream>
#include <vector>
#include <algorithm>
#include <WINDOWS.H>
#include <SET>
using namespace std;
/************************************************************************/
/*最近点  基本的分治思想 
                                                                     */
/************************************************************************/

struct Point{int x; int y;};

typedef vector<Point> PointsVec;

//#define min(a,b) a<b?a:b 

void PrintVector(const vector<Point>& src)
{
	vector<Point>::const_iterator it;
	for (it=src.begin();it!=src.end();it++)
	{
		printf("%d  %d\n",it->x, it->y);
	}
}

bool operator==(const Point& p1, const Point& p2){
	return p1.x==p2.x && p1.y==p2.y;
}

bool operator<(const Point& p1, const Point& p2){
	return p1.x<p2.x;
}
bool operator>(const Point& p1, const Point& p2){
	return p1.x>p2.x;
}

int distance(Point p1, Point p2){
	return (p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y);
}

int merge(PointsVec Sy, int miniDist, int midValue){
	PointsVec::iterator it;
	for (it=Sy.begin();it!=Sy.end();)
	{
		if(it->x <= midValue+miniDist && it->x >= midValue-miniDist)
			it++;
		else
			it = Sy.erase(it);
	}

	for (int i=0; i<Sy.size(); i++)
	{
		for (int j=i+1; j<=i+7&&j<Sy.size();j++)
		{
			int tmp = distance(Sy[i],Sy[j]);
			if(tmp<miniDist)
					miniDist = tmp;
		}
	}
	return miniDist;
}

int ClosestPoints(PointsVec Sx, PointsVec Sy){
	int miniDist = INT_MAX;
	if(Sx.size()<=3){
		for (int i=0; i<Sx.size(); i++)
		{
			for (int j=i+1; j<Sx.size(); j++)
			{
				int tmp = distance(Sx[i],Sx[j]);
				if(tmp<miniDist)
					miniDist = tmp;
			}
		}
	}else{
		int midIndex = (Sx.size()-1)/2;
		int midValue = Sx[midIndex].x;
		PointsVec Sxl, Sxr, Syl, Syr;
		for (int i=0; i<Sx.size(); i++){
			if(Sx[i].x<=midValue)
				Sxl.push_back(Sx[i]);
			else
				Sxr.push_back(Sx[i]);

			if(Sy[i].x<=midValue)
				Syl.push_back(Sy[i]);
			else
				Syr.push_back(Sy[i]);
		}
		int minL = ClosestPoints(Sxl,Syl);
		int minR = ClosestPoints(Sxr,Syr);
		miniDist = min(minL,minR);
		miniDist = merge(Sy,miniDist,midValue);
	}
	return miniDist;
}

int NativeClosePoints(PointsVec srcPoints){
	int miniDist = INT_MAX;
	for (int i=0; i<srcPoints.size(); i++)
	{
		for (int j=i+1; j<srcPoints.size(); j++)
		{
			int tmp = distance(srcPoints[i],srcPoints[j]);
			if(tmp<miniDist)
				miniDist = tmp;
		}
	}
	return miniDist;
}

bool cmp1(Point a, Point b){
	return a.x < b.x;
}

bool cmp2(Point a, Point b){
	return a.y < b.y;
}

int main(int argc, char* argv[])
{
	PointsVec srcPoint,Sy;

	srand(GetTickCount());

	set<Point> points;
	for (int i=0; i<100; i++)
	{
		Point pt;
		pt.x = rand()%10;
		pt.y = rand()%10;
		points.insert(pt);//.push_back(pt);
	}
	for (set<Point>::iterator it=points.begin(); it!=points.end(); it++)
	{
		srcPoint.push_back(*it);
	}

	Sy = srcPoint;
	sort(srcPoint.begin(), srcPoint.end(), cmp1);
	//PrintVector(srcPoint);
	sort(Sy.begin(), Sy.end(), cmp2);
	//PrintVector(Sy);
	cout << ClosestPoints(srcPoint,Sy) << endl;
	cout << NativeClosePoints(srcPoint) << endl;

	return 0;
}

2.凸包及其直径

// 凸包.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include <iostream>
#include <VECTOR>
#include <algorithm>

using namespace std;

struct Point{int x;int y;};

bool operator==(const Point& a, const Point& b){
	return a.x==b.x && a.y ==b.y;
}

bool compare(const Point& p1, const Point& p2)
{
	int temp = p1.x*p2.y-p1.y*p2.x;
	if(temp<0)
		return false;
	else if(temp==0 && p1.x*p1.x+p1.y*p1.y>p2.x*p2.x+p2.y*p2.y)
		return false;
	else
		return true;

}

void PrintVector(const vector<Point>& src)
{
	vector<Point>::const_iterator it;
	for (it=src.begin();it!=src.end();it++)
	{
		printf("%d  %d\n",it->x, it->y);
	}
}

bool isLeft(Point p1, Point p2, Point p3)
{
	Point pt1, pt2;
	pt1.x = p1.x - p2.x;
	pt1.y = p1.y - p2.y;

	pt2.x = p3.x - p2.x;
	pt2.y = p3.y - p2.y;
	
	int temp = pt2.x*pt1.y - pt1.x*pt2.y;
	if(temp > 0)
		return true;
	else
		return false;
}

/************************************************************************/
/*计算凸包：主要解决两个问题
1.找到基准点之后，要按照每个点和基准点组成的向量和X轴横坐标的夹角，升序排序
这里判断的依据是例如pi,pj和p0    则(pi-p0)*(pj-p0)>0   则pi的夹角小于pj   (pi-p0)*(pj-p0)<0   则pi的夹角大于pj
(pi-p0)*(pj-p0)=0 说明三点在一条直线上，只取离基准点最远的那个点，其他点可以忽略
2.第二个难点就是判断是否左转   例如pi,pi+1和pi+2三个点
(p[i+1])-p[i])*(p[i-1]-p[i])>0   加入p[i+1]
否则 弹出p[i]   继续回溯
                                                                     */
/************************************************************************/
void CalcConvexHull(vector<Point>& srcPoints, vector<Point>& desPoints){
	if(srcPoints.size()<3)
		return;
	Point ptBase = srcPoints.front();
	vector<Point>::iterator it;
	for (it=srcPoints.begin(); it!=srcPoints.end();it++)
	{
		if(it->y < ptBase.y)
			ptBase = *it;
		else if(it->y==ptBase.y && it->x<ptBase.x)
			ptBase = *it;
	}

	for (it=srcPoints.begin(); it!=srcPoints.end();)
	{
		it->x -= ptBase.x;
		it->y -= ptBase.y;
		it++;
	}

	sort(srcPoints.begin(),srcPoints.end(),compare);
	PrintVector(srcPoints);

	it = srcPoints.begin();
	desPoints.push_back(*it);
	it++;
	desPoints.push_back(*it);
	it++;
	desPoints.push_back(*it);
	it++;
	for (;it!=srcPoints.end();it++)
	{
		Point p1, p2, p3;
		p1 = desPoints.at(desPoints.size()-2);
		p2 = desPoints.at(desPoints.size()-1);
		p3 = *it;
		while(!isLeft(p1,p2,p3)){
			desPoints.pop_back();
			p1 = desPoints.at(desPoints.size()-2);
			p2 = desPoints.at(desPoints.size()-1);
		}
		desPoints.push_back(p3);
	}
	printf("******************************\n");
	PrintVector(desPoints);
}

int cross(Point p1, Point p2, Point p3){
	Point pt1, pt2;
	pt1.x = p1.x - p2.x;
	pt1.y = p1.y - p2.y;
	
	pt2.x = p3.x - p2.x;
	pt2.y = p3.y - p2.y;
	
	int temp = pt2.x*pt1.y - pt1.x*pt2.y;
	return temp;
}

int distance(Point p1, Point p2){
	return (p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y);
}

int _max(int a, int b){
	return a>b?a:b;
}

/************************************************************************/
/*旋转卡壳法 求凸多边形的直径，基本思想就是对每一条边，遍历每一个顶点到该边的两个断点的距离；
距离该边最远的顶点到两个端点的距离也最远，又因为顶点到边的距离是单峰函数，所以当下一个顶点小于
当前顶点时就退出，由于按逆时针方向输入各条边，所以最大的顶点也随逆时针反向旋转，不需要每次都遍历一遍
所有顶点，使用叉积来代表顶点到边的距离                                                                     */
/************************************************************************/
int rotating_calipers(vector<Point> &ch)
{
	ch.push_back(ch[0]);
	int size = ch.size();
	int p=1, result=0;
	for (int i=0; i<size-1; i++)
	{
		while(cross(ch[p+1],ch[i], ch[i+1])>cross(ch[p],ch[i], ch[i+1]))
			p = (p+1)%(size-1);
		result = _max(result,_max(distance(ch[p],ch[i]), distance(ch[p+1],ch[i+1])));
	}
	return result;
}

int main(int argc, char* argv[])
{
	Point p1, p2, p3, p4, p5, p6;
	p1.x = 1; p1.y=1;
	p2.x = 2; p2.y=5;
	p3.x = 4; p3.y=1;
	p4.x = 4; p4.y = 3;
	p5.x = 2; p5.y = 2;
	p6.x = 3; p6.y = 2;
	vector<Point> src, result;
	src.push_back(p1);
	src.push_back(p2);
	src.push_back(p3);
	src.push_back(p4);
	src.push_back(p5);
	src.push_back(p6);
	CalcConvexHull(src,result);
	cout << rotating_calipers(result) << endl;
	return 0;
}

参考资料：

算法导论

http://www.cnblogs.com/DreamUp/archive/2010/09/16/1828131.html
http://www.cnblogs.com/Booble/archive/2011/04/03/2004865.html
http://www.cppblog.com/staryjy/archive/2010/09/25/101412.html