计算几何基础入门详解

计算几何基础

1、最开始要提一下的还是精度问题，这在计算几何中可以说非常关键！！！

const double eps=1e-8; //我一般都用eps表示，控制在1e-8左右，当然有些题会有它的要求
const double PI=acos(-1.0); //π
int sgn(double x){ //判断一个数是否大于0
	if(fabs(x)<eps) return 0; //等于0
	if(x<0) return -1;	//小于0
	else return 1; //大于0
}

2、点的代码表示
用一个结构体来表示点，我们一般也把点积、叉积的重载操作符定义在点的结构体里，对之后写代码非常方便，而且有效，自己也可以自定义去重载一些运算符
对于点旋转某个角度后得到的坐标，本人不给予证明，感兴趣的同学请自行证明，或者百度 证明过程
关于点旋转这个问题让我想到我做过的一道最小椭圆覆盖的题，给出本人博客链接，感兴趣的同学可以去看一下
博客链接

struct Point{ 
	double x,y; //横纵坐标
	Point(){}
	Point(double _x,double _y) {x=_x;y=_y;}
	Point operator -(const Point &b)const {return Point(x-b.x,y-b.y);} //减 
	double operator ^(const Point &b)const {return x*b.y-y*b.x;} //叉积 
	double operator *(const Point &b)const {return x*b.x+y*b.y;} //点积 
	void transXY(double B){ //求点旋转B弧度后的坐标(90度为π/2弧度)
		double tx=x,ty=y;
		x=tx*cos(B)-ty*sin(B); 
		y=tx*sin(B)+ty*cos(B);
	}
};

3、向量（既然讲到了点就顺便提一下）
向量是什么我就不解释了，大家中学时期都已经学过了，我只讲在程序中怎么去表示它
定义：

typedef Point Vector;

模长：就是向量的长度，记为|A|

double length(const Vector& v){
    return sqrt(v.x*v.x+v.y*v.y);
}

3、直线与线段

struct Line{
	Point s,e; //直线用两个点表示 
	Line(){}
	Line(Point _s,Point _e) {s=_s;e=_e;}
	pair<int,Point>operator &(const Line &b)const{ //判断两直线的关系 
		Point res=s;
		if(sgn((s-e)^(b.s-b.e))==0){
			if(sgn((s-b.e)^(b.s-b.e))==0) return make_pair(0,res); //两直线重合 
			else return make_pair(1,res); //两直线平行 
		}
		double t=((s-b.s)^(b.s-b.e))/((s-e)^(b.s-b.e));
		res.x+=(e.x-s.x)*t;
		res.y+=(e.y-s.y)*t;
		return make_pair(2,res); //两直线有交点，交点为res 
	}
};

顺便写了一个用法，大家可以自己动手试一试

int main(){
	Line lin[100];
	lin[1]=Line(Point(0,0),Point(1,1));
	lin[2]=Line(Point(0,9),Point(9,0)); 
	pair<int,Point>pr=lin[1]&lin[2];
	Point p=pr.second;
	if(pr.first==0) cout<<"重合"<<endl;
	else if(pr.first==1) cout<<"平行"<<endl;
	else cout<<"交点: ("<<p.x<<","<<p.y<<")"<<endl; 
	return 0;
}

4、点积与叉积的运用
点积很好理解，其实就是两个点的横坐标相乘+纵坐标相乘
我重点讲一下叉积：
叉积有几种比较常见的运用，比如
4.1、判断点在直线的哪一侧

通过这张图很清晰的可以看到顺时针叉积<0，逆时针叉积大于0
如果p在线段左侧 向量pa^pb>0；
如果p在线段右侧 向量pa^pb<0；

double operator ^(const Point &b)const{//叉积 
	return x*b.y-y*b.x;
} 
double operator *(const Point &b)const{//点积 
	return x*b.x+y*b.y;
} 

4.2、求多边形的面积

为了让大家更加熟悉叉积的操作我选了两道简单题
第一题：判断点在直线的哪一侧 本人博客
第二题：用叉积求多边形的面积 本人博客

讲到这里，计算几何最最最基础的部分大概是讲完了
现在我们来讲一下计算几何中的一些操作

5、两个点之间的距离
这个大家应该是都会的，不多说直接贴一下代码

//运用重载运算符 
double dist(Point a,Point b){
	return sqrt((a-b)*(a-b));
}
//未运用重载运算符 
double dist(Point a,Point b){
	return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}

6、判断两条线段是否相交
直接看代码会比较抽象，大家可以用笔画一画，就很容易理解代码的含义
有些题目会考一些严格意义上的相交，我们就要把前四行的=符号去掉
我觉得有必要把这题找出来让大家练习一下 本人博客

bool inter(Line l1,Line l2){
	return
	max(l1.s.x,l1.e.x)>=min(l2.s.x,l2.e.x) &&
	max(l2.s.x,l2.e.x)>=min(l1.s.x,l1.e.x) &&
	max(l1.s.y,l1.e.y)>=min(l2.s.y,l2.e.y) &&
	max(l2.s.y,l2.e.y)>=min(l1.s.y,l1.e.y) &&
	sgn((l2.s-l1.e)^(l1.s-l1.e))*sgn((l2.e-l1.e)^(l1.s-l1.e))<=0 &&
	sgn((l1.s-l2.e)^(l2.s-l2.e))*sgn((l1.e-l2.e)^(l2.s-l2.e))<=0;
}

备注：这是去年我学长的讲课笔记，欢迎参考
判断线段相交核心代码（前四行判断是否在坐标轴上有重合，后两行判断是否跨立）：
判断线段相交：（思想）
① 首先两条线段在X轴和Y轴上的投影必须必须有重合。
② 线段的两个端点必须在另一个线段的两侧（相互跨立）
所以假设固定一条线段，把另一条线段的两个点与线段做叉积操作，如果结果同号就说明两个点在同侧（（C-A）X（C-B））X （（D-A）X（D-B）） > 0 （同侧）
（此处的X 为叉乘的意思）

7、判断直线与线段是否相交
判断线段与直线是否相交：（思想）
① 线段的两个端点必须在直线的两侧，与判断线段是否相交相似（相互跨立）
为什么不用判断投影是否重合？？
答: 直线与线段在X轴和y轴上当且仅当直线与线段平行时才会仅仅在一条坐标轴上投影重合，但是当直线与线段平行时，线段与直线相互跨立是必定不存在的，由此得出结论：当判断线段与直线是否相交时只需判断是否跨立即可。

bool line_inter_seg(Line l1,Line l2) {
	return sgn((l2.s-l1.e)^(l1.s-l1.e))*sgn((l2.e-l1.e)^(l1.s-l1.e)) <= 0;
}

8、点到直线的距离
result表示点到直线上距离最近的点的坐标
证明过程如下：

通过上面截图，大家可以非常容易理解代码的含义

Point PointToLine(Point P,Line L){
	Point result;
	double t=((P-L.s)*(L.e-L.s))/((L.e-L.s)*(L.e-L.s));
	result.x=L.s.x+(L.e.x-L.s.x)*t;
	result.y=L.s.y+(L.e.y-L.s.y)*t;
	return result;
}

9、求点到线段最近的点，也就相当于求点到线段的最短距离
与上一点的证明过程类似
由点到直线的距离的计算我们可以得到关系比t，但与直线不同的是，这个点可能在坐标轴上与线段没有重合，即点到线段的距离就是点与线段的两个端点距离的最小值。此处对t作了一个特殊判断，当t∈[0,1]时，说明P点到线段所在直线的最短距离的点在线段上，就可以按照点到直线的距离函数的方法进行计算。当t不属于这个区间的时候，直接返回线段端点与P点较近的点即可。
请大家自己在证明一遍！！！

Point NearestPointToLineSeg(Point P,Line L){
	Point result;
	double t=((P-L.s)*(L.e-L.s))/((L.e-L.s)*(L.e-L.s));
	if(t>=0&&t<=1){
		result.x=L.s.x+(L.e.x-L.s.x)*t;
		result.y=L.s.y+(L.e.y-L.s.y)*t;
	}	
	else{	
		if(dist(P,L.s)<dist(P,L.e)) result=L.s;
		else result=L.e;
	}
	return result;
}

10、求多边形的面积
计算多边形面积代码（要求参数中点的集合必须是多边形顶点的顺时针或逆时针）

double CalcArea(Point p[],int n){
	double res=0;
	for(int i=0;i<n;i++) res+=(p[i]^p[(i+1)%n])/2;
	return fabs(res);
}

核心思想： 将多边形切割成多个小三角形，根据叉积的性质（叉积的值为两向量组成的平行四边形的面积），求出各个小三角形的面积相加。此处看似是直接对点进行遍历，其实是利用了坐标原点（0，0），p[i]^p[(i+1)%n 实际上是两个以原点为出发点的向量的叉积，再除二以后即为以向量为两条边的三角形的面积。

11、判断点是否在线段上

bool OnSeg(Point P,Line L){
	return
	sgn((L.s-P)^(L.e-P))==0 &&
	sgn((P.x-L.s.x)*(P.x-L.e.x))<=0 &&
	sgn((P.y-L.s.y)*(P.y-L.e.y))<=0;
}

核心思想： 当P点在线段上时（设线段上两个点为s和e点），P点与s和e点分别构成的向量应当是叉积为0的，同时P点的x与y应当介于s和e的x与y之间，满足所有条件即可证明点P在线段上。

12、判断点是否在凸多边形内
原理极其简单，就是我们传说中的暴力啦，自己看代码吧

int inConvexPoly(Point a,Point p[],int n){
	for(int i=0;i<n;i++){
		if(sgn((p[i]-a)^(p[(i+1)%n]-a))<0) return -1; //在凸多边形外 
		else if(OnSeg(a,Line(p[i],p[(i+1)%n]))) return 0; //在凸多边形边界上 
	}
	return 1; //在凸多边形内 
}

13、判断点是否在多边形内

int inPoly(Point p,Point poly[],int n){
	int cnt;
	Line ray,side;
	cnt=0; ray.s=p;
	ray.e.y=p.y; 
	ray.e.x=-100000000000.0; //构造一条一点p为起始点的射线 
	for(int i=0;i<n;i++){ //枚举边 
		side.s=poly[i];
		side.e=poly[(i+1)%n];
		if(OnSeg(p,side)) return 0; //在边界上 
		//如果平行轴则不用考虑 
		if(sgn(side.s.y-side.e.y)==0) continue;
		if(OnSeg(side.s,ray)){
			if(sgn(side.s.y-side.e.y)>0) cnt++;
		}
		else if(OnSeg(side.e,ray)){
			if(sgn(side.e.y-side.s.y)>0) cnt++;
		}
		else if(inter(ray,side)) cnt++;
	}
	if(cnt%2==1) return 1; //在多边形内 
	return -1; //在多边形外 
}

核心思想： 在欧式空间内，以一个点为顶点出发作一条射线，当射线穿越多边形时，可划分为两种状态，进入多边形和穿出多边形。当射线与多边形的接触点为偶数时（包括0），则说明点在多边形外，反之则在多边形内部。 此处有三个状态需要进行特殊处理：
① 射线的顶点在多边形上。
② 射线经过了多边形的一个顶点。
③ 多边形的边与射线平行。
函数中cnt变量保存的是射线与多边形的接触点的数量。同时为避免重复运算，当射线经过多边形的一个顶点时，仅当经过的点是循环中线段side的两个端点中纵坐标y较大的点时才算有效。（注： 函数中ray保存的是以p为顶点方向为负方向的射线，side保存的是多边形的一条边，是一个线段）

14、判断多边形是否为凸多边形

bool isconvex(Point poly[],int n){
	bool s[3];
	memset(s,false,sizeof(s));
	for(int i=0;i<n;i++){
		s[sgn((poly[(i+1)%n]-poly[i])^(poly[(i+2)%n]-poly[i]))+1]=true;
		if(s[0] && s[2]) return false;
	}
	return true;
}

核心思想： 当一个多边形为凸多边形时，以poly[i]为出发点到poly[i+1]和ploy[i+2] 构成的两个向量的叉积的符号应当是保持一致的。

15、将线段往左移动一定距离

void change(Point a,Point b,Point &c,Point &d,double p){//将线段ab往左移动距离p
    double len=dist(a,b);
    double dx=(a.y-b.y)*p/len;
    double dy=(b.x-a.x)*p/len;
    c.x=a.x+dx; c.y=a.y+dy;
    d.x=b.x+dx; d.y=b.y+dy;
}