计算几何
[TOC]
二维几何基础
1.输入时计算坐标一般是实数,在写程序时一般用精度较高的double类型,不用float类型;double类型读入时用lf%格式,输出是用f%格式;
2.浮点数运算会产生精度误差,所以需要设置一个eps(偏差值),一般取1e-8;
3.判断一个浮点数是否等于0,不能直接用 == 0来判断,需要用sgn()函数判断是都小于eps;判断;
判断两个浮点数是否相等时,也不能直接用 == 判断是否相等,而是要用dcmp()函数判断是否相等;
fabs(x)就是返回x的绝对值
const double pi = acos(-1, 0); //高精度圆周率
const double eps = 1e-8; //偏差值,有时用1e-10
int sgn(double x){
if(fabs(x) < eps) return 0;//判断x是否等于0
else return x<0? -1:1; //判断x是正数还是负数
}
点和向量
1.点
二维平面中的点坐标(x,y);
struct point{
double x,y;
point(){}
point(double x,double y):x(x),y(y){}//构造函数方便赋值
};
2.两点之间的距离
把两点看成直角三角形的斜边
1.用库函数hypot()计算
hypot函数头文件math.h或cmath
hypot(a,b)的返回值为double类型，相当于sqrt(aa+bb)
double distance(point A,point B){
return hypot(A.x - B.x,A.y - B.y);
}
2.用sqrt()函数计算
double distance(point A,Point B){
return sqrt((A.x - B.x) *(A.x - B.x) + (A.y - B.y) * (A.y - B.y));
}
3.向量运算
在struct point中,对向量运算==重载==运算符;返回一个结构体;
typedef point vector;//向量和点的表达方式一样,所以改名区分
point operator - (point B) { return point(x - B.x,y - B.y); }
point operator + (point B) { return point(x + B.x,y + B.y); }
point operator * (double k) { return point(x * k,y * k); }
point operator / (double k) { return point(x / k,y / k); }
bool operator == (vector a, vector b) { return sgn(a.x - b.x) == 0 && sgn(a.y - b.y) == 0; }
点积和叉积
1.点积
A . B = |A| |B| cos;如果已知向量A和向量B,则不需要角度;
double dot(vector a, vector b) { return a.x * b.x + a.y * b.y; }
点积的作用
1.判断A与B的夹角是钝角还是锐角;
2.求向量A的长度 / 向量A长度的平方
double len(vector A){ return sqrt(dot(A,A));}
double len(vector A){ return dot(A,A);}
3.求向量A和B的夹角大小
double angle(vector A,vector B){ return acos(dot(A,B) / len(A) / len(B)); }
2.叉积
A x B = |A| |B| sin;角度表示向量A旋转到向量B所经过的夹角;
可以用右手定则进行判断正负;向量具有先后顺序,不同的顺序有不同的正负答案
double cross(vector A,vector B){ return A.x * B.y - A.y * B.x ;}
1.判断向量A,B的方向关系;
A x B > 0,B在A的逆时针; A x B < 0,B在A的顺时针;
A x B = 0,共线,方向相同或相反;
2.计算两向量构成的平行四边形的有向面积
3个点A,B,C,以A为公共点,得到向量B - A,C - A,所构成的面积如下
double area2(point A,point B,point C){ return cross(B - A,C - A); }
以B或C为公共点构成的平行四边形面积相同,正负不同;
3.计算3个点构成的三角形的面积
等于三个点构成的平行四边形的一半
double area1(point A,point B,point C){ return cross(B - A,C - A) / 2; }
4.向量旋转
求向量A逆时针旋转角度为rad后的向量
vector rotate(vector A,double rad){ return vector(A.x * cos(rad) - A.y * sin(rad),A.x * sin(rad) + A.y * cos(rad)); }
特殊情况是旋转90度
逆时针旋转90度 : rotate(A, pi /2),返回 vector(-A.y, A.x);
顺时针旋转90度 : rotate(A, - pi /2),返回 vector(A.y, - A.x);
求单位法向量 ,即逆时针旋转 90度,然后取单位值
vector normal(vector A){ return vector(- A.y / len(A), A.x / len(A)); }
5.用叉积检查两个向量是否平行或重合
bool parallel(vector A,vector B){ return sgn(cross(A, B)) == 0; }
点和线
1.一条直线可以用线上面两个点来表示;
- 点向式表示直线 : p = p0 + vt;由一个点p0 和方向向量 v决定;p0是直线上任意一点(x,y); v是方向向量,给定两个点A,B; 那么v = B - A;
如果t无限制,p就是直线;
如果t在[0,1]内,p就是A,B之间的线段;
如果p大于0,p就是射线;
- 斜截式 y = kx + b;
- 普通式 ax + by + c = 0;
//点向式
struct Line{
point p1,p2; //默认p是起点
vector v; //由p,q确定的方向向量
Line(){}
Line(point a,point b){ //构造函数
p1=a,p2=b,v=p1-p2;
}
point end(double t){ //点P=p+v*t
return p+v*t;
}
point spos(){ return p1; }//线段起点
point tpos(){ return p2; }//线段终点
double length(){ return sqrt(distance(p1,p2)); }//线段长度
void print(){
printf("Line:(%lf,%lf)->(%lf,%lf)\n",p1.x,p1.y,p2.x,p2.y);
}
};
//斜截式
struct Line{
point p1, p2;
Line(){}
Line(point p1,point p2) : p1(p1),p2(p2){}
Line(point p,double angle){
p1 = p;
if(sgn(angle - pi / 2) == 0) {p2 = (p1 + point(0, 1)) ;}
else{p2 = (p1 + point(1, tan(angle))) ;}
}
};
//普通式
struct Line{
point p1, p2;
Line(){}
Line(point p1,point p2) : p1(p1),p2(p2){}
Line(double a,double b,double c){
if(sgn(a) == 0){
p1 = point(0, -c / b);
p2 = point(1, -c / b);
}
else if(sgn(b) == 0){
p1 = point(- c / a, 0);
p2 = point(- c / a, 1);
}
else{
p1 = point(0, -c / b);
p2 = point(1, (-c - a) / b);
}
}
};
线段可以用直线表示,加以区分
typedef Line segment;
2. 点和直线的位置关系
用直线上两点p1 和 p2 与点p构成两个向量,用叉积的正负判断方向,从而得到位置关系;
int point_line_relation(point p,Line v){
int c = sgn(cross(p - v.p1,v.p2 - v.p1));
if(c < 0) return 1;//左边
if(c > 0) return 2;//右边
else return 0;//线上
}
3. 判断点是否在线上
先用叉积判断是否共线, 然后用点积看p和v的两个端点产生的角是否为180°;
bool point_on_seg(point p ,Line v){
return sgn(cross(p - v.p1,v.p2 - v.p1)) == 0 && sgn(dot(p - v.p1,p - v.p2)) <= 0;
}
4.点到直线的距离
已知点p和直线v(p1,p2),求p到v的距离;首先用叉积先求p,p1,p2构成的平行四边形的面积,然后用面积除以平行四边形的底边长,也就是线段(p1,p2)的长度,就得到了平行四边形的高,即p到直线的距离
double dis_point_line(point p,Line v){
return fabs(cross(p - v.p1,v.p2 - v.p1)) / distance(v.p1,v.p2);
}
5. 点在直线上的投影
已知直线上两点p1,p2以及直线外的一点p,求投影点p0;
k = |p0 - p1| / |p2 - p1|k就是线段p0p1 和 p2p1长度的比值
(p - p1) . (p2 - p1) = |p- p1| * |p2 - p1| * cos = |p0 - p1| *|p2 - p1|
| p0 - p1 | = ((p - p1) . (p2 - p1)) / | p2 - p1 |
k = |p0 - p1| / |p2 - p1| = ((p - p1) . (p2 - p1)) / | p2 - p1 |*| p2 - p1 |;
p0 = p1 + k * (p2 - p1) = p1 + ((p - p1) . (p2 - p1)) / | p2 - p1 |*| p2 - p1 | * (p2 - p1);
point point_line_proj(point p,line v){
double k = dot(v.p2 - v.p1,p - v.p1) / len2(v.p2 - v.p1);
return v.p1 + (v.p2 - v.p1) * k;
}
6.点关于直线的对称点
先求一个点p对一条直线v的镜像点,先求点p在直线上的投影q,再求对称点p’
point point_line_symmetry(point p,Line v){
point q = point_line_proj(p,v);
return point(2 * q.x - p.x;2 * q.y - p.y);
}
7.点到线段的距离
点p到线段AB的距离等于在以下三段距离中去最小
1.从p出发做AB的垂线,如果交点在AB线段上,这个距离就是最小值
2.p到A的距离 3.p到B的距离
double dis_point_seg(point p,segment v){
if(sgn(dot(p-v.p1,v.p2-v.p1))<0||sgn(dot(p-v.p2,v.p1-v.p2))<0)
return min(distance(p,v.p1),distance(p,v.p2));
return dis_point_line(p,v);
}
8.两条直线的位置关系
int Line_relation(Line v1,Line v2){
if(sgn(cross(v1.p2 - v1.p1,v2.p2 - v2.p1)) == 0){
if(point_line_relation(v1.p1,v2) == 0) return 1;//重合
else return 0;//平行
}
return 2;//相交
}
9.两条直线的交点
point cross_point(point a,point b,point c,point d){//Line 1:ab,Line 2:cd
double s1 = cross(b - a,c - a);
double s2 = cross(b - a,d - a);
return point(c.x * s2 - d.x * s1,c.y * s2 - d.y * s1) / (s2 - s1);//要保证不共线和平行
}
10.判断两个线段是否相交
利用叉积正负来判断;如果一条线段的两个端点在一条线段的两侧,那么这两个端点与另一线段产生的两个叉积正负相反,也就是说两个叉积相乘为负.
bool cross_segment(point a,point b,point c,point d){//Line 1:ab,Line 2:cd
double c1 = cross(b - a,c - a),c2 = cross(b - a,d - a);
double d1 = cross(d - c,a - c),d2 = cross(d - c,b - c);
return sgn(c1) * sgn(c2) < 0 && sgn(d1) * sgn(d2) < 0;
}
11.求两条线段的交点
先判断是否相交,若相交,问题转化为直线交点问题
多边形
判断点在多边形内部
1.转角法思想的一种实现:
以点p为起点引一条水平线,检查与多边形每条边的相交情况
检查以下3个参数:c = cross(p - j,i - j); u = i.y - p.y; v = j.y - p.y
叉积c用来检查p点在线段ij的左侧还是右侧;u,v用来检查经过p的水平线是否击穿过线段ij;
用num计数
if(c > 0 && u < 0 && v >= 0) num ++;
if(c < 0 && u >= 0 && v < 0) num –;
多边形的形状是由各个顶点的排列顺序决定的
int point_in_polygon(point pt,point *p,int n){
for(int i = 0;i < n;i ++){
if(p[i] == pt) return 3; //在多边形的顶点上
}
for(int i = 0;i < n;i ++){
Line v = Line(p[i], p[(i+1)%n]);
if(point_on_seg(pt,v)) return 2; //在多边形的边上
}
int num = 0;
for(int i = 0;i < n;++ i){
int j = (i + 1) % n;
int c = sgn(cross(pt - p[j],p[i] - p[j]));
int u = sgn(p[i].y - pt.y);
int v = sgn(p[j].y - pt.y);
if(c > 0 && u < 0 && v >= 0) num ++;
if(c < 0 && u >= 0 && v < 0) num --;
}
return num != 0; //1:点在内部, 0: 点在外部
}
2.pnpoly算法
假设多边形的坐标存放在一个数组里，首先我们需要取得该数组在横坐标和纵坐标的最大值和最小值，根据这四个点算出一个四边型，首先判断目标坐标点是否在这个四边型之内，如果在这个四边型之外，那可以跳过后面较为复杂的计算，直接返回false。
if (p.x < minX || p.x > maxX || p.y < minY || p.y > maxY) {
return false;
// 这个测试都过不了。。。直接返回false；
}
核心算法:
参数nvert 代表多边形有几个点。浮点数testx, testy代表待测试点的横坐标和纵坐标，vertx,verty分别指向储存多边形横纵坐标数组的首地址。
int pnpoly (int nvert, float *vertx, float *verty, float testx, float testy) {
int i, j, c = 0;
for (i = 0, j = nvert-1; i < nvert; j = i++) {
if ( ( (verty[i]>testy) != (verty[j]>testy) ) &&
(testx < (vertx[j]-vertx[i]) * (testy-verty[i]) / (verty[j]-verty[i]) + vertx[i]) )
c = !c;
}
return c;
}
3.射线法
以该点为起点引一条射线，与多边形的边界相交奇数次，说明在多边形的内部。
double dists(point p,point a,point b)//点到线段的距离
{
if (a==b) return len(p-a);
vector v1=b-a,v2=p-a,v3=p-b;
if (dcmp(dot(v1,v2))<0) return len(v2);
else if (dcmp(dot(v1,v3))>0) return len(v3);
return fabs(cross(v1,v2))/len(v1);
}
int pointin(point p,point* a,int n)
{
int wn=0,k,d1,d2;
for (int i=1;i<=n;i++)
{
if (sgn(dists(p,a[i],a[(i+1-1)%n+1]))==0) return -1;//判断点是否在多边形的边界上
k = sgn(cross(a[(i+1-1)%n+1]-a[i],p-a[i]));
d1 = sgn(a[i].y-p.y);
d2 = sgn(a[(i+1-1)%n+1].y-p.y);
if (k>0&&d1<=0&&d2>0) wn++;
if (k<0&&d2<=0&&d1>0) wn--;
}
if (wn) return 1;
else return 0;
}
多边形的面积
以原点为中心划分三角形,然后求多边形的面积
double polygon_area(point *p,int n){
double area = 0;
for(int i = 0;i < n;i ++){
area += cross(p[i],p[(i+1) % n]);
}
return area / 2;
}
求多边形的重心
hdu 1115
point polygon_center(point *p,int n){
point ans(0,0);
if(polygon_area(p,n) == 0) return ans;
for(int i = 0;i < n;i ++){
ans = ans + (p[i] + p[(i + 1) % n]) * cross(p[i],p[(i + 1) % n]);
}
return ans / polygon_area(p, n) / 6;
}
