计算几何
二维几何基础
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(a x a+b x b)
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中,对向量运算重载运算符;返回一个结构体;
struct Point {
double x, y;
Point() {}
Point(double x, double y) : x(x), y(y) {}
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); }//放大k倍
Point operator/(double k) { return Point(x / k, y / k); }//缩小k倍
bool operator==(Point B) { return sgn(x - B.x) == 0 && sgn(y - B.y) == 0; }
};
点积和叉积
1.点积
A . B = |A| |B| cos;如果已知向量A和向量B,则不需要角度;
typedef Point Vector;
double Dot(Vector A, Vector B) { return A.x * B.x + A.y * B.y;}
点积的作用
1.判断A与B的夹角是钝角还是锐角;
2.求向量A的长度 / 向量A长度的平方
double vector_length(Vector A) { return sqrt(Dot(A, A));}
double vector_length_square(Vector A) { return Dot(A, A);}
3.求向量A和B的夹角大小
double Angle(Vector A, Vector B) {
return acos(Dot(A, B) / vector_length(A) / vector_length(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 Area(Point A, Point B, Point C) {
return Cross(B - A, C - A);
}
以B或C为公共点构成的平行四边形面积相同,正负不同;
3.计算3个点构成的三角形的面积
等于三个点构成的平行四边形的一半
double Area(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));
}
求向量A顺时针旋转角度为rad后的向量
Point rotate(Point a, double angle){
return Point(a.x * cos(angle) + a.y * sin(angle), -a.x * sin(angle) * cos(angle));
}
特殊情况是旋转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 / vector_length(A), A.x / vector_length(A));
}
5.用叉积检查两个向量是否平行或重合
bool Parallel(Vector A, Vector B) {
return sgn(Cross(A, B)) == 0;
}
6.ToLeftTest
//判断折线bc是不是向ab的逆时针方向（左边）转向凸包构造时将会频繁用到此公式
bool ToLeftTest(Point a, Point b, Point c){
return Cross(b - a, c - 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;
//点向式(根据参数t来控制)
struct Line {
Point v, p;
Line(Point v, Point p) : v(v), p(p) {}
Point point(double t){
return v + (p - v) * t;
}
};
struct Line {
Point p1, p2;
Line() {}
//根据端点确定直线
Line(Point p1, Point p2) : p1(p1), p2(p2) {}
//根据一个点和倾斜角angel确定直线，0 <= angel < pi
Line(Point p, double angel) {
p1 = p;
if (sgn(angel - pi / 2) == 0)p2 = (p1 + Point(0, 1));
else p2 = p1 + Point(1, tan(angel));
}
//ax + by + c = 0
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构成两个向量,用叉积的正负判断方向,从而得到位置关系;
//点和直线关系:1 在左侧;2 在右侧;0 在直线上。根据视线从p1向p2看的左右
int Point_line_relation(Point p, Line v) {
int c = sgn(Cross(p - v.p1, v.p2 - v.p1));
if (c < 0)return 1; //1：p在v的左边
if (c > 0)return 2; //2：p在v的右边
return 0; //0：p在v上
}
3. 判断点是否在线上
先用叉积判断是否共线, 然后用点积看p和v的两个端点产生的角是否为180°;
// 点和线段的关系：0 点p不在线段v上；1 点p在线段v上。
bool on_segment(Point p, Point a, Point b) {
return sgn(Cross(p - a, p - b)) == 0 && sgn(Dot(p - a, p - b)) <= 0;
}
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) / vector_length_square(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;//1 重合
else return 0;//0 平行
}
return 2; //2 相交
}
9.两条直线的交点
//求两直线ab和cd的交点
//调用前要保证两直线不平行或重合
//叉积为零则平行或重合
Point Cross_point(Point a, Point b, Point c, Point d) { //Line1:ab, Line2: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.判断两个线段是否相交
利用叉积正负来判断;如果一条线段的两个端点在一条线段的两侧,那么这两个端点与另一线段产生的两个叉积正负相反,也就是说两个叉积相乘为负.
//两线段是否相交：1 相交，0不相交
bool Cross_segment(Point a, Point b, Point c, Point d) {//Line1:ab, Line2: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 --;
多边形的形状是由各个顶点的排列顺序决定的
//判断点和任意多边形的关系: 3 点上; 2 边上; 1 内部; 0 外部
int Point_in_polygon(Point pt, Point *p, int n) { //点pt，多边形Point *p
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.射线法
以该点为起点引一条射线，与多边形的边界相交奇数次，说明在多边形的内部。
int pointin(Point p,Point *a,int n) {
int wn = 0;
for (int i = 0;i < n;i ++) {
if (sgn(on_segment(p, a[i], a[(i+1)%n]))==0) return -1;//判断点是否在多边形的边界上
int k = sgn(Cross(a[(i+1-1)%n+1]-a[i], p-a[i]));
int d1 = sgn(a[i].y-p.y);
int 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) { //Point *p表示多边形。从原点开始划分三角形
double area = 0;
for (int i = 0; i < n; i++)
area += Cross(p[i], p[(i + 1) % n]);
return area / 2; //面积有正负，不能简单地取绝对值
}
求多边形的重心
hdu 1115
//多边形重心
//将多边形三角剖分，算出每个三角形重心（三角形重心是3点坐标平均值）
//对每个三角形有向面积求加权平均值
Point Polygon_center(Point *p, int n) { //求多边形重心。Point *p表示多边形。
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.;
}
浙公网安备 33010602011771号