计算几何基础

2018-03-01  本文已影响0人  fruits_

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

这里采用向量的解法。有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版的总结 */

上一篇下一篇

猜你喜欢

热点阅读