计算几何基础及模板

1.先来点对于精度要求的辅助函数

double PI = acos(-1);
double INF = 1e20;
double EPS = 1e-6;

bool IsZero(double x) {
	return - EPS < x && x < EPS;
}	//	是否为0 

int dcmp(double x) {
    if(fabs(x) < EPS) return 0;
    else return x < 0 ? -1 : 1;
} 	// 判断是正是负 

 

2.矢量、点、直线的定义以及四则运算

#define Point Vector  // 点和向量的表示法都是一样的

struct Vector{
	double x,y;
	Vector(double x=0,double y=0) : x(x),y(y) {}
}; 

Vector operator + (Vector p,Vector q) {
	return Vector(p.x + q.x, p.y + q.y);
}	//	加 

Vector operator -(Vector p, Vector q) {
	return Vector(p.x - q.x, p.y - q.y);
}	//	减 

Vector operator *(double k, Vector p) {
	return Vector(k * p.x, k * p.y);
}    //  数乘

double operator *(Vector p, Vector q) {
	return p.x * q.x + p.y * q.y;
} 	//   点乘

struct Line{
	Point a,b;
	Line(Point a,Point b) : a(a),b(b) {}
}; //  直线的定义 由两个点构成

 

3.点积的功能

求同向还是异向；（>0同向，=0垂直，<0异向）

求投影；

求出投影后用勾股定理求点到直线距离；

double dot(Vector p, Vector q) {
	return p.x*q.x+p.y*q.y;
} 	//	点积

double project(Vector p, Vector n) {
	return dot(p, unit(n));  
} 	//	p在n上的投影长度 

 

4.叉积的功能

求面积(三角形面积等于任意两边求叉积/2)；

求顺时针方向还是逆时针方向；（即a->b是顺时针转还是逆时针转）

若a 逆时针旋转小于180 度可到b, 则结果为正，否则结果为负

判断是否在半平面上

注意使用使用叉积的时候外面一定要加上()，也就是要写成(p^q)

double operator ^(Vector p, Vector q) {
	return p.x * q.y - q.x * p.y;
}	//	叉乘 

double Cross(Vector p, Vector q) {
	return p.x * q.y - q.x * p.y;
}	//	叉积  

 

5.一些其他向量/点之间的运算

1.给两个向量，求围成的三角形面积

2.求两点之间的距离

3.求一个点旋转alpha角度得到的点

double area(Vector p, Vector q) {
	return (p^q)/2.0;
}

double dist(Point p, Point q) {
	return length(p - q);
}	//	两点之间的距离

Point rotate(Point b, Point a,double alpha)  {
	Vector p = b - a;
	return Point(a.x + (p.x * cos(alpha)- p.y * sin(alpha)),a.y + (p.x * sin(alpha) + p.y * cos(alpha)) );
}	//   b绕a旋转alpha得到c点，顺时针旋转，用旋转矩阵理解 

 

6.点与直线之间的一些操作

1.求点到直线的距离

2.判断点在向量AB的那一边

3.判断点是否在线段AB上

4.求点P在直线L上的垂足（用于求对称点）

5.求过点P关于L的垂线

//1.求点到直线的距离

//2.判断点在向量AB的那一边

//3.判断点是否在线段AB上

//4.求点P在直线L上的垂足

//5.求过点P关于L的垂线

double dist(Point p, Line l) {
	return fabs((p - l.a) ^ (l.b - l.a)) / length(l.b-l.a);
}   //  点到直线距离 


int sideOfLine(Point p, Point a, Point b)
{  
	double result = (b - a) ^ (p - a);
	if (IsZero(result)) return 0;  //p 在 AB上
	if (dcmp(result) > 0) return 1;  //p在AB左侧
	return -1;  //p在AB右侧
}	//判断p在直线 AB的哪一侧


bool OnSegment (Point p, Point a1, Point a2) {
    if (dcmp(Cross(a1-p,a2-p))) return 0;
    else 
	if (dcmp(p.x-min(a1.x,a2.x))>=0 && dcmp(p.x-max(a1.x,a2.x))<=0 &&
 dcmp(p.y-min(a1.y,a2.y))>=0 && dcmp(p.y-max(a1.y,a2.y))<=0 )  return 1;
    else return 0;
}     //   点是否在线段上


Point foot(Point p, Line l) {
	return l.a + project(p - l.a, l.b - l.a)*unit(l.b - l.a);
}	//  p点在l上的垂足


Line Vertical(Point p, Line l) {
	return Line(p,p + (rotate(l.b, l.a, PI / 2) - l.a));
}	//  求过p点的l的垂线 

 

7.直线与直线之间的操作

1.求两个直线的交点

2.求两个直线的夹角

Point intersect(Line l, Line m, string msg) {
	double x = area(m.a - l.a, l.b - l.a);
	double y = area(l.b - l.a, m.b - l.a);
	if (IsZero(x + y)) {
		if (IsZero(dist(l.a, m))) msg = "重合";
		else msg = "平行";
		return NULL;
	}
	return m.a + x / (x + y) * (m.b - m.a);
}  //  两条直线的交点，然后按照比例相加减即可 


double angle(Line l, Line m) {
	return acos( fabs(project(l.b - l.a, m.b - m.a)) / length(l.b - l.a));
}    //   两条线的夹角 

8.求多边形面积

使用有向面积计算

按逆时针方向顺次为各边指定方向

对于每条边AB，累加(A^B)/2

 

9.判断点是否在多边形内部

从这个点引一条射线，统计多边形穿过这条线正反多少次，设其在

逆时针穿过加一 顺时针穿过减一

最后wn!=0的话说明在内部

#include<cstdio>
#include<iostream>
#include<cstring>
#include<cmath>
#include<algorithm>
#define Point Vector 
using namespace std;

/*
为了简略篇幅
这边省略了一些基本函数
到时候看自己代码库就行
*/

Point poly[10000];
Point p;
int n,i; 

/*-----------------------------------------------------------------------------*/
int is_Point_in_Poly(Point p, Point *poly, int n) 
{
    int wn = 0;
    for(int i = 1; i <= n; i++)
    {
        if (OnSegment(p, poly[i], poly[(i+1)%n])) return 0;    //  注意这边不可以使用SideOfLine 
        int k = dcmp(Cross(poly[(i+1)%n]-poly[i], p-poly[i]));
        int d1 = dcmp(poly[i].y-p.y);
        int d2 = dcmp(poly[(i+1)%n].y-p.y);
        if (k>0 && d1<=0 && d2>0) wn++;
        if (k<0 && d2<=0 && d1>0) wn--;
    }
    if (wn!=0)  return 1;
    else return -1;
}

/*-----------------------------------------------------------------------------*/
int main()
{
	scanf("%d",&n);
	for (i=1; i<=n; i++) cin>>poly[i].x>>poly[i].y;
	cin>>p.x>>p.y;	
	int ans=is_Point_in_Poly(p,poly,n);
	if (ans>0) cout<<"内部"<<endl;
	if (ans==0) cout<<"边上"<<endl;
	if (ans<0) cout<<"外部"<<endl;	
} 

 

10.凸包

直观来看，把点集 S 中所有点钉在一个木板上，

用一个橡皮筋框起所有点来，一撒手，橡皮筋就表示出了 S 的凸包

任意两点的连线都在凸包内

凸包算是一个点集了

凸包的求法 基于水平序的Andrew算法（具体详解见刘汝佳p272）

1、先把所有点按x排序，x相同按y排序

2、然后依次加点，当新点在凸包左边就继续，否则就继续删除最新加入的点直到新点在左边

3、重复这一过程直到碰到最右边的点，这样求出了下凸包，然后反过来从p[n]做一次就求出了上凸包

#include<cstdio>
#include<iostream>
#include<cstring>
#include<cmath>
#include<algorithm>
#define Point Vector 
using namespace std;

/*
为了简略篇幅
这边省略了一些基本函数
到时候看自己代码库就行
*/

Point p[10000];
Point ch[10000];
int n,m,i;

bool cmp(Point p,Point q) {
	if (p.x<q.x) return 1;
	if (p.x==q.x && p.y<q.y) return 1;
	return 0;
}


/*-------------------------------------------------------------------------------*/
int ConvexHull(Point *p, int n, Point*ch )	//Andrew算法
{
    sort(p+1,p+n+1,cmp);
    int m = 0;
    for (int i = 1; i <= n; i++)
    {
        while (m>1 &&  sideOfLine(p[i],ch[m-2],ch[m-1])<=0 ) m--;  
 // 或者写成 dcmp( Cross(ch[m-1]-ch[m-2], p[i]-ch[m-2]) )<=0
        ch[m++]=p[i];
    }
    int k = m;
    for (int i = n-1; i >= 1; i--)
    {
        while (m>k &&  sideOfLine(p[i],ch[m-2],ch[m-1])<=0 ) m--;
        ch[m++]=p[i];
    }
    if (n>1) m--;
    return m;
}

int main()
{
	scanf("%d",&n);
	for (i=1; i<=n; i++) cin>>p[i].x>>p[i].y;
	m=ConvexHull(p,n,ch);
	cout<<"凸包"<<endl; 
	for (i=0; i<m; i++) cout<<ch[i].x<<' '<<ch[i].y<<endl;	
}