计算几何初步西电acm笔记1

Ridiculous · 2020-07-14 19:50:24 · 个人记录

计算几何初步

点积

大家高中学的一定都很好，不在赘述。

叉积

叉积：

\vec a=(x_0,y_0),\vec b=(x_1,y_1)

cross=\vec a \times \vec b=|a|*|b|*sin\theta=x0y1-y0x1

叉积的简单用处:

通过叉积的正负判断两个向量的方向。cross>eps则，\vec b位于\vec a逆时针方向，cross<-eps位于顺时针方向，-eps<cross<eps则两向量共线。
求三角形面积。S=|cross|\div2
求简单多边形的有向面积：对于多边形ABCDE...（字母为简单多边形的顶点，按照顺时针或逆时针的方向排序）。有S=(\overrightarrow{OA}\times\overrightarrow{OB}+\overrightarrow{OB}\times\overrightarrow{OC}...)\div2

有向面积指计算的过程中叉积出的结果有正有负。当多边形的顶点按照逆时针排序时，计算出的面积为正，当多边形的顶点按照顺时针排序时，计算出的面积为负。手动简单模拟可知其原理，并且该公式对凸多边形与凹多边形均适用。
与点积配合求向量夹角。有\theta=abs(atan2(\vec a\times \vec b,\vec a\cdot\vec b))

题目：圆的k次交，hihocoder 1429：

图片来自adkdycoin神犇

将相交的区间提取出来，根据区间被覆盖的次数来进行计算。

圆的1次交的情况

r^2\times (\theta-sin\theta)\div2

圆的2次交的情况

线段交

具体内容移步网络资源。

static bool is_intersected(const Segment2& seg1, const Segment2& seg2)//快速排斥实验+跨立实验
    {
        if (max(seg1.st.x, seg1.ed.x) >= min(seg2.st.x, seg2.ed.x) &&
            max(seg2.st.x, seg2.ed.x) >= min(seg1.st.x, seg1.ed.x) &&
            max(seg1.st.y, seg1.ed.y) >= min(seg2.st.y, seg2.ed.y) &&
            max(seg2.st.y, seg2.ed.y) >= min(seg1.st.y, seg2.ed.y) &&
            Vector2(seg1.st, seg2.st).cross(Vector2(seg1.st, seg1.ed)) * Vector2(seg1.st, seg1.ed).cross(Vector2(seg1.st, seg2.ed)) >= -eps &&
            Vector2(seg2.st, seg1.st).cross(Vector2(seg2.st, seg2.ed)) * Vector2(seg2.st, seg2.ed).cross(Vector2(seg2.st, seg1.ed)) >= -eps)
            return true;
        return false;
    }
    static Point2 get_intersection(const Segment2& seg1, const Segment2& seg2)//使用前判断是否相交,使用相似来证明
    {
        double S1 = Vector2(seg1.st, seg1.ed).cross_abs(Vector2(seg1.st, seg2.ed));
        double S2 = Vector2(seg1.st, seg1.ed).cross_abs(Vector2(seg1.st, seg2.st));
        Point2 ret;
        ret.x = (S1 * seg2.st.x + S2 * seg2.ed.x) / (S1 + S2);
        ret.y = (S1 * seg2.st.y + S2 * seg2.ed.y) / (S1 + S2);
        return ret;
    }

题目：UVA393

https://www.luogu.com.cn/problem/UVA393

极角排序

按照极角进行排序，从而让我们可以对其进行各种操作（扫描，dp，分治等）

https://www.luogu.com.cn/problem/UVA11696

UVA11696，题目大意：给定一些点和圆，保证圆不相交，点不落在圆内。如果两个点的连线没有经过圆，则两点连通。求连通块的个数。

二维凸包

简单的说，在平面上给出一些点，从这些点中选出一些点构成一个凸多边形，使得所有点都落在凸多边形以内。则此凸多边形被称为凸包。

凸包可以在O(nlogn)的时间内求出。求凸包的算法有graham扫描法和Andrew算法。前者使用极角排序，后者使用坐标排序。对于Andrew算法来说：

首先把所有点以坐标y为第一关键字，坐标x为第二关键字，进行排序。
之后进行两次扫描，每一次扫描会得到凸包的一半。我们使用栈来保存当前状态找到的凸包。我们按照排序的顺序依次访问所有点。
初始状态我们取两个点加入栈内。
设我们逆时针生成一个凸包，当前栈顶两个点分别为A,B(B为栈顶)，新加入的点为C。我们通过判断\overrightarrow {AB}\times \overrightarrow{BC}的正负来判断加入当前点后栈中的点围成的图形是否为凸。如果加入点会使其变凹，则移除栈顶点，并再次判断。若栈内点的个数不足两个，不再进行判断（新加入的点一定不会构成非凸的图形)。访问所有点，我们可得到右侧的一个凸包。
以排序后的最后一个点为起点，按照排序后点集中相反的顺序遍历，重复4中的操作，得到左侧凸包。
合并两个凸包，去掉最下方的多余点。

vector<Point2> Andrew(vector<Point2> p, string mode = "default")//若排除三点共线（旋转卡壳）使用mode special
{
    if (p.size() <= 2)
        return p;
    double mode_num = 0;
    if (mode == "special")
        mode_num = -eps;
    sort(p.begin(), p.end(), [](const Point2& a, const Point2& b) {return a.y < b.y || (a.y == b.y && a.x < b.x); });
    int top = -1, tmp;
    vector<Point2> st;
    top = 1;
    st.push_back(p[0]), st.push_back(p[1]);
    for (int i = 2; i < p.size(); i++)
    {
        double t = Vector2(st[top - 1], st[top]).cross(Vector2(st[top], p[i]));
        while (top > 0 && Vector2(st[top - 1], st[top]).cross(Vector2(st[top], p[i])) > mode_num)
        {
            top--;
            st.pop_back();
        }
        top++;
        st.push_back(p[i]);
    }
    tmp = top;
    top++;
    st.push_back(p[p.size() - 2]);
    for (int i = p.size() - 3; i >= 0; i--)
    {
        while (top > tmp && Vector2(st[top - 1], st[top]).cross(Vector2(st[top], p[i])) > mode_num)
        {
            top--;
            st.pop_back();
        }
        top++;
        st.push_back(p[i]);
    }
    st.pop_back();
    return st;
}

poj1696 :space ant

平面上给定50个点，是最左侧的某个点作为起点，蚂蚁在点与点之间直线移动，且路径只能沿逆时针方向且路径不能相交

It can not turn right due to its special body structure.
It leaves a red path while walking.
It hates to pass over a previously red colored path, and never does that.

贪心求以当前点为中心，极角最小的点，选择其为下一个移动的点。

实际上我们构造了一个又一个凸包，而这能覆盖到所有的点。

图片引自https://blog.csdn.net/u013480600/article/details/40125055

poj1259 最大空凸包

dp[i][j]=area(O,i,j)+max(dp[j][k]),\overrightarrow {kj}\times\overrightarrow{ji}>0

表示以ji,iO为凸包最后两条边时，凸包的最大面积。

考虑优化，优化转移过程。记录下每次转移过程中需要用到的信息。

g[i][j]=max(dp[i][k],k\in[0,j])

t,dp[i][j]=area(O,i,j)+g[j][k]

特判Oij共线的情况，此时不更新dp数组。

半平面交

半平面交，故名思意，其为一些被直线分隔而成的半平面的交集。

对于一系列直线分隔而成的半平面，其描述如下：

\begin{equation} \left\{ \begin{array}{l} A_1x+B_1y+C_1\ge0\\ A_2x+B_2y+C_2\ge0\\ ......\\ A_nx+B_ny+C_n\ge0 \end{array} \right. \end{equation}

符合上述不等式组的所有点都是半平面交得到的集合。

注：以下所有图片均来自于oiwiki

通常情况下，竞赛题目中的半平面交的面积往往存在，此时直线围成的部分构成了一个凸多边形。

很自然的可以想到，我们可以用类似于求凸包的方式来处理半平面交，同时我们令半平面均为有向直线（向量）的左侧。

首先我们对直线按照逆时针排序，我们使用双端队列来保存当前的直线和交点。

半平面交的求取过程中会遇到以下几种情况。

case 1：队尾点被当前直线约束，该种情况类似于graham扫描法求取凸包的情况，弹出队尾点与队尾直线。

case 2：新加入直线对队首点产生约束。

case 3：处理完所有直线后，判断队首直线对队尾点是否产生约束。实际上，直线之间是互相约束的（也就是半平面的交集），而我们每次加入直线时都考虑的是该直线是否会约束前面加入的直线，而没有考虑前面加入的直线是否约束该直线。通过这一步，我们将去掉被前面直线所约束的队尾直线。

case 1与case 2的判断顺序不能调换，否则会出现：

代码：

vector<Line2> halfplane_intersect(vector<Line2>line)
{
    //判断是否被约束，即点是否在直线的右侧。
    auto check = [&](const Point2& p, const Line2& a)
    {
        if (Vector2(a.st, p).cross(Vector2(a.st, a.ed)) > 0)
            return true;
        return false;
    };
    //极角排序
    sort(line.begin(), line.end());
    deque<Line2>que;
    deque<Point2>ans;
    que.push_back(line[0]);
    for (int i = 1; i < line.size(); i++)
    {
        if (abs(line[i].A * line[i - 1].B - line[i].B * line[i - 1].A) < eps)//直线平行时取靠左侧的直线
            continue;
        while (!ans.empty() && check(ans.back(), line[i]))//判断是否会和最后加入的几个点发生冲突
        {
            ans.pop_back();
            que.pop_back();
        }
        while (!ans.empty() && check(ans.front(), line[i]))//判断是否会和最开始加入的点冲突
        {
            ans.pop_front();
            que.pop_front();
        }
        Point2 p = Line2::get_intersection(line[i], que.back());
        ans.push_back(p);
        que.push_back(line[i]);
    }
    while (!ans.empty() && check(ans.back(), que.front()))//判断最后加入的点是否会有不合法现象
    {
        ans.pop_back();
        que.pop_back();
    }
    ans.push_back(Line2::get_intersection(que.back(), que.front()));
    return vector<Line2>(que.begin(), que.end());
}

题目：

https://www.luogu.com.cn/problem/UVA1396

UVA1396

http://acm.pku.edu.cn/JudgeOnline/problem?id=3384

poj 3384 Feng Shui

用两个大小相等的圆覆盖一个多边形（圆必须在多边形内，可以相交，问最多能覆盖多边形的面积。

旋转卡壳

如何求得一个凸包的直径（最远点对的距离）？

旋转卡壳！

旋转卡壳用两条平行线夹住我们的凸包，其中一条线与凸包的边相交，另一条线卡在凸包的另一侧。旋转这两条平行线的同时并计算距离，就可以在O(n)的时间内求出凸包的直径。

下图转自https://www.jvruo.com/archives/79/对旋转卡壳的讲解。

当与边重合的直线确定后，我们可以发现，另一条直线经过的点是离该直线最远的点。同时，因为此多边形是一个凸包，点到直线的距离是一个凸函数，这使得我们可以顺着直线的旋转方向枚举另外一个点，时间复杂度O(n)

代码：

double rotate_calipers(vector<Point2>p)
{
    double ans = INT64_MIN;
    int q = 1;
    //q为对踵点
    //i为枚举的边上的点
    for (int i = 0; i < p.size() - 1; i++)
    {
        while (Triangle(p[q], p[i + 1], p[i]).area < Triangle(p[(q + 1) % p.size()], p[i + 1], p[i]).area)
            q = (q + 1) % p.size();
        ans = max(ans, max(Point2::get_dis(p[i], p[q]), Point2::get_dis(p[i + 1], p[q])));
    }
    return ans;
}

假题

有些题目看起来像计算几何，实际上并没有什么关系。

https://www.luogu.com.cn/problem/UVA1411

UVA1411 Ants

题意：给定一些黑点白点(n\le100)，要求一个黑点连接一个白点，并且所有线段都不相交

杂项

以下内容不一定对所有人适用，请酌情参考，即使我再极力注意这些情况后也几乎每次都会犯下列错误的其中一种，在此将其列出。

精度控制。因为浮点数计算带来的误差，计算几何需要考虑到精度问题。通常来讲，多次使用三角函数（三角函数的导数在某些位置为0）与一些奇怪的除法运算（除一个特别小的数或特别大的数)，会引起精度误差。在做一些题目时，往往需要考虑精度问题。同时，有的题目也会给出精度判别，如点到圆心的距离与半径的长度之差不大于1e-6则判断为点在圆上。此时，需要根据具体情况调整程序中eps的大小。
特判情况。特判会给程序带来巨大的不稳定性。枚举各种情况往往是一件心烦且容易出错的事情。然而在计算几何的题目的正解中也往往会充斥着大量的特判。常见的直线无斜率，角度跨边界([-\pi,\pi]或[0,2\pi])等。需要时刻留意边界条件。
算法的正确性。很多题目需要你在草稿纸上画图来思考问题，这时需要注意你所思考的情况是否包含到题目的所有情况，你所设计的算法是否能够对每一种情况都适用。可以试一下POJ 2826，著名讨论题。
注意类型。请仔细检查编译器的警告选项，当其给出类似于：“从double转换为int，可能丢失数据”等类似的警告时，请仔细检查代码，判断是否存在强制转换带来的bug。类似的，当你做除法时，请仔细注意是否存在类似于下方的情况。
```
double d;
int i0,i1;
(i0/i1)*d;//你所写的
(1.0*i0/i1)*d//你想要的
```
注意板子的正确性。如果你和我一样，觉得自己用的代码一定要自己造轮子。那么伴随着这样的做法，你将经历~~我已经经历了~~一系列WA，RE，TLE等。当你过了一道题时，往往并不能验证板子的通用性。因为精度，题目条件等诸多原因，板子的鲁棒性并不能很好的得到验证。同时自己写的板子很可能伴随着大常数，代码啰嗦等问题。请通过大量的做题~~然而我并没有做到~~来验证自己所写代码的正确性。

刷题清单

看神仙的题单吧

https://blog.csdn.net/qq_34731703/article/details/53769806

想不开的同学

网上有很多计算几何的资源，也有很多模板。如果你想不开或懒于搜索，下方是我自己的代码。很多函数甚至尚未验证，即使验证了也有可能出现上述描述的问题，请谨慎取用。

https://paste.ubuntu.com/p/xDBSmSbVyB/

计算几何初步 西电acm笔记1