计算几何（基础部分）

计算几何（基础部分）

• 计算几何
• • 注意
  • 点与向量模板
  • 点与线
  • • 直线
    • • 定义
    • 线段与射线
    • 点与线模板
  • 多边形
  • • 三角形
    • • 三角形面积
      • 三角形四心
  • 普通多边形（通常按照逆时针储存所有顶点）
  • • 判断点在凸多边形内
    • Pick定理
    • 普通多边形模板
  • 圆
  • • 圆模板

计算几何

注意

• 输出时一定要小心不要输出 − 0 -0−0，比如

  double a = -0.000001;
  printf("%.4f", a);
  

• 使用反三角函数时，要注意定义域的范围，比如，经过计算 x = 1.000001

  double x = 1.000001;
  double acx = acos(x);
  //可能会返回runtime error
  
  //此时我们可以加一句判断
  double x = 1.000001;
  if(fabs(x - 1.0) < eps || fabs(x + 1.0) < eps)
  	x = round(x);
  double acx = acos(x);
  
  

点与向量模板

struct Point{
    double x, y;
    Point(double x = 0, double y = 0):x(x),y(y){}
};
typedef Point Vector;
Vector operator + (Vector A, Vector B){
    return Vector(A.x+B.x, A.y+B.y);
}
Vector operator - (Point A, Point B){
    return Vector(A.x-B.x, A.y-B.y);
}
Vector operator * (Vector A, double p){
    return Vector(A.x*p, A.y*p);
}
Vector operator / (Vector A, double p){
    return Vector(A.x/p, A.y/p);
}
bool operator < (const Point& a, const Point& b){
    if(a.x == b.x)
        return a.y < b.y;
    return a.x < b.x;
}
const double eps = 1e-6;
int sgn(double x){
    if(fabs(x) < eps)
        return 0;
    if(x < 0)
        return -1;
    return 1;
}
int dcmp(double x, double y){
    if(fabs(x - y) < eps)
        return 0;
    if(x > y)
        return 1;
    return -1;
}
bool operator == (const Point& a, const Point& b){
    if(sgn(a.x-b.x) == 0 && sgn(a.y-b.y) == 0)
        return true;
    return false;
}
double Dot(Vector A, Vector B){//点乘
    return A.x*B.x + A.y*B.y;
}
double Length(Vector A){//模长
    return sqrt(Dot(A, A));
}
double Angle(Vector A, Vector B){//向量角
    return acos(Dot(A, B)/Length(A)/Length(B));
}
double Cross(Vector A, Vector B){//叉乘
    return A.x*B.y-A.y*B.x;
}
double Area2(Point A, Point B, Point C){//有向平行四边形面积
    return Cross(B-A, C-A);
}
Vector Rotate(Vector A, double rad){//rad为弧度 且为逆时针旋转的角
    return Vector(A.x*cos(rad)-A.y*sin(rad), A.x*sin(rad)+A.y*cos(rad));
}
Vector Normal(Vector A){//向量A左转90°的单位法向量
    double L = Length(A);
    return Vector(-A.y/L, A.x/L);
}
bool ToLeftTest(Point a, Point b, Point c){//判断bc是不是ab的左转向
    return Cross(b - a, c - b) > 0;
}

点与线

直线

定义

点向式:x₀ + y₀ + v_x t + v_y t = 0,即直线可以用直线上的一个点P₀和方向向量v表示:P = P₀ + v*t 其中t为参数

线段与射线

利用带参数限制的直线点向式方程表示

点与线模板

struct Line{//直线定义
    Point v, p;
    Line(Point v, Point p):v(v), p(p) {}
    Point point(double t){//返回点P = v + (p - v)*t
        return v + (p - v)*t;
    }
};
//计算两直线交点
//调用前需保证 Cross(v, w) != 0
Point GetLineIntersection(Point P, Vector v, Point Q, Vector w){
    Vector u = P-Q;
    double t = Cross(w, u)/Cross(v, w);
    return P+v*t;
}
//点P到直线AB距离公式
double DistanceToLine(Point P, Point A, Point B){
    Vector v1 = B-A, v2 = P-A;
    return fabs(Cross(v1, v2)/Length(v1));
}//不取绝对值，得到的是有向距离
//点P到线段AB距离公式
double DistanceToSegment(Point P, Point A, Point B){
    if(A == B)
        return Length(P-A);
    Vector v1 = B-A, v2 = P-A, v3 = P-B;
    if(sgn(Dot(v1, v2)) < 0)
        return Length(v2);
    if(sgn(Dot(v1, v3)) > 0)
        return Length(v3);
    return DistanceToLine(P, A, B);
}
//点P在直线AB上的投影点
Point GetLineProjection(Point P, Point A, Point B){
    Vector v = B-A;
    return A+v*(Dot(v, P-A)/Dot(v, v));
}
//判断p点是否在线段a1a2上
bool OnSegment(Point p, Point a1, Point a2){
    return sgn(Cross(a1-p, a2-p)) == 0 && sgn(Dot(a1-p, a2-p)) <= 0;
}
//判断两线段是否相交
bool SegmentProperIntersection(Point a1, Point a2, Point b1, Point b2){
    double c1 = Cross(a2-a1, b1-a1), c2 = Cross(a2-a1, b2-a1);
    double c3 = Cross(b2-b1, a1-b1), c4 = Cross(b2-b1, a2-b1);
    //if判断控制是否允许线段在端点处相交，根据需要添加
    if(!sgn(c1) || !sgn(c2) || !sgn(c3) || !sgn(c4)){
        bool f1 = OnSegment(b1, a1, a2);
        bool f2 = OnSegment(b2, a1, a2);
        bool f3 = OnSegment(a1, b1, b2);
        bool f4 = OnSegment(a2, b1, b2);
        bool f = (f1|f2|f3|f4);
        return f;
    }
    return (sgn(c1)*sgn(c2) < 0 && sgn(c3)*sgn(c4) < 0);
}

多边形

三角形

三角形面积

​ 利用两条边叉积除以二取绝对值

$$\begin{gathered} S=\sqrt{p(p-a)(p-b)(p-c)},p=\frac{(a+b+c)}{2} \\ S=\frac{absinC}{2} \end{gathered}$$

三角形四心

**外心：**三边中垂线交点，到三角形三个顶点距离相同

**内心：**角平分线的交点，到三角形三边的距离相同

**垂心：**三条高线的交点

**重心：**三条中线的交点，到三角形三顶点距离的平方和最小的点，三角形内到三边距离之积最大的点

普通多边形（通常按照逆时针储存所有顶点）

判断点在凸多边形内

只需要判断点是否在所有边的左边（按逆时针顺序排列的顶点集）ToLeftTest（P,A,B)

Pick定理

皮克定理是指一个计算点阵中顶点在格点上的多边形面积公式该公式可以表示为
$$2S=2a+b-2$$
其中a表示多边形内部的点数，b表示多边形边界上的点数，S 表示多边形的面积。

常用形式
$$S=a+\frac{b}{2}-1$$

常用计算

给你多边形的顶点，问多边形内部有多少点
$$a=S-\frac{b}{2}+1$$

普通多边形模板

//多边形有向面积
double PolygonArea(Point* p, int n){//p为端点集合，n为端点个数
    double s = 0;
    for(int i = 1; i < n-1; ++i)
        s += Cross(p[i]-p[0], p[i+1]-p[0]);
    return s;
}
//判断点是否在多边形内，若点在多边形内返回1，在多边形外部返回0，在多边形上返回-1
int isPointInPolygon(Point p, vector<Point> poly){
    int wn = 0;
    int n = poly.size();
    for(int i = 0; i < n; ++i){
        if(OnSegment(p, poly[i], poly[(i+1)%n])) return -1;
        int k = sgn(Cross(poly[(i+1)%n] - poly[i], p - poly[i]));
        int d1 = sgn(poly[i].y - p.y);
        int d2 = sgn(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;
    return 0;
}

圆

圆模板

struct Circle{
    Point c;
    double r;
    Circle(Point c, double r):c(c), r(r) {}
    Point point(double a){//通过圆心角求坐标
        return Point(c.x + cos(a)*r, c.y + sin(a)*r);
    }
};

//求圆与直线交点
int getLineCircleIntersection(Line L, Circle C, double& t1, double& t2, vector<Point>& sol){
    double a = L.v.x, b = L.p.x - C.c.x, c = L.v.y, d = L.p.y - C.c.y;
    double e = a*a + c*c, f = 2*(a*b + c*d), g = b*b + d*d - C.r*C.r;
    double delta = f*f - 4*e*g;//判别式
    if(sgn(delta) < 0)//相离
        return 0;
    if(sgn(delta) == 0){//相切
        t1 = -f /(2*e);
        t2 = -f /(2*e);
        sol.push_back(L.point(t1));//sol存放交点本身
        return 1;
    }
    //相交
    t1 = (-f - sqrt(delta))/(2*e);
    sol.push_back(L.point(t1));
    t2 = (-f + sqrt(delta))/(2*e);
    sol.push_back(L.point(t2));
    return 2;
}

//求两圆的公切线
// a[i] 和 b[i]存放第 i 条公切线与 圆A 和 园B 的交点
// 返回值为切线的条数 如果没有切线返回值为-1
int getTangents(Circle A, Circle B,vector<Point> &a, vector<Point> &b){
    int cnt = 0;
    // 以A为半径更大的那个圆进行计算
    if(A.r < B.r) return getTangents(B, A, b, a);
    double d2 = Dot(A.c-B.c, A.c-B.c);  // 圆心距平方
    double rdiff = A.r - B.r;		// 半径差
    double rsum = A.r + B.r;		//半径和
    if(d2 < rdiff * rdiff) return 0; 	// 情况1，内含,没有公切线
    Vector AB = B.p - A.p;				// 向量AB，其模对应圆心距
    double base = atan2(AB.y, AB.x);		// 求出向量AB对应的极角
    if(d2 == 0 && A.r == B.r) return -1;// 情况3，两个圆重合，无限多切线
    if(d2 == rdiff * rdiff){ 			// 情况2，内切，有一条公切线
        a.push_back(A.point(base));			
        b.push_back(B.point(base));
        cnt++;
        return 1;
    }
    // 求外公切线
    double ang = acos((A.r - B.r) / sqrt(d2)); //求阿尔法
    // 两条外公切线
    a.push_back(A.point(base+ang)); b.push_back( B.point(base+ang)) ; cnt++;
    a.push_back(A.point(base-ang)); b.push_back(B.point(base-ang)) ; cnt++;
    if(d2 == rsum * rsum){  // 情况5，外切，if里面求出内公切线
        a.push_back(A.point(base)); b.push_back(B.point(pi+base)); cnt++;
    }
    else if(d2 > rsum * rsum){	//情况6，相离，再求出内公切线
        db ang = acos((A.r + B.r) / sqrt(d2));
        a.push_back(A.point(base + ang)); b.push_back(B.point(pi+base+ang));cnt++;
        a.push_back(A.point(base - ang)); b.push_back(B.point(pi+base-ang));cnt++;
    }
    // 此时，d2 < rsum * rsum 代表情况 4 只有两条外公切线
    return cnt;
}

//求两圆交点
void IntersectionOf2Circles(Circle c1, Circle c2, Point &P1, Point &P2)
{
	float a1, b1, R1, a2, b2, R2;
	a1 = c1.c.x;
	b1 = c1.c.y;
	R1 = c1.r;

	a2 = c2.c.x;
	b2 = c2.c.y;
	R2 = c2.r;

	//
	double R1R1 = R1*R1;
	double a1a1 = a1*a1;
	double b1b1 = b1*b1;

	double a2a2 = a2*a2;
	double b2b2 = b2*b2;
	double R2R2 = R2*R2;

	double subs1 = a1a1 - 2 * a1*a2 + a2a2 + b1b1 - 2 * b1*b2 + b2b2;
	double subs2 = -R1R1 * a1 + R1R1 * a2 + R2R2 * a1 - R2R2 * a2 + a1a1*a1 - a1a1 * a2 - a1*a2a2 + a1*b1b1 - 2 * a1*b1*b2 + a1*b2b2 + a2a2*a2 + a2*b1b1 - 2 * a2*b1*b2 + a2*b2b2;
	double subs3 = -R1R1 * b1 + R1R1 * b2 + R2R2 * b1 - R2R2 * b2 + a1a1*b1 + a1a1 * b2 - 2 * a1*a2*b1 - 2 * a1*a2*b2 + a2a2 * b1 + a2a2 * b2 + b1b1*b1 - b1b1 * b2 - b1*b2b2 + b2b2*b2;
	double sigma = sqrt((R1R1 + 2 * R1*R2 + R2R2 - a1a1 + 2 * a1*a2 - a2a2 - b1b1 + 2 * b1*b2 - b2b2)*(-R1R1 + 2 * R1*R2 - R2R2 + subs1));
	if(sgn(subs1)!=0)//分母不为0
	{
		P1.x = (subs2 - sigma*b1 + sigma*b2) / (2 * subs1);
		P2.x = (subs2 + sigma*b1 - sigma*b2) / (2 * subs1);
	
		P1.y = (subs3 + sigma*a1 - sigma*a2) / (2 * subs1);
		P2.y = (subs3 - sigma*a1 + sigma*a2) / (2 * subs1);
	}

}

//求两圆相交的面积
//通过计算两个圆相交所构成的两个扇形面积和减去其构成的筝形的面积
double AreaOfOverlap(Point c1, double r1, Point c2, double r2){
    double d = Length(c1 - c2);
    if(r1 + r2 < d + eps)
        return 0.0;
    if(d < fabs(r1 - r2) + eps){
        double r = min(r1, r2);
        return pi*r*r;
    }
    double x = (d*d + r1*r1 - r2*r2)/(2.0*d);
    double p = (r1 + r2 + d)/2.0;
    double t1 = acos(x/r1);
    double t2 = acos((d - x)/r2);
    double s1 = r1*r1*t1;
    double s2 = r2*r2*t2;
    double s3 = 2*sqrt(p*(p - r1)*(p - r2)*(p - d));
    return s1 + s2 - s3;
}

__EOF__

本文作者：waby
本文链接：https://www.cnblogs.com/waby/p/15857129.html
关于博主：评论和私信会在第一时间回复。或者直接私信我。
版权声明：本博客所有文章除特别声明外，均采用 BY-NC-SA 许可协议。转载请注明出处！
声援博主：如果您觉得文章对您有帮助，可以点击文章右下角【推荐】一下。您的鼓励是博主的最大动力！