计算几何基础2

BobPang · 2023-04-11 23:29:54 · 个人记录

写在前面

建议不要看这篇文章，因为：

这是本juruo的第二篇blog
~~我自己还没学完呢~~
这原来只是我为了提升自己水平写的（众所周知给别人讲解可以吸收90%的芝士）

因此，请划走，不要浪费你宝贵的时间。

如果你真要看这篇文章，请先阅读这个，本文是它的续集。

线段与直线

直线有关

直线的正交判定

判断两条直线是否正交有很多方法，这里讲述在上一章的基础上最简单的方法。

首先，判断两条直线是否正交可以转化为判断两个向量是否正交。我们知道，\cos 90^\circ=0，且 \overrightarrow{a}\cdot\overrightarrow{b}=\left|\overrightarrow{a}\right|\left|\overrightarrow{b}\right|\cos\theta。由于 \left|\overrightarrow{a}\right|,\left|\overrightarrow{b}\right|\neq0，所以 \overrightarrow{a}\cdot\overrightarrow{b}=0\Longleftrightarrow\cos\theta=0\Longleftrightarrow\theta=90^\circ\Longleftrightarrow \overrightarrow{a},\,\overrightarrow{b} 正交。即，两个向量正交等价于它们的内积为 0 。

接下来思考如何实现把判断直线是否正交转化为判断向量是否正交。为了达到这个目标，我们需要先计算 \overrightarrow{AB},\ A(a_x,a_y),\,B(b_x,b_y) 。这很简单，只需将 B 向下平移 a_x，向左平移 a_y 即可。以下是代码实现：

Vector getVector(Point a,Point b)
{
    return makeVector(b.x-a.x,b.y-a.y);
}

接着是判断两个向量是否正交的函数：

bool ifOrthogonal(Vector a,Vector b)
{
    return equals(dot(a,b),0.0);
}

最后是实现直线与向量的转化并判断两直线是否正交的函数：

bool ifOrthogonal(Line a,Line b)
{
    Vector p,q;
    p=getVector(a.p1,a.p2);
    q=getVector(b.p1,b.p2);
    return ifOrthogonal(p,q);
}

直线的平行判定

直线的平行判定与正交判定差别很小。唯一的不同是，正交看的是向量的数量积，而平行则需要看向量积。\sin 0^\circ=\sin 180^\circ=0，且 \left|\overrightarrow{a}\times\overrightarrow{b}\right|=\left|\overrightarrow{a}\right|\left|\overrightarrow{b}\right|\sin\theta。由于 \left|\overrightarrow{a}\right|,\left|\overrightarrow{b}\right|\neq0，所以 \left|\overrightarrow{a}\times\overrightarrow{b}\right|=0\Longleftrightarrow\sin\theta=0\Longleftrightarrow\theta=0^\circ 或 180^\circ\Longleftrightarrow\overrightarrow{a},\overrightarrow{b} 平行。即，两个向量平行等价于它们的外积的模为 0 。

代码实现如下：

bool ifParallel(Vector a,Vector b)
{
    return equals(cross(a,b),0.0);
}
bool ifParallel(Line a,Line b)
{
    Vector p,q;
    p=getVector(a.p1,a.p2);
    q=getVector(b.p1,b.p2);
    return ifParallel(p,q);
}

映射

我们可以通过向量来较为简单地求点 P(p_x,p_y) 在直线 P_1P_2 上的映射。

首先，我们设向量 \overrightarrow{h}=\overrightarrow{P_1P},\,\overrightarrow{s}=\overrightarrow{P_1P_2},\,\left<\overrightarrow{h},\overrightarrow{s}\right>=\theta，点 P 在直线 P_1P_2 上的垂足为 x，P_1x=t 。

我们可以发现，t=\left|\overrightarrow{h}\right|\cos\theta 。又因 \overrightarrow{h}\cdot\overrightarrow{s}=\left|\overrightarrow{h}\right|\left|\overrightarrow{s}\right|\cos\theta，所以 t=\dfrac{\overrightarrow{h}\cdot\overrightarrow{s}}{\left|\overrightarrow{s}\right|} 。但是 t 只是一个长度，用 t 不能求出 x 的位置。因此我们要借助 \overrightarrow{s} ，再算出 t 与 \left|\overrightarrow{s}\right| 的比例，并用 \overrightarrow{s} 乘上这个比例，就能得到 \overrightarrow{P_1x} 从而进一步得到 x 。

按照上面的步骤，首先我们算出比例，\dfrac{t}{\left|\overrightarrow{s}\right|}=\dfrac{\overrightarrow{h}\cdot\overrightarrow{s}}{\left|\overrightarrow{s}\right|^2}，并发现 \overrightarrow{P_1x}=\overrightarrow{s}\dfrac{\overrightarrow{h}\cdot\overrightarrow{s}}{\left|\overrightarrow{s}\right|^2} 。因此，x=P_1+\overrightarrow{s}\dfrac{\overrightarrow{h}\cdot\overrightarrow{s}}{\left|\overrightarrow{s}\right|^2} 。

计算映射点的代码如下。

Point project(Line s,Point p)
{
    Vector P=makeVector(p.x,p.y),P1=makeVector(s.p1.x,s.p1.y),P2=makeVector(s.p2.x,s.p2.y);
    double ratio=dot(P-P1,P2-P1)/(P2-P1).norm();
    Vector x=P1+(P2-P1)*ratio;
    return makePoint(x.x,x.y);
}

映像

~~之前的东西都太难了，这次来点简单的，放松一下。~~

映像，就是指以直线 P_1P_2 为对称轴，翻转点 P 到 P' 。

看这个图，我们可以发现，求映像的方法很简单：先求出映射点，再把从点 P 到映射点的向量乘2，就结束了。代码实现如下：

Point reflect(Line s,Point p)
{
    Point pro=project(s,p);
    Vector Vpro=makeVector(pro.x,pro.y),P=makeVector(p.x,p.y);
    Vector x=P+(Vpro-P)*2;
    return makePoint(x.x,x.y);
}

线段有关

CCW

CCW(Counter_Clockwise)，实现很简单却非常重要。CCW可以求出三个点的位置关系。

具体地，你需要把 P_0,P_1,P_2 三个点分成以下五类：

方便起见，我们设 \overrightarrow{P_0P_1}=\overrightarrow{a},\,\overrightarrow{P_0P_2}=\overrightarrow{b} 。各种情况的判断方法如下。

外积的大小 \left|\overrightarrow{a}\times\overrightarrow{b}\right|>0 时，\overrightarrow{b} 位于 \overrightarrow{a} 的逆时针方向。
外积的大小 \left|\overrightarrow{a}\times\overrightarrow{b}\right|<0 时，\overrightarrow{b} 位于 \overrightarrow{a} 的顺时针方向。
不符合1，2时，说明 P_2 位于直线 P_0P_1 上。由于 \cos\theta 在 90^\circ<\theta\leq180^\circ 为负，因此在 \overrightarrow{a}\cdot\overrightarrow{b}<0 时，P_2 位于线段 P_0P_1 后方，即属于第3种情况。
不符合1，2，3时，说明 \overrightarrow{a} 与 \overrightarrow{b} 方向相同。此时只需比较它们的模，就能知道到底该属哪种情况了。当 \left|\overrightarrow{a}\right|<\left|\overrightarrow{b}\right| 时，属于第4种情况。
不符合1，2，3，4时，只能是第5种情况。

ccw() 函数可以返回三个点的位置。

int ccw(Point p0,Point p1,Point p2)
{
    Vector a=getVector(p0,p1),b=getVector(p0,p2);
    if(cross(a,b)>EPS)
        return 1;
    if(cross(a,b)<-EPS)
        return -1;
    if(dot(a,b)<-EPS)
        return 2;
    if(a.length()<b.length()-EPS)
        return -2;
    return 0;
}

判断线段相交

"线段相交"被定义为两条线段有重合部分。因此，两条线段共用一个端点、重合等情况均被视为相交。

给定4个点 A,B,C,D，容易想到判断线段 AB 与 CD 是否相交的方法如下：

线段 AB 与线段 CD 相交 \Longleftrightarrow 线段 AB 与直线 CD 相交 \&\& 直线 AB 与线段 CD 相交

如何判断线段和直线是否相交？如果 C、D 在直线 AB 同侧，则线段 CD 与直线 AB 不相交，反之亦然。这样就可以用ccw解决了：

bool intersect(Segment s1,Line s2)
{
    return (ccw(s2.p1,s2.p2,s1.p1)*ccw(s2.p1,s2.p2,s1.p2)<=0);
}

再回到原来的线段与线段的问题：

bool intersect(Segment s1,Segment s2)
{
    Line S1,S2;
    S1.p1=s1.p1;S1.p2=s1.p2;
    S2.p1=s2.p1;S2.p2=s2.p2;
    return (intersect(s1,S2) && intersect(s2,S1));
}

线段的交点

现在需要求线段 AB、线段 CD 的交点。简单起见，我们设 \overrightarrow{AB}=\overrightarrow{a},\,\overrightarrow{CD}=\overrightarrow{b},\,\overrightarrow{AC}=\overrightarrow{h},\,\overrightarrow{AD}=\overrightarrow{s}，AB 和 CD 交于 x。

首先，我们可以求出 C 在直线 AB 上的高 d1 。d1 也可以看成点 C 在 \overrightarrow{a},\overrightarrow{h} 所组成的平行四边形上的高。这样一来，就可以利用外积求出 d1：

d1=\dfrac{\left|\overrightarrow{a}\times\overrightarrow{h}\right|}{\left|\overrightarrow{a}\right|}

点 D 在直线 AB 上的高 d2 同理：

d2=\dfrac{\left|\overrightarrow{a}\times\overrightarrow{s}\right|}{\left|\overrightarrow{a}\right|}

接下来，就可以求出 Cx:xD=t=\dfrac{d1}{d1+d2}

所以交点 x=C+\overrightarrow{b}\times t。

计算线段交点的代码如下：

Point crossPoint(Segment s1,Segment s2)
{
    Point A=s1.p1,B=s1.p2,C=s2.p1,D=s2.p2;
    Vector a=getVector(A,B),b=getVector(C,D),h=getVector(A,C),s=getVector(A,D);
    double d1=abs(cross(a,h)),d2=abs(cross(a,s));
    double t=d1/(d1+d2);
    Vector ans=makeVector(C.x,C.y)+b*t;
    return makePoint(ans.x,ans.y);
}

注意上述代码中，d1 是 \left|\overrightarrow{a}\times\overrightarrow{h}\right| 而非 \dfrac{\left|\overrightarrow{a}\times\overrightarrow{h}\right|}{\left|\overrightarrow{a}\right|} 。这是因为我们只需要算出 t，而 t=\dfrac{d1}{d1+d2} ，d1、d1+d2 都是分母为 \left|\overrightarrow{a}\right| 的数，因此可以在计算 d1 时就把 \left|\overrightarrow{a}\right| 约去。d2 同理。

Update

2023-4-11 写完上一章，写完“直线的平行判定”。

2023-4-13 写完“映像”。

2023-4-14 写完“ccw”。