算法学习笔记之计算几何--三角形，多边形与圆

Introduction

相比前两节这一节的内容较少设计到算法知识，不过多边形与圆是集合中十分基础的研究对象，ACM竞赛中也经常涉及，所以在这里将常见的相关知识总结一下。

三角形

三角形的面积

直接使用公式：

S△ABC=|12AB→×AC→|

// 三角形面积double TriangleArea(Point a, Point b, Point c) {    return Det(b - a, c - a) / 2.0;}

判断点在三角形内

判断点P是否在三角形△ABC内主要有两种方法：

1. 用求凸包时运用过的叉积法
2. 面积法

第一种不在赘述，只要顺次求一遍叉积看结果是否为同一符号即可。第二种的话若点P在三角形△ABC内则有S△ABC=S△PAB+S△PBC+S△PCA。

三角形的心

三角形的四心总结如下：

名称 定义 图示 备注 内心 三内角角平分线交点 三角形内切圆圆心 外心 三边垂直平分线交点 三角形外接圆圆心 垂心 三条高线交点 重心 三条中线交点 中线被重心以1:2的比例划分（靠近角的一端较长）

三角形的四心有一首打油诗：

内心全靠角平分，
外心中点垂线伸，
垂心垂直画三高，
形心角连线中心。

求重心

求重心坐标十分简单，即为三角形三顶点坐标的平均值。

设三角形三顶点为A,B,C，重心坐标公式为：

G=13(∑xi,∑yi)

// 三角形重心Point Centroid(Point a, Point b, Point) {    return (a + b + c) / 3.0;}

求内心

求内心稍为麻烦。先来引入一个引理：

引理：

E为AB上一点，F为AC上一点，若AE:EB=m:l，AF:FC=n:l，则BF与CE的交点I的坐标为(lxA+mxB+nxCl+m+n,lxA+mxB+nxCl+m+n)

可以这么理解，为了使重心落在BF与CE的交点上，我们需要在A点上放上质量为l的重物，在B点上放上质量为m的重物，在C点上放上质量为n的重物，这样根据重心的公式就能得到交点的坐标。

AE:EB=b:a，AF:FC=c:a，根据面的引理可以得出内心计算公式：

I=(axA+bxB+cxCa+b+c,ayA+byB+cyCa+b+c)

// 三角形内心Point Incentre(Point a, Point b, Point c) {    double a_edge = (b - c).Norm();    double b_edge = (a - c).Norm();    double c_edge = (a - b).Norm();    return (a * a_edge + b * b_edge + c * c_edge) / (a_edge + b_edge + c_edge);}

求外心

过三角形各顶点的圆叫做三角形的外接圆，外接圆的圆心即三角形外心，外心到三顶点距离相等。

外心的计算公式为:

// 三角形外心Point Circumcentre(Point a, Point b, Point c) {    double a1 = b.x - a.x, b1 = b.y - a.y, c1 = (a1 * a1 + b1 * b1) / 2.0;    double a2 = b.x - a.x, b2 = b.y - a.y, c2 = (a1 * a1 + b1 * b1) / 2.0;    double d = a1 * b2 - a2 * b1;    return a + Point(c1 * b2 - c2 * b1, a1 * c2 - a2 * c1) / d;}

求垂心

可以根据垂心与外心，重心的关系得出（三角形的垂心、重心和外心共线，且重心在垂心和外心连线的三等分点处，这里不予证明）。

// 三角形垂心Point Orthocentre(Point a, Point b, Point c) {    return Centroid(a, b, c) * 3.0 - Circumcentre(a, b, c) * 2.0;}

多边形

点在多边形内

要判断点是否在多边形内，可以该点向右引一条射线，然后看与多边形的交点数。考虑到复杂的非凸多边形情况，还有一些特殊情况需要考虑，我们需要的是射线真真正正的穿过了多边形的边且没有被重复计数。

最一般的情况是这样：

特殊情况的第一种是点在多边形上，这种情况下无论是否把点在多边形上算作穿过多边形的边都会陷入两难的境地，因此需要对这种情况进行特殊处理。

无论是点在多边形的边上还是在多边形的点上都属于这种情况，对于这种情况我们需要特判。

第二种是射线经过多边形的顶点：

第三种情况是射线经过多边形的边：

对于后两种情况解决方法是我们在判断一条线段是否被射线穿过时，我们只要规定当且仅当一个端点完全位于射线的上方/下方，另一个端点位于射线的下方/上方或位于射线上。这样仅仅是经过了顶点而没有实际穿过便会被计数为0或2，不影响结果。而经过了顶点且实际穿过了边的便会计数为1.与边重叠的会被计数为0或2.而对于点位于边上的情况特殊处理即可。

// 判断点是否在多边形内，包括在边上bool IsInPolygon(Point p, vector<Point> &polygon) {  bool in_polygon = false;  for (int i = 0; i < polygon.size() - 1; i++) {    // 判断点是否和顶点重合    if (DoubleCmp(p.x - polygon[i].x) == 0 && DoubleCmp(p.y - polygon[i].y) == 0) {      return true;    }    if (DoubleCmp(p.x - polygon[i + 1].x) == 0 && DoubleCmp(p.y - polygon[i + 1].y) == 0) {      return true;    }    // 判断线段的端点是否位于点的两侧    if ((DoubleCmp(polygon[i].y - p.y) >= 0 && DoubleCmp(polygon[i + 1].y - p.y) < 0) || (DoubleCmp(polygon[i].y - p.y) < 0 && DoubleCmp(polygon[i + 1].y - p.y) >= 0)) {      double x = polygon[i].x + (polygon[i + 1].x - polygon[i].x) / (polygon[i + 1].y - polygon[i].y) * (p.y - polygon[i].y); // 线段上与射线y坐标相同的点的x坐标      // 如果点在边上      if (DoubleCmp(x - p.x) == 0) {        return true;      }      // 如果射线穿过边界      if (DoubleCmp(p.x - x) < 0) {        in_polygon = !in_polygon;      }    }  }  return in_polygon;}

多边形有向面积

对于凸多边形，从一个顶点出发，将多边形分割为n−2个三角形即可。

double ConvexPolygonArea(vector<Point> &points) {  double area = 0;  for (int i = 1; i < points.size() - 1; i++) {    area += Det(points[i] - points[0], points[i + 1] - points[0]) / 2;  }  return area;}

其实对于任意多边形上述的方法也是正确的。为什么呢？因为我们计算的是有向面积，多算的部分最终会正负相消。如下图，假设A为基点，我们在计算的过程中计算的是12(AB→×AC→+AC→×AD→+AD→×AE→)，由于C在D的顺时针方向，所以面积为负，正好减掉了多出来的部分。这个求任意多边形的面积的公式称作surveyor公式，这篇文章有详细证明。

double PolygonArea(vector<Ponit> &points) {  double area = 0;  for (int i = 1; i < points.size() - 1; i++) {    area += Det(points[i] - points[0], points[i + 1] - points[0]) / 2;  }  return area;}

圆

通过圆心和半径可以确定圆。

// 圆结构体struct Circle {    Point center;    double r;    Circle(Point center = Point(0, 0), double r = 0) : center(center), r(r) {}    // 根据圆心角求坐标点    Point GetPoint(double a) {      return Point(center.x + cos(a) * r, center.y + sin(a) * r);    }};

直线与圆的交点

设直线上两点A和B，圆的圆心为C，半径为r。求直线与圆的交点的解析方法是解方程组，设交点P为A+t(B−A)，代入圆的方程整理后得(at+b)2+(ct+d)2=r2，进一步化简得et2+ft+g=0，其中e=a2+c2，f=2∗(a∗b+c∗d)，g=(b2+d2−rr)。然后用熟悉的解一元二次方程的方法即可解得。

代码如下：

// 直线与圆交点vector<Point> LineCircleIntersect(Line line, Circle circle) {  vector<Point> res;  double a = line.v.x, b = line.p.x - circle.center.x, c = line.v.y, d = line.p.y - circle.center.y;  double e = a * a + c * c, f = 2 * (a * b + c * d), g = b * b + d * d - circle.r * circle.r;  double delta = f * f - 4 * e * g;  if (DoubleCmp(delta) == 0) {    res.push_back(line.GetPoint(-f / 2 * e));  }  else if (DoubleCmp(delta) > 0) (    res.push_back(line.GetPoint(-f - sqrt(delta)));    res.push_back(line.GetPoint(-f + sqrt(delta)));  )  return res;}

圆与圆的交点

设两圆圆心分别为C1和C2，半径为r1和r2，圆心距为d，回顾初等几何知识可知两圆的位置关系有以下几种情况：

1. d = 0,内含且为同心圆
2. 0 < d < fabs(r1 - r2),内含
3. d = fabs(r1 - r2),内切
4. d = fabs(r1 + r2),外切
5. fabs(r1 - r2) < d < r1 + r2,相交
6. d > r1 + r2,外离

从下图可以很容易的看出可以根据C1P1和C1P2的圆心角算出交点的坐标，根据余弦定理，我们可以算出∠C2C1P1的值，根据向量C1C2→的极角a，加减∠C2C1P1就可以得到C1P1和C1P2的极角。

计算向量的极角（建议做成向量结构体的方法）：

double Angle(Vector v) {  return atan2(v.y, v.x);}

两圆交点的代码如下：

// 两圆交点vector<Point> CircleIntersect(Circle c1, Circle c2) {  vector<Point> res;  if (DoubleCmp(fabs(c1.r - c2.r) - d) <= 0 && DoubleCmp(c1.r + c2.r - d) >= 0) {    double d = (c1.center - c2.center).Norm();    double angle_c1c2 = (c2.center - c1.center).Angle(); // 向量c1c2的极角    double angle_c2c1p1 = acos((c1.r * c1.r + d * d - c2.r * c2.r) / (2 * c1.r * d));    Point p1 = c1.GetPoint(angle_c1c2 - angle_c2c1p1);    Point p2 = c1.GetPoint(angle_c1c2 + angle_c2c1p1);    res.push_back(p1);    if (p1 != p2) {      res.push_back(p2);    }  }  return res;}

过定点作圆的切线

代码如下：

vector<Vector> LineCircleTangent(Point p, Circle c) {  vector<Vector> res;  Vector u = c.center - p;  double dis = u.Norm();  // p在圆上，只有一条切线  if (DoubleCmp(dis - c.r) == 0) {    res.push_back(u.Rotate(pi / 2));  } else {    double angle = asin(c.r / dis);    res.push_back(u.Rotate(angle));    res.push_back(u.Rotate(-angle));  }  return res;}

两圆的公切线

代码如下：

vector<Point> CircleTangent(Circle a, Circle b) {}