三维向量夹角计算_夜深人静写算法（四） 计算几何入门

一、前言

“

最近听过感触最深的一句话：思想能够到达的地方比光速还要快！月球到了！太阳到了！你能知道的就是你能到达的，所以我们才要扩充我们的知识面，你要知道你才能到啊！
  太阳之外是什么？不知道了？到不了了，原因是什么？因为无知！因为我们把时间和精力花在了喝酒、吃饭、K歌。
  如果你有足够大的勇气去说再见，放弃一切娱乐，从现在开始每天学习，生活一定会奖励你一个新的开始。人生最难的是改变自己。凤凰涅槃浴火重生的人都很了不起。

”

二、计算几何基本概念

• 计算几何是计算机科学的一个分支，以往的解析几何，是用代数的方法，建立坐标系去解决问题，但是很多时候需要付出一些代价，比如精度误差，而计算几何更多的是从几何角度，用向量的方法来尽量减少精度误差，例如：将除法转化为乘法、避免三角函数等近似运算 等等。

1、浮点数精度

1)double 代替 float

• c++ 中 double 的精度高于 float，对精度要求较高的问题，务必采用 double；

2)浮点数判定

• 由于浮点数(小数)中是有无理数的，即无限不循环小数，也就是小数点后的位数是无限的，在计算机存储的时候不可能全部存下来，一定是近似的存储的，所以浮点数一定是存在精度误差的(实际上，就算是有理数，也是存在误差的，这和计算机存储机制有关，这里不再展开，有兴趣可以参见我博客的文章：C++ 浮点数精度判定)；
• 两个浮点数是否相等，可以采用两数相减的绝对值小于某个精度来实现：

const double eps = 1e-8;bool EQ(double a, double b) {    return fabs(a - b) }

• 并且可以用一个三值函数来确定某个数是零、大于零还是小于零：

int threeValue(double d) {    if (fabs(d)         return 0;    return d > 0 ? 1 : -1;}

3)负零判定

• 因为精度误差的存在，所以在输出的时候一定要注意，避免输出 -0.00：

    double v = -0.0000000001;    printf("%.2lf\n", v);

• 避免方法是先通过三值函数确定实际值是否为0，如果是0，则需要取完绝对值后再输出：

    double v = -0.0000000001;    if(threeValue(v) == 0) {        v = fabs(v);    }    printf("%.2lf\n", v);

4)避免三角函数、对数、开方、除法等

• c++ 三角函数运算方法采用的是 CORDIC算法，一种利用迭代的方式进行求解的算法，其中还用到了开方运算，所以实际的算力消耗还是很大的，在实际求解问题的过程中，能够避免不用就尽量不用。
• 除法运算会带来精度误差，所以能够转换成乘法的也尽量转换为乘法运算。

2、点和向量

1)定义

• 本章为计算几何入门级，所以只介绍二维的情况，对于二维的点，只需要两个浮点数表示即可，点的定义如下：

typedef double Type;class Point2D {private:    Type x, y;public:    Point2D(Type _x, Type _y) : x(_x), y(_y) {}    ...};

• 两点相减就变成了向量，向量支持平移(可以将起点平移到坐标系原点(0,0))，所以向量的程序运算和点几乎无差别，可以直接用 typedef 定义：

typedef Point2D Vector2D;

• 数学表示如下：

2)四则运算

• 点(向量)的四则运算 加减乘除 如下：

加

• 实现的是两个向量相加，代码如下：

Point2D Point2D::operator+(const Point2D& other) const {    return Point2D(x + other.x, y + other.y);}

• 如图二-2-1所示，红色向量为两向量相加后的向量； 图二-2-1

减

• 实现的是两个向量相减，代码如下：

Point2D Point2D::operator-(const Point2D& other) const {    return Point2D(x - other.x, y - other.y);}

• 如图二-2-2所示，红色向量为两向量相减后的向量； 图二-2-2

乘

• 向量和向量的乘法分为点乘和叉乘，这个后面再展开；
• 这里的乘法为 向量 和 数字 的乘法，代码实现如下：

Point2D Point2D::operator *(const double &k) const {    return Point2D(x * k, y * k);}

• 如图二-2-3所示，向量乘法就是将向量扩大()或缩短() 倍； 图二-2-3

除

• 向量和向量的之间不存在除法；
• 向量 和 数字 的除法，可以转换成向量乘法(乘上 的倒数)。

2、向量的模

• 两点之间的欧几里得距离，为两点组成向量的模，代码实现如下：

double Vector2D::len() {    return sqrt(x*x + y*y);}

• 调用的时候，可以直接调用两点之间的减法，再调用 成员函数，实现如下：

Point2D a(0, 2), b(3, 4);double dist = (b - a).len();

• 优化：在单纯进行向量长度比较的时候，往往只需要知道哪个更长(或者更短)，

3、标准化

• 有时候为了方便计算，会把向量转换成 单位向量(模为 1 的向量)，此所谓标准化。实现如下：

Point2D Vector2D::normalize() {    double l = len();    if (threeValue(l)) {        x /= l, y /= l;    }    return *this;}

• 由于存在零向量，所以标准化过程中，在进行除法的时候需要先用三值函数判断向量的模是否为0，避免引起除零错误；

4、点乘

• 点乘(又叫点积、数量积、内积)，表示的是两个向量每一维相乘的加和，结果是个数字，二维的情况如下：

• 代码实现如下：

Type Point2D::operator*(const Point2D& other) const {    return x*other.x + y*other.y;}

• (当然，点积同样适用于三维、甚至多维的情况)
• 点乘可以用来求两个向量的夹角，公式如下：

• 从以上公式可以看出，点乘的结果和夹角的关系：
• 1)点乘为0，夹角为90度；
• 2)点乘为1，夹角为0度；
• 3)点乘为-1，夹角为180度；
• 特殊的，任何向量和零向量的点乘都为 0；
• 下文会介绍 点乘 在线段判交中的应用；

5、叉乘

• 叉乘(又叫叉积、向量积，外积)，表示的是两个向量经过运算后，和这两个向量垂直的向量(和点乘不同，它的结果还是一个向量)，如下：

• 其中
• 当两个向量都在 的平面上时，，退化为二维的情况，则有
• 代码实现如下(因为 方向已经确定，为垂直原向量，所以返回值可以定义为一个数值类型)：

Type Point2D::X(const Point2D& other) const {    return x*other.y - y*other.x;}

• 在二维计算几何中，叉乘可以用来判断两个点是否在一个向量的同一侧，如图二-5-1所示：

图二-5-1

• 对于向量 ，假设有任意一个点 ，对原点到这个点做一个向量 。那么有叉乘：

• 1)，则 在 左侧(参考点A)；

• 2)，则 和 共线(参考点B)；

• 3)，则 在 右侧(参考点C)；

• 所以，只要两个点分别和对应向量做叉乘，再将结果相乘，如果值大于0，则代表同侧；小于零代表异侧；

• 叉乘还代表了两个向量组成的平行四边形的面积，有公式：

• 同样也代表了两个向量组成的三角形的面积的两倍。

6、旋转

• 通过 叉乘 和 点乘 联立方程组，可以求出一个点绕着原点旋转 角度后的点的位置；
• 令旋转前的点的位置为 ，旋转 角度后的点的位置 ，则有：

• 两个方程，两个未知数，求解得到：

• 表示成矩阵的形式为：

• 代码实现如下：

Point2D Point2D::turn(double theta) {    double costheta = cos(theta);    double sintheta = sin(theta);    return Point2D(        x*costheta - y*sintheta,         x*sintheta + y*costheta    );}

三、计算几何基础算法

1、线段判交

• 判断两个线段是否相交，在计算几何算法中有着广泛的应用，两个线段的相交情况有如下三种：

图三-1-1 * 这里用枚举来代表三种线段的相交关系：

enum SegCrossType {    SCT_NONE = 0,          // 不相交    SCT_CROSS = 1,         // 相交于内部    SCT_ENDPOINT_ON = 2,   // 其中一条线段的端点在另一条上};

• 并且定义线段的数据结构如下：

class Segment2D {    Point2D s, t;public:    ...};

1)相交于内部

• 首先对于两个线段，选取其中一个线段作为向量，另外一个线段的两端点作为测验点，如果这两个测验点在向量两边，那么可以肯定，这两个线段一定是不相交的(检测是否在向量两边可以用上文提到的叉乘)，如图三-1-2所示： 图三-1-2
• 然后，选取另外一个线段向量作为向量重复同样的测验；
• 如果两次测验都满足在向量两边，则这两个线段必然有交点；
• 这种测验被称为 跨立测验， c++ 代码实现如下：

Type Segment2D::cross(const Point2D& p) const {    return (p - s).X(t - s);}bool Segment2D::lineCross(const Segment2D& other) const {    return threeValue(cross(other.s)) * threeValue(cross(other.t)) == -1;}

• cross实现的是点 在向量 的哪边；
• lineCross判断两个点是否在向量的两边(通过取三值函数后，相乘为 -1 来实现)；

图三-1-3

• 最后，需要两次跨立测验都成立，才算满足线段相交的条件，如图三-1-3所示；

2)相交于端点

• 当一个线段的端点 P 在另一个线段 (S, T) 上时，只有三种情况：
• 1)P = S;
• 2)P = T；
• 3)∠SPT = 180°；
• 利用叉乘 和 点乘 进行组合判断即可，代码实现如下：

bool Segment2D::pointOn(const Point2D& p) const {    // 满足两个条件：    //  1.叉乘为0，    (p-s)×(t-s) == 0    //  2.点乘为-1或0，(p-s)*(p-t) <= 0    return cross(p) == 0 && (p - s)*(p - t) <= 0;}

• 然后对四个点分别进行端点相交判断；

3)线段判交实现

• 两次 跨立测验 和 四次 端点测验，实现如下：

SegCrossType Segment2D::segCross(const Segment2D& other) {    if (this->lineCross(other) && other.lineCross(*this)) {        // 两次跨立都成立，则必然相交与一点        return SCT_CROSS;    }    // 任意一条线段的某个端点是否在其中一条线段上，四种情况    if (pointOn(other.s) || pointOn(other.t) ||        other.pointOn(s) || other.pointOn(t)) {        return SCT_ENDPOINT_ON;    }    return SCT_NONE;}

2、多边形面积

• 多边形的面积可以通过将多边形进行三角剖分后，计算每个三角形的面积加和后得到，如图三-2-1所示： 图三-2-1
• 求三角形面积采用叉乘来实现；
• C++代码实现如下：

struct Polygon {    int n;    Point2D p[MAXP];    double area();};double Polygon::area() {    double sum = 0;    p[n] = p[0];    for (int i = 0; i         sum += p[i].X(p[i + 1]);    return sum / 2;}

• 这里的点是按照 进行极角排序的，并且保证 ；
• 如果点是按照顺时针排序，或者凹多边形的情况，求出来的面积可能是负数，如果进行输出的时候需要取下绝对值；

3、凸多边形判定

• 对于一个凸多边形而言，任意一条边作为向量，其它所有的顶点一定是在他的同侧的，所以可以利用这个性质来判断一个多边形是否是凸多边形； 图三-3-1

// 是否凸多边形bool Polygon::isConvex() {    bool s[3] = { false, false, false };    p[n] = p[0], p[n + 1] = p[1];    for (int i = 0; i         s[threeValue((p[i + 1] - p[i]) * (p[i + 2] - p[i])) + 1] = true;        // 叉乘有左有右，肯定是凹的        if (s[0] && s[2]) return false;    }    return true;}

4、点在多边形内判定

• 通过从无限远的地方引入一个随机点，然后和判定点做一条线段，把多边形的每条边和这条线段做判交检测，如果相交得到的交点数量，如果是奇数，说明在多边形内；如果是偶数，说明不在多边形内；
• 射线法同样也适用于凹多边形；

5、多边形点的逆时针转换

• 如果点是按照顺时针排布的，那么利用叉乘计算出来的面积为负数，根据这个特点可以将多边形的点进行反序，代码如下：

// 转成逆时针顺序void Polygon::convertToCounterClockwise() {    if (area() >= 0) {        return;    }    for (int i = 1; i <= n / 2; ++i) {        Point2D tmp = p[i];        p[i] = p[n - i];        p[n - i] = tmp;    }}

四、计算几何算法和应用

• 介于篇幅，以下算法只简单进行应用介绍，具体的算法实现将在后续章节进行详细展开。

1、凸包

• 凸包可想象为一条刚好包著所有点的橡皮圈，它一定是一个凸多边形。
• 凸包一般用来进行点集筛选，比如求解点集中的 最大面积三角形、最大周长三角形，因为这些点一定是在点集的凸包上，可以先进行凸包筛选降低点集数量，减少算法运行时间。

2、旋转卡壳

• 旋转卡壳一般配合凸包，用来求 凸多边形直径(点集的最远点对)、宽、凸多边形间最大距离、凸多边形间最小距离、最小面积外接矩形、最小周长外接矩形等等；

3、半平面交

• 半平面交一般可以用来求解以下问题：
• 1、两个多边形的交、并。
• 2、多边形的核：求解一个区域，使得这个区域内的任意点可以看到给定图形的任意角落。
• 3、求可以放进凸多边形的最大内接圆。
• 4、一些线性规划问题。

4、圆和多边形交

• 将多边形进行三角剖分以后，分别和圆进行求交，转换成三角形和圆相交的问题。

5、最小包围球

• 最小覆盖圆的三维情况，利用四点确定一个球体，然后判断其它点是否在这个球体中，用随机算法进行优化。

6、模拟退火

• 模拟退火是一种概率算法，在某个解空间内寻找最优解。比如求解三角形和圆的最大相交面积，寻找最近的最远点等等。