美文网首页
计算几何基础

计算几何基础

作者: fruits_ | 来源:发表于2018-03-01 01:36 被阅读0次

    判断点是否在线段上、判断两条线段是否相交

    这里采用向量的解法。有2个概念:向量的内积和外积。内积又称为点积dot product,公式即 a·b = |a||b|cosΘ。

    两个向量a = [a1, a2,…, an]
    和b = [b1, b2,…, bn]的点积定义为:
    a·b=a1b1+a2b2+……+anbn。 
    

    外积又称为叉积,定义为 |a ×b| = |a|·|b|·sin<a,b> 在三维的向量中, 利用三阶行列式,写成


    在这里我们考虑平面的二维情况。对于向量a和b的内积:
    a · b > 0 说明向量a与b夹角范围[0, 90)
    a · b = 0 说明a与b正交
    a · b < 0 说明a与b夹角范围(90°, 180°]
    对于平面上三点A、B、C,有向量AB和向量AC
    AB × AC > 0 说明AC在AB的逆时针方向
    AB × AC < 0 说明AC在AB的顺时针方向
    AB × AC = 0 说明ABC三点共线,相互的方向不知

    对于二维向量p1=(x1,y1)和p2=(x2,y2),有内积p1·p2=x1x2+y1y2 有外积p1×p2=x1y2-y1x2
    判断点q是否在线段p1-p2上(注意 这里的p1-p2是真的点相减,从而形成一个以源点为一端 减的结果坐标为一端的点,这样一个坐标的信息足以描述一条线段):
    先利用外积 (p1-q)×(p2-q)=0判断q是否在 直线 p1-p2上,再用内积(p1-q)·(p2-q)<=0判断点q是否落在线段p1-p2上。前一个外积作用很显然,后一个主要是看点q在该线段是不是这样子的:
    [p1]------[q]-----[p2] 即q在p1 p2中,如果是这样,那么内积就该因为夹角是180度而为负,如果q与这两个端点之一重合,那么内积又为0,所以,判断条件是 " <= 0"

    判断两线段相交

    思路是看两 直线 的交点是否在两条线段上,有变量t,直线p1-p2上的点可表示为p1+t(p2-p1) ,交点也在q1-q2上,所以根据该交点和q1 q2共线有
    (q2-q1) × (p1 + t(p2-p1)-q1) = 0
    由上面的公式可以把t给表示出来,再把交点代入①式有
    p1 + t(p2-p1)
    = p1 + (p2-p1) * ((q2-q1)×(q1-p1)) / ((q2-q1)×(p2-p1))
    另外由于当两线段平行时,虽然对应的直线没有交点,但是可能有公共点,这个可以通过看线段的端点是否在另一条线段上来判断。
    得模板

    const double EPS = 1e-10;
    
    double add(double a, double b) {
        if (abs(a + b) < EPS * (abs(a) + abs(b))) return 0;
        return a + b;
    }
    
    struct P {
        double x, y;
        P(double _x = 0, double _y = 0) : x(_x), y(_y) {}
        P operator + (P p) { return P(add(x, p.x), add(y, p.y)); }
        P operator - (P p) { return P(add(x, -p.x), add(y, -p.y)); }
        P operator * (double d) { return P(x * d, y * d); }
        double dot(P p) { return add(x * p.x, y * p.y); }
        double det(P p) { return add(x * p.y, -y * p.x); }
    };
    
    // 判断p是否在线段p1-p2上
    bool on_seg(P p1, P p2, P q) { return (p1 - q).det(p2 - q) == 0 && (p1 - q).dot(p2 - q) <= 0; }
    // 计算直线p1-p2与直线q1-q2的交点 直线!
    P intersection(P p1, P p2, P q1, P q2) {
        return p1 + (p2 - p1) * ((q2 - q1).det(q1 - p1) / (q2 - q1).det(p2 - p1));
    }
    
    凸包问题

    什么是凸包? 解释在这
    这里介绍基于平面扫描的Graham扫描算法。首先把点集按x->y坐标的字典序升序排列,排序后的第一个点和最后一个点必然在凸包点集之中,它们之间的部分分成上下两条链分别求解。在构造过程中的凸包末尾加上新顶点后,可能会破坏凸性,要做的是把凸的部分的点从末尾去除。这里给出书上的一个图帮助理解:


    那怎么判断新加的点会不会破坏凸性呢?利用内积。已知内积利用右手定则可以判断方向的正负。
    [假设内积是c,则若坐标系是满足右手定则的,当右手的四指从a以不超过180度的转角转向b时,竖起的大拇指指向是c的方向] 或者干脆知道上面总结的,对于直线关于内积正负的相对位置关系。
    已经在凸包序列中的,按先后顺序为qs[k - 2], qs[k - 1], 现在要加入ps[i]

    以构造下链为例,qs[k-2], qs[k-1]已经组成了凸包的一部分,现在考虑加入ps[i],但是发现加入ps[i]后凸性被破坏,这里认为被破坏,即点qs[k-1]凹了下去,亦即向量(qs[k-1],qs[k-2])和向量(qs[k-1],ps[i])的内积>=0,即后者在前者的逆时针方向, 同理考虑构造上链计算和判别是完全一样的。那么结合上面点的结构体可得计算凸包的代码:
    // 字典序比较
    bool cmp_x(const P& p, const P& q) {
        if (p.x != q.x) return p.x < q.x;
        return p.y < q.y;
    }
    // 求凸包 其中ps是所有点的集合
    vector<P> convex_hull(P* ps, int n) {
        sort(ps, ps + n, cmp_x);
        int k = 0;              // 凸包顶点数量k
        vector<P> qs(2 * n);    // 构造中的凸包点集
        // 构造凸包的下侧
        for (int i = 0; i < n; ++i) {
            while (k > 1 && (qs[k - 1] - qs[k - 2]).det(qs[k - 1] - ps[i]) >= 0) --k;
            qs[k++] = ps[i];
        }
        // 构造凸包上侧
        for (int i = n - 2, t = k; i >= 0; --i) {
            while (k > t && (qs[k - 1] - qs[k - 2]).det(qs[k - 1] - ps[i]) >= 0) --k;
            qs[k++] = ps[i];
        }
        qs.resize(k - 1);
        return qs;
    }
    

    /* 参考挑战程序设计竞赛第2版的总结 */

    相关文章

      网友评论

          本文标题:计算几何基础

          本文链接:https://www.haomeiwen.com/subject/ospwxftx.html