计算几何基础
判断点是否在线段上、判断两条线段是否相交
这里采用向量的解法。有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版的总结 */