计算几何学习笔记

终于会一点计算几何啦!

  几何是计算机科学的一个大板块。因为计算机非常不适合各种开根号的浮点运算,所以我们一般不在编程时使用解析几何处理问题,而是用计算几何。而计算几何的数学基础就是线性代数。本文章很适合对线性代数了解并不深刻的同学学习计算几何。

向量知识

  我们知道在数学中,通常将向量的起点定为坐标原点,用终点表示一个向量。这是一个向量和一个点在坐标轴中都用$(x,y)$

  为了方便将点和向量的概念区分开来,在算法竞赛中,我们通常这样定义点和向量。

struct Point {
  double x, y;
  Point(double _x = 0, double _y = 0) : x(_x), y(_y) {}
};

typedef Point Vector;

向量的加减乘除

  需要注意的是两个点$B-A$表示的是一个$\overrightarrow{AB}$

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);
}

向量和点的大小关系

  这个标题可能有点误导人,很显然的,向量和点没有办法比较大小。只有向量和向量或者点和点之间才行。因为计算机对浮点数的处理很不准确。所以我们定义一个三态函数$dcmp(x)$,来判断$x$的正负。点和点之间的大小关系是按照字典序的。

const double EPS = 1e-10;
bool dcmp(double x) {
  if (fabs(x) <= EPS) {
    return 0;
  } else {
    return x < 0 ? -1 : 1;
  }
}
bool operator == (const Point &A, const Point &B) {
  return !dcmp(A.x - B.x) && !dcmp(A.y - B.y);
}
bool operator < (const Point &A, const Point &B) {
  if (A.x == B.x) {
    return A.y < B.y;
  } else {
    return A.x < B.x;
  }
}

向量的点积及其相关应用

  向量和向量之间有两种乘法,一种是我们即将介绍的点积,还有一种是叉积。从几何意义上来理解点积会很容易。比如我们现在有两个向量$\vec A$和$\vec B$,我们将两个向量的点积认为是将$\vec B$投影到$\vec A$上,然后再将它们的长度相乘。我们用$\theta$表示这两个向量之间的夹角后,可以得到一个式子:

    $dot(\vec a, \vec b) = \left|A\right| × \left|A\right| × cos(\theta)$

  再通过一些线性代数的知识,我们发现,其实两个向量的点积,就是一个$1×2$的矩阵乘上一个向量。所以我们就可以写出下面的代码。

double dot(Vector A, Vector B) {
  return A.x * B.x + A.y * B.y;
}

点积求角度和长度

  再根据之前的定义,我们就可以很轻易地求求出$\theta$的大小和向量的长度(注意这里的$\theta$是用弧度制表示的)

double length(Vector A) {
  return sqrt(dot(A, A));
}
double angle(Vector A, Vector B) {
  return acos(dot(A, B) / length(A) / length(B));
}

向量的叉积及其相关应用

  和点积类似,叉积也有它的几何意义,而且更容易理解。就是求两个向量组成的平行四边形的有向面积。我们用两个向量组成一个$2 × 2$的矩阵。根据行列式的定义我们知道,他们的有向面积就是他们的行列式。然后根据
$det(\begin{bmatrix} a & b \\ c & d \end{bmatrix}) = ad - bc$,我们就可以写出两个向量叉乘的式子:

    $cross(\vec A, \vec B) = det(\begin{bmatrix} A_x & A_y \\ B_x & B_y \end{bmatrix}) = A_x × B_y - A_y × B_x$

  于是我们就有了下面的代码

double cross(Vector A, Vector B) {
  return A.x * B.y - A.y * B.x;
}

叉积求距离

  通过上面介绍的叉积,我们可以求出一个点到直线的距离。直线我们通常用两个点表示,设这两个点分别是$A$和$B$。再设点$P(x, y)$。所以$P$到$\overrightarrow{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$上的投影点为$Q$,$P$到线段$AB$距离为$D$。那么距离有三种情况:

$D=\begin{cases}
P到B距离 & Q在\overrightarrow{AB}上\\
P到A距离 & Q在\overrightarrow{BA}上 \\
P到直线AB距离 & Q在直线AB上
\end{cases}$

然后根据分类讨论的情况,我们可以用如下代码实现:

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 (dcmp(dot(V1, V2) < 0)) {
    return length(V2);
  } else if (dcmp(dot(V1, V3) > 0)) {
    return length(V3);
  } else {
    return fabs(cross(V1, V2)) / length(V1); 
  }
}

  上面求距离的方法虽然都出现了$P$的投影点$Q$,但我们都通过巧妙的方法避免了求$Q$,然而$Q$在某些时候是有用的,所以我们现在来考虑如何求它。设直线$AB=A+t\vec v$($v$为$\overrightarrow{AB}$),则可以设$Q=A+t’*\vec v$。因为$PQ$垂直于$AB$,所以我们可以得出式子:

    $dot(v, P - (A + t’*\vec v)) = 0$

  展开可得:

    $dot(v, P - A) - t’(dot(v, v)) = 0$

  解得:

    $t=\frac{dot(v, P - A)}{dot(v, v)}$

Point getLineProjection(Point P, Point A, Point B) {
  Vector V = B - A;
  return A + V * (dot(V, P - A) / dot(V, V));
}

叉积求多边形面积

  求面积的方法很简单,就是对每两个下标相邻的点求叉积。我们考虑这为什么是对的。首先,在凸多边形,这显然是成立的。但如果这是是一个凹多边形,那我们好像会多求一些面积。我们可以发现其实这些面积是负的(回忆一下有向面积的定义),所以这些并没能够使多边形面积变大。也就是我们求出来的面积还是多边形的面积。综上,这种方法在任何多边形中都适用。代码如下:

double polygonArea(Point *P, int n) {
  double area = 0;
  for (int i = 2; i < n; ++i) area += cross(P[i] - P[1], P[i + 1] - P[1]);
  return area;
}

直线求交

求两条线段交点

Point getLineIntersection(Point P, Vector V, Point Q, Vector W) {
  Vector U = P - Q;
  double t = cross(W, U) / cross(V, W), t2 = cross(V, U) / cross(V, W); 
  return P + V * t; // Q + W * t2
}

判断两条线段关系

  设两条线段的端点分别为$A_1,A_2$和$B_1,B_2$。只要每条线段的两个端点都在另一条线段的两侧就说明这两天线段相交。代码如下:

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);
  return dcmp(C1 * C2) < 0 && dcmp(C3 * C4) < 0;
}

判断点和线段的关系

bool onSegment(Point P, Point A1, Point A2) {
  return !dcmp(cross(A1 - P, A2 - P)) && dcmp(dot(A1 - P, A2 - P)) < 0;
}

向量旋转

  就是将终点旋转。

Vector rotate(Vector A, double rad) {
  return Vector(A.x * cos(rad) - A.y * sin(rad), A.x * sin(rad) + A.y * cos(rad));  
}

例题

UVA11178

题意

  给定一个三角形的三个顶点坐标。分别作出这个三角形每个内角的三等分线,它们相交能得到一个小的等边三角形。请求出小等边三角形的三个顶点坐标。
ZaM0AJ.jpg

题解

  按题意模拟即可。因为要求的三个顶点是对称的,所以我们只要能求出一个顶点的坐标,剩下两个就可以用相同的方法求得。我们考虑如何求点$D$坐标。我们可以先用点积求出$AB$和$CD$的夹角$\theta$,然后将$BC$旋转$\frac{\theta}{3}$得到$BD$,再用相同的方法得到$CD$,然后求$D$点坐标就变成了求$BD$和$CD$的交点。

代码

// Author: 23forever:  
#include <bits/stdc++.h>
using namespace std;

inline int read() {
  int x = 0, w = 1;
  char c = ' ';

  while (!isdigit(c)) {
    c = getchar();
    if (c == '-') w = -1;
  }

  while (isdigit(c)) {
    x = (x << 1) + (x << 3) + (c ^ 48);
    c = getchar();
  }

  return x * w;
}

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);
}

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;
}

Vector rotate(Vector A, double rad) {
  return Vector(A.x * cos(rad) - A.y * sin(rad), A.x * sin(rad) + A.y * cos(rad));  
}

Point getLineIntersection(Point P, Vector V, Point Q, Vector W) {
  Vector U = P - Q;
  double t = cross(W, U) / cross(V, W), t2 = cross(V, U) / cross(V, W); 
  return P + V * t; // Q + W * t2
}

Point A, B, C, D, E, F;

void init() {
  A.x = read(), A.y = read();
  B.x = read(), B.y = read();
  C.x = read(), C.y = read();
}

Point calc(Point A, Point B, Point C) {
  Vector V1 = A - B, V2 = C - B;
  double a1 = angle(V1, V2);
  V2 = rotate(V2, a1 / 3);

  Vector V3 = A - C, V4 = B - C;
  double a2 = angle(V3, V4);
  V4 = rotate(V4, -a2 / 3);

  return getLineIntersection(B, V2, C, V4);
}

int main() {
#ifdef forever23
  freopen("test.in", "r", stdin);
  //freopen("test.out", "w", stdout);
#endif
  int t = read();

  while (t--) {
    init();

    D = calc(A, B, C);
    E = calc(B, C, A);
    F = calc(C, A, B);
    printf("%f %f %f %f %f %f\n", D.x, D.y, E.x, E.y, F.x, F.y);
  }

  return 0;
}

总结

  计算几何使计算机拥有了更强大的图形处理能力,也使出题人能更多的毒瘤题(逃。但计算几何码量大,容易写挂,所以从学会到拿分是一个要经历大量的练习。之后我可能会再学一下凸包和半平面交,至少在考场上能拿个暴力分吧。

本博客所有文章均采用 CC BY-NC-SA 4.0 许可协议,转载请注明出处
本文链接:https://23forever.com/2019/07/04/Computational-geometry/