计算几何-基础

作者: 阿臻同学 | 来源:发表于2019-05-08 23:26 被阅读0次

一个基于坐标系的几何画图网页，对于计算几何，这个网页用来画图、打草稿啥的挺好用。

向量

点积（内积）

线性代数中的通用公式

对于 $n$ 维的向量 $\vec a$ 和 $\vec b$ ， $a_i$ 和 $b_i$ 表示向量 $\vec a$ 和 $\vec b$ 在第 $i$ 个维度上的分量。

$\begin{bmatrix}a_1 \\\\a_2\\\\\ldots\\\\a_{n-1} \\\\a_n\end{bmatrix} \cdot \begin{bmatrix}b_1 \\\\b_2\\\\\ldots\\\\b_{n-1} \\\\b_n\end{bmatrix} = a_1b_1+a_2b_b+\ldots+a_{n-1}b_{n-1}+a_nb_n$

亦可化简为： $\vec a\cdot\vec b=\sum_{i=1}^{n}a_ib_i$

对于二维向量 $\vec a(x_1,y_1)$ 和 $\vec b(x_2,y_2)$

设向量 $\vec a$ 和 $\vec b$ 的模分别为 $|\vec a|$ 和 $|\vec b|$ ，其夹角为 $\theta$ ，则有

$\begin{align} \vec a\cdot\vec b &= |\vec a||\vec b|cos \theta \\\\ \theta &= arccos(\frac{\vec a\cdot\vec b}{|\vec a||\vec b|}) \end{align}$

$\vec a在\vec b上的投影：\frac{\vec a\cdot\vec b}{|\vec b|}$
向量自身的点积为其长度的平方。
点积与方向的关系：

$\vec a\cdot\vec b$	$\theta$	$\vec a$ 和 $\vec b$ 的方向
<0	钝角	反向
=0	直角	垂直
>0	锐角	同向

叉积（外积）

线性代数中的通用公式

涉及到行列式展开计算，这里用三维举例。
$\begin{bmatrix}x_1 \\\\y_1\\\\z_1\end{bmatrix} \times \begin{bmatrix}x_2 \\\\y_2\\\\z_2\end{bmatrix} = \begin{bmatrix}x_1y_2-x_2y_1 \\\\x_1z_2-x_2z_1\\\\z_1y_2-z_2y_1\end{bmatrix}$

对于二维向量 $\vec a(x_1,y_1)$ 和 $\vec b(x_2,y_2)$

$\vec a\times\vec b= x_1y_2-x_2y_1$

叉积可以用来求出法向量方向。
叉积为向量围成的平行四边形的面积（带符号）。
若 $\vec a\times\vec b>0$ ，则表示 $\vec a$ 可通过逆时针方向旋转至 $\vec b$ 方向，小于0则为顺时针，等于0则表示两个向量共线。

单位向量

对于向量 $\vec{AB}$ ，其单位向量 $\ unit(\vec{AB})=\frac{\vec{AB}}{|\vec{AB}|}$

投影

对于向量 $\vec{AB}$ 、 $\vec{AC}$ ， $\vec{AC}$ 在 $\vec{AB}$ 上的投影为 $\vec{AC} \cdot (\vec{AB}的单位向量)$

例题

POJ 2007：Scrambled Polygon

代码模板

// 二维向量
struct Vector
{
    double x, y;

    Vector ()
    {}

    Vector (double x, double y) : x(x), y(y)
    {}

    // 点积
    double operator* (Vector v)
    {
        return x * v.x + y * v.y;
    }

    // 叉积
    double operator^ (Vector v)
    {
        return x * v.y - y * v.x;
    }

    // 乘法 
    Vector operator* (double d)
    {
        return Vector(x * d, y * d);
    }

    // 加法
    Vector operator+ (Vector v)
    {
        return Vector(x + v.x, y + v.y);
    }

    // 减法
    Vector operator- (Vector v)
    {
        return Vector(x - v.x, y - v.y);
    }

    // 模
    double length ()
    {
        return sqrt((*this) * (*this));
    }

    // 单位化
    Vector unit ()
    {
        return (*this) * (1.0 / length());
    }

};

点

点与向量的表示方法相同，都可以通过两个坐标来表示，直接重用向量即可，无需重复定义。

绕点旋转

点绕点旋转，其实可以直接理解成 $\vec{AB}$ 旋转到 $\vec{AC}$ ，其中点 $A$ 为旋转的中心。

可以先将 $A$ 平移至原点，旋转完成后，再平移回原来位置。

设:

原向量 $\vec{AB}$ 与x轴夹角为 $\beta$ ，表示为 $(x,y)$
旋转后的向量 $\vec{AC}$ 与x轴夹角为 $\alpha+ \beta$ ,表示为 $(x^,,y^,)$
$|\vec{AB}|$ 为 $r$

有如下推导：

$\begin{align} x &= r\cdot \cos \beta \\\\ y &= r\cdot \sin \beta \\\\ x^, &= r \cdot \cos{(\alpha+\beta)}\\ &= r \cdot (\cos\alpha\cdot \cos\beta-\sin\alpha\cdot \sin\beta) \\ &= x\cdot \cos\alpha-y\cdot \sin\alpha \\\\ y^, &= r \cdot \sin{(\alpha+\beta)}\\ &= r \cdot (\sin\alpha\cdot \cos\beta+\cos\alpha\cdot \sin\beta) \\ &= x\cdot \sin\alpha+y\cdot \cos\alpha \end{align}$

代码模板

#define Point Vector

    // 作为Vector成员函数
    Vector rotate (double alpha)
    {
        return Vector(x + (x * cos(alpha) - y * sin(alpha)),
                      y + (x * sin(alpha) + y * cos(alpha)));
    }

例题

POJ 2624：4th Point

直线

两点连成线，用两个点相减即可表示。

点-线距离

对于点A、B所在直线 $\overline {AB}$ 和点P，P点到直线 $\overline {AB}$ 的距离有：

设 $P$ 到 $\overline {AB}$ 的距离为 $\ dist$
$\begin{aligned} S_{\vec {PA}、\vec{PB}所构成的平行四边形}&= \vec {PA}\times\vec {PB}\\\ &=|\vec {AB}|\times dist\qquad \mbox{平行四边形面积公式}\\\ \Rightarrow dist &= \frac{\vec {PA}\times\vec {PB}}{|\vec {AB}|} \end{aligned}$

点-线位置

对于点A、B所在直线 $\overline {AB}$ 和点P，利用前面介绍的叉积与向量方向关系的性质，构造向量 $\vec {PA}$ ，通过 $\vec {PA} \times \vec{AB}$ 来判断P和 $\overline {AB}$ 的位置关系：
$P和\overline {AB}的位置关系= \begin{cases} P在\overline {AB}左侧， &\mbox{$\vec {PA} \times \vec{AB}>0$}\\ P在\overline {AB}上， &\mbox{$\vec {PA} \times \vec{AB}=0$}\\ P在\overline {AB}右侧， &\mbox{$\vec {PA} \times \vec{AB}<0$}\\ \end{cases}$

过点做垂线

对于点A、B所在直线 $\overline {AB}$ 和点P：

垂线：将直线所在向量 $\vec {AB}$ 旋转 $\frac{\pi}{2}$ ，并平移至经过点 $P$ 即可。

垂足：利用点积求出 $\vec{AP}$ 在 $\vec{AB}$ 上投影的长度 $\ l$ ， $A+l\cdot (\vec{AB}单位向量)$ 即为垂足的位置。

直线-直线的位置

对于点 $A$ 、 $B$ 所在直线 $\overline {AB}$ 和C、D所在直线 $\overline {CD}$ ：

如果二者不相交（平行或相等），则 $A$ 和 $B$ 在 $\overline {CD}$ 同一侧，其叉积值和符号都相同，相减为 $0$ 。反之则不为零。
$\begin{align} S_{ACD} &= \vec{AC}\times\vec{AD} \\ S_{BCD} &= \vec{BC}\times\vec{BD} \\ S_{ACD}-S_{BCD} &= \begin{cases} 0 &,\mbox{$\vec{AB}\parallel\vec{CD}$或$\vec{AB}=\vec{CD}$}\\ 非0 &,\mbox{$\vec{AB}\nparallel\vec{CD} $}\\ \end{cases} \end{align}$

如果二者不相交，则再通过 $A$ 到 $\overline{CD}$ 的距离来判断二者是否相等。距离为 $0$ 则相等，否则平行。
若二者相交，通过相似三角形的比例关系，可求得则其交点位置为：

$A+\frac{S_{ACD}}{S_{ACD}+A_{BCD}}\cdot\vec{AB}$

代码

// 二维直线
struct Line
{
    Point x, y;

    // 通过两点构造直线
    Line (Point x, Point y) : x(x), y(y)
    {}

    // 点和方向向量来构造线
    static Line makeLine (Point x, Vector v)
    {
        return Line(x, x + v);
    }

    bool operator== (const Line &rhs) const
    {
        return x == rhs.x &&
               y == rhs.y;
    }

    bool operator!= (const Line &rhs) const
    {
        return !(rhs == *this);
    }

    // 线长度
    double length ()
    {
        return (y - x).length();
    }

    // 点到该直线的距离
    double dist (Point p)
    {
        return fabs((x - p) ^ (y - p)) / length();
    }

    // #define EPS 1e-6
    // 判断点和直线的位置
    int side (Point p)
    {
        double result = (y - x) ^(p - x);
        if (fabs(result) < EPS)
            return 0;   // 在线上
        else if (result > 0)
            return 1;   // 左侧
        else
            return -1;  // 右侧
    }

    // 过点做垂线
    Line vertical (Point p)
    {
        return makeLine(p, (y - x).rotate(PI / 2));
    }

    // 垂足
    Point foot (Point p)
    {
        Vector self = y - x;
        return x + self.unit() * self.project(p - x);
    }

    Point intersect (Line l)
    {
        double s1 = ((x-l.x) ^(l.y - l.x)) / 2 ;
        double s2 = ((l.y-l.x) ^(y - l.x)) / 2;
        if (fabs(s1 + s2) < EPS)
        {
            if (l.dist(x))
                return l.x;   // 重合
            else
                return l.y;   // 平行
        } else
            return x + (y - x) * (s1 / (s1 + s2)); // 交点
    }

};

线段

点-线段位置

对于点 $P$ 和线段 $AB$ ，构造向量 $\vec{PA}、\vec{AB}$ ，求叉积有：
$\vec{PA}\times\vec{AB}= \begin{cases} 0 &,\mbox{$P$和$AB$在同一条直线上(共线)}\\ 非0 &,\mbox{$P$在线段$AB$外}\\ \end{cases}$
对于共线的情况，再根据两个端点的坐标判断， $P$ 是否在线段 $AB$ 上，即满足以下情况的，即为点 $P$ 在线段 $AB$ 上：
$\begin{cases} \min(A.x,B.x) &\leq P.x &\leq\max(A.x,B.x) \\ \min(A.y,B.y) &\leq P.y &\leq\max(A.y,B.y) \\ \end{cases}$

直线-线段位置

对于直线 $AB$ 和线段 $CD$ ，二者位置关系的判断：

构造向量 $\vec{CA}、\vec{DA}$ ，并求叉积 $\vec{CA}\times\vec{AB}、\vec{CB}\times\vec{AB}$
若叉积都不为，且符号相同，则可认为点都在向量同一侧，或二者重合。
- 如果一个端点到直线的距离为 $0$ ，则二者重合。
- 否则，二者不相交。
否则二者相交。交点位置求法同直线-直线交点位置。

线段-线段位置

规范相交：两条线段恰有一个不是端点的公共点。

对于直线 $AB$ 和线段 $CD$ ，利用叉积的符号判断是否为规范相交，即任一条线段的两个端点都分属另一条线段的两侧，满足以下条件：
$\begin{cases} (\vec{AB}\times\vec{AC})&与&(\vec{AB}\times\vec{AD})&异号 \\ (\vec{CD}\times\vec{CA})&与&(\vec{CD}\times\vec{CB})&异号 \\ \end{cases} \\\\ \Rightarrow \begin{cases} (\vec{AB}\times\vec{AC})\cdot(\vec{AB}\times\vec{AD})&<&0 \\ (\vec{CD}\times\vec{CA})\cdot(\vec{CD}\times\vec{CB})&<&0 \\ \end{cases} \\\\$

如不满足以上条件，再判断两条线段的端点是否在另一条线段上，若是则为非规范相交，否则二者不相交。

例题

POJ 3304：Segments

POJ 1039：Pipe

代码

// 二维线段
struct Seg
{
    Point x, y;

    Seg (const Point &x, const Point &y) : x(x), y(y)
    {}

    // 判断点是否在线段内
    bool pointIn (Point p)
    {
        double cross = (x - p) ^(y - x);
        if (fabs(cross) < EPS)
            return 0;
        if (fabs(p.x - min(x.x, y.x)) < EPS &&
            fabs(p.y - min(x.y, y.y)) < EPS &&
            fabs(max(x.x, y.x) - p.x) < EPS &&
            fabs(max(x.y, y.y) - p.y) < EPS)
            return 1;
        return 0;
    }

    Point intersect (Line l)
    {
        int s1 = l.side(x);
        int s2 = l.side(y);
        if (s1 * s2 > 0)
        {
//            if (l.dist(x) < EPS)
//                return l.y; // 重合
//            else
            return Point(EPS / 10, EPS / 10);// 不相交
        } else if (s1 == 0 && s2 == 0)
        {
            return l.x;// 重合
        } else if (s1 == 0)   // 一个端点在直线上
        {
            return x;
        } else if (s2 == 0)  // 另一个端点在直线上
        {
            return y;
        } else
        {
            return Line(x, y).intersect(l); // 交点
        }
    }

    // 判断线段相交
    bool isIntersected (Seg s)
    {
        Point p1 = x;
        Point p2 = y;
        Point p3 = s.x;
        Point p4 = s.y;
        double s1 = (p3 - p1) ^(p4 - p1);
        double s2 = (p4 - p2) ^(p3 - p2);
        if (((p2 - p1) ^ (p3 - p1) * ((p2 - p1) ^ (p4 - p1))) < -EPS &&
            (((p4 - p3) ^ (p1 - p3)) * ((p4 - p3) ^ (p2 - p3))) < -EPS)
            return 1;   // 规范相交
        if (s.pointIn(p1) || s.pointIn(p2) || pointIn(p3) || pointIn(p4))
            return 1;   // 非规范相交

        return 0;
    }

    bool operator== (const Seg &rhs) const
    {
        return x == rhs.x &&
               y == rhs.y;
    }

    bool operator!= (const Seg &rhs) const
    {
        return !(rhs == *this);
    }

};

多边形

凸多边形

百度百科：凸多边形是一个内部为凸集的简单多边形。凸多边形（Convex Polygon）指如果把一个多边形的所有边中，任意一条边向两方无限延长成为一直线时，其他各边都在此直线的同旁，那么这个多边形就叫做凸多边形，其内角应该全不是优角，任意两个顶点间的线段位于多边形的内部或边上。

就是没有凹角的多边形，特点是所有边所在的直线都不会穿过其他边。

多边形面积

对多边形进行三角剖分：

设原点为 $O$
按逆时针方向为个条边指定方向
对于每条边 $\vec{AB}$ ，累加 $\frac{\vec{OA}\times\vec{OB}}{2}$

点-多边形位置

设多边形 $P$ 的顶点序列为 $\{p_1,p_2, \ldots ,p_n\}$ 、待判定的点 $q$ 。

射线法

算法思路

过 $q$ 做水平射线 $l$ ：

若 $l$ 与 $P$ 的边界不相交，则 $\ l$ 在 $P$ 的外部。
若相交，则再分析交点数的奇偶性。

交点数点相对于多边形的位置

奇数内部

偶数外部
需要注意的是，当多边形不是凸多边形时，需要考虑特殊情况：
1. 射线恰好经过 $P$ 的某个顶点
2. 射线恰好与 $P$ 的某条边重合

交点数	点相对于多边形的位置
奇数	内部
偶数	外部

算法特点

运算速度快
精度高
特殊情况较多

转角法

算法思路

沿多边形走一圈，累计绕点旋转了多少角度。（需要保证边的方向一致）

角度	点相对于多边形的位置
$0$	外部
$\pi$	点在多边形的边上
$\pm2\pi$	内部

假设向量 $\vec{a}、\vec{b}$ 其角度 $\ \theta\$ 通过点积公式可以求出:
$\theta = \arccos{\frac{\vec{a}\cdot\vec{b}}{|\vec{a}|\cdot|\vec{b}|}}$
（仅适用于凸包）或使用累乘叉乘，判断点和边的位置关系：

当点位于所有边向量的同一侧时（叉乘值同号），其在多边形内部。
如果没有出现在同一侧，则为外部（叉乘值异号）。
在多边形上的情况是可以直接判断出来的（叉乘为0）。

算法特点

几乎没有特殊情况
精度低、速度较慢（需要用到反三角函数、开方等）

求凸包（凸壳）

给定一个平面点集，要求找到一个最小的凸多边形，满足点集中所有点都在凸多边形内部或边上。

稳定凸包

对于已有的凸包中，无法通过在平面中添加一个点获取到更大的凸包，则该已有的凸包称为稳定凸包。

特性：凸包的每条边上，都有大于等于 $3$ 个的顶点。

具体可见例题【POJ 1228 Grandpa's Estate】。

半平面交

半平面

一条直线将一个平面分为2个半平面，直线是有向的，设直线方向左边的平面为我们研究的半平面（包含直线）。

半平面交

被（多个）半平面包含的点的集合。

性质：半平面交的结果是一个凸区域。

求交

对于每一条多边形的每一条边 $\vec{AB}$ 、表示半平面的直线 $\vec{l}$ ：

对多边形按照逆时针排序。
如果边 $\vec{AB}$ 的起点 $A$ 位于半平面 $l$ 中，则将点 $A$ 加入结果集。
求边 $\vec{AB}$ 与 $l$ 的交点，并将结果加入结果集。（即使点 $A$ 不在半平面中，也需要这步操作）

增量法

构造一个足够大的矩形。
依次用半平面和该矩形求交。

时间复杂度： $O(n^2)$

分治法（二分）

将 $n$ 个半平面分成 $2$ 个大小近似相等的集合。
递归地构造凸多边形区域 $A$ 与 $A'$ 。
合并 $A$ 与 $A'$ 。

时间复杂度： $O(n \cdot \log{n})$

例题

POJ 1569:Myacm Triangles
POJ 1113-Wall
POJ 1228-Grandpa's Estate

代码

// 二维多边形
struct Polygon
{
    // 点集
    vector<Point> points;

    Polygon (const vector<Point> &points) : points(points)
    {}

    Polygon ()
    {}

    /**
     * 根据已排序的点集来求半个凸壳（上半个或下半个）
     * @param pts 已排序的点集
     * @return 半个凸包中的点，有序，逆时针序
     */
    deque<Point> getHalfConvexHull (const vector<Point> &pts) const
    {
        deque<Point> deq;
        int si = pts.size();
        int i = 2;

        deq.pb(pts[0]);
        deq.pb(pts[1]);
        while (i < si)
        {
            Point nt = pts[i];
            ++i;
            while (deq.size() >= 2)
            {
                Point p2 = deq[deq.size() - 2], p1 = deq.back();
                Vector b = p2 - p1, a = p1 - nt;
                double d = (b ^ a);
                if (d <= 0)
                    deq.pop_back();
                else
                    break;
            }

            deq.pb(nt);
        }
        return deq;
    }

    /**
     * 水平序扫描法
     * 求该多边形的凸包
     * @return 凸包的点集，有序，逆时针序
     */
    vector<Point> getConvexHull ()
    {
        vector<Point> pts(points);
        vector<Point> con;


        // 2个点无法构成凸包
        if (pts.size() <= 2)
        {
            return con;
        }

        // 下凸壳
        sort(pts.begin(), pts.end());
        deque<Point> temp = getHalfConvexHull(pts);
        for (int i = 0; i < temp.size(); ++i)
        {
            con.pb(temp[i]);
        }
        temp.clear();


        // 上凸壳
        sort(pts.begin(), pts.end(), greater<Point>());
        temp = getHalfConvexHull(pts);
        for (int i = 1; i < temp.size() - 1; ++i)
        {
            con.pb(temp[i]);
        }

        return con;
    }
};

计算几何-基础

向量

点积（内积）

线性代数中的通用公式

对于二维向量和

叉积（外积）

线性代数中的通用公式

对于二维向量和

单位向量

投影

例题

代码模板

点

绕点旋转

代码模板

例题

直线

点-线距离

点-线位置

过点做垂线

直线-直线的位置

代码

线段

点-线段位置

直线-线段位置

线段-线段位置

例题

代码

多边形

凸多边形

多边形面积

点-多边形位置

射线法

算法思路

算法特点

转角法

算法思路

算法特点

求凸包（凸壳）

稳定凸包

半平面交

半平面

半平面交

求交

增量法

分治法（二分）

例题

代码

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读

对于二维向量 $\vec a(x_1,y_1)$ 和 $\vec b(x_2,y_2)$

对于二维向量 $\vec a(x_1,y_1)$ 和 $\vec b(x_2,y_2)$