计算几何

二维几何基础

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.一条直线可以用线上面两个点来表示;

  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就是射线;

  1. 斜截式 y = kx + b;
  2. 普通式 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;

image-20220120205516703

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.;
}
posted on 2022-01-27 15:26  沐白白  阅读(179)  评论(0)    收藏  举报