【笔记篇】最良心的计算几何学习笔记（一）

世界以痛吻我，我却报之以歌。

开新坑...
虽然不知道这坑要填多久...
文章同步上传到github...
有想看的可以去看看→_→

*温馨提示:

看本文之前请务必学习或回顾数学-必修2的解析几何和数学-必修4的平面向量有关内容...
本文中的代码基本都是口胡的, 不过还是用了IDE, 只保证能过编译, 不保证正确OvO...

一浮点数相关...

众所周知, 浮点数的运算是有误差的...所以有时候会出现一些很蛋疼的事情.. 比如两个该相等的东西会得出不相等的结果... 为了避免这种事情的发生, 我们要允许一些误差... 我们定义

const double eps=1e-9;

当然这里如果要求精度的话也可以开long double... 而且这个1e-9也是不一定的...
有些题目开1e-7 1e-8 1e-10都是不一定的... (诅咒所有卡精度的出题人吃方便面只有调料包)
这个eps有什么用呢...
这样我们就可以自定义一下实数的比较...
然而double作为基本变量类型是不能重载==运算符的, 这样就要写个函数, 而且像作者这么zz的人很可能会忘记调用OvO

inline int dcmp(const double &a){
  	if(fabs(a)<eps) return 0;   //绝对值不超过eps即视为0(卡精度的来源)
	return a<0?-1:1;			//a<0 返回-1 a>0 返回1
}

这样我们关于double的东西就基本写完了..

二平面向量相关

看数学-必修4了没? 没有快去看...
根据平面向量基本定理, 每个向量(vec a)都可以写成(x*vec{e_1}+y*vec{e_2})的形式...
这里的(vec{e_1},vec{e_2})可以是任意的, 但是为了方便我们通常取正交基底(就是(vec{e_1}perpvec{e_2})啦~)
每个向量都可以由数对((x,y))唯一确定...
这样我们就可以这么定义一个向量:

struct vec{
	double x,y;
};

很显然, 这个也可以表示一个点.
所以我们可以加一句这个, 当然也可以不加..

typedef vec point;

在此基础上, 我们就来定义各类运算.
1.为了方便, 我们可以在vec结构体里加一个构造函数...

struct vec{
    double x,y;
    vec(int p,int q):x(p),y(q){}
    vec(){x=y=0.0;} //默认构造函数还是要有的..
};

2.向量的模?
表示为(|vec v|), 就是向量的长度啦~

inline double len(const vec &a){
    return sqrt(a.x*a.x+a.y*a.y);
}

3.比较基础的, 向量的四则运算.
这个由定义显然

vec operator +(const vec &a,const vec &b){
    return vec(a.x+b.x,a.y+b.y);
}
vec operator -(const vec &a,const vec &b){
    return vec(a.x-b.x,a.y-b.y);
}
vec operator *(const vec &a,const double &b){
    return vec(a.x*b,a.y*b);
}
vec operator *(const double &b,const vec &a){ //变量类型不同, 不能直接用交换律..
    return vec(a.x*b,a.y*b);
}
vec operator /(const vec &a,const double &b){
    return vec(a.x/b,a.y/b);
}

4.然后就是点积(数量积 ~~内积标积~~)和叉积(向量积 ~~外积矢积~~)了..
不过似乎并不考虑叉积的方向..(大概是因为二维计算几何的缘故吧..)
(不过好像叉积的定义是"矢量叉积定义为由(0,0)、p1、p2和p1+p2所组成的平行四边形的带符号的面积"?)

double operator ^(const vec &a,const vec &b){ //点积就用^就行了 反正应该是用不到异或的..
    return a.x*b.x+a.y*b.y;
}
double operator *(const vec &a,const vec &b){
    return a.x*b.y-a.y*b.x;
}

点积的话根据必修四, 可以判断垂直.. 若(vec pcdotvec q=0), 则(vec pperpvec q).

叉积的一个重要的作用就是判断向量的顺逆时针关系..

若(vec p imesvec q>0), 则(vec p)在(vec q)的顺时针方向;
若(vec p imesvec q<0), 则(vec p)在(vec q)的逆时针方向;
若(vec p imesvec q=0), 则(vec pparallelvec q), 但可能同向也可能反向.
有一些小小的注意事项(虽然不一定有用)..
据我猜测计算几何中没什么用的(但数学考试有用啊←_←)
(vec acdotvec b=|a||b|cos heta, vec a imesvec b=|a||b|sin heta)
点积是满足交换律的(显然), 而叉积是不满足交换律的(也显然), 不过叉积满足反交换律,
即(vec a imesvec b=-vec b imesvec a)
点积和叉积都满足加法的分配律.
显然不能连续点积, 叉积不满足结合律.
其他的似乎更没啥用而且我也看不懂想研究的可以点上方传送门去baidu看..

三 [跟多边形没关系的]各种常见的简单的判断

1.两点之间距离公式先放在这里

inline double dis(const vec &a,const vec &b){
    return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}

其实写多了之后可能会更喜欢写len(a-b)...

2.两个向量的顺逆时针关系见↑↑↑

3.折线段的拐向判断
其实跟上面没有太大区别..
两条有公共顶点线段(p_0p_1)和(p_1p_2)(因为是折线啊),
很显然我们可以表示两个向量(vec{v_1}=p_2-p_0, vec{v_2}=p_1-p_0)然后求(vec{v_1} imesvec{v_2})...

若(vec{v_1} imesvec{v_2}>0), 则折线往右拐;
若(vec{v_1} imesvec{v_2}<0), 则折线往左拐;
若(vec{v_1} imesvec{v_2}=0), 则(p_0,p_1,p_2)三点共线.

4.判断点(q)是否在以(p_0,p_1)为对角顶点的矩形(不加说明则默认矩形边与坐标轴平行)内
这个简单啊只要分别判断(q)的横纵坐标比(p_0,p_1)的最小值大, 最大值小(当然是可以相等的)即可..

inline bool ptInRect(const vec &q,const vec &p0,const vec &p1){
    return min(p0.x,p1.x)<=q.x&&q.x<=max(p0.y,p1.y)&&min(p0.y,p1.y)<=q.y&&q.y<=max(p0.y,p1.y);
}

5.判断点(q)是否在线(p0p1)上
根据2. , 很显然的可以判断((q-p_0) imes(p_1-p_0)==0)
但这样就够了吗? 并不是... 因为(q)可能在延长线or反向延长线上...
所以我们要保证(q)在以(p_0,p_1)为对角顶点的矩形内...(你以为我为什么要写3. ?)

inline bool ptOnSeg(const vec &q,const vec &p0,const vec &p1){
    return (p1-p0)*q==0&&ptInRect(q,p0,p1);
}

这样就可以咯..

6.判断线段(p_0p_1)和直线(q_0q_1)是否相交
如果线段与直线相交, 那么线段跨立直线..

首先如图,我们可以得出(vec{v_1}=p_0-q_0,vec{v_2}=q_1-q_0,vec{v_3}=p_1-q_0).
由图可以看出, 我们要让(p_0)和(p_1)放到直线(q_0q_1)的两侧, 就要让(vec{v_1} imesvec{v_2})与(vec{v_2} imesvec{v_3})同号(或者有一方为0)
所以有((vec{v_1} imesvec{v_2})*(vec{v_2} imesvec{v_3})geqslant0), 也就是(((p_0-q_0) imes(q_1-q_0))*((q_1-q_0) imes(p_1-q_0))geqslant0).

inline bool segCutLine(const vec &p0,const vec &p1,const vec &q0,const vec &q1){
    return ((p0-q0)*(q1-q0))*((q1-q0)*(p1-q0))>=0;
}

7.判断两条线段(p_0p_1,q_0q_1)是否相交

这个我们可以分为两步

快速排斥实验
两个线段的端点为对角顶点的两个矩形如果不相交, 那么这两条线段显然不相交.
跨立实验
明确一点: 如果两线段相交, 则两线段必然相互跨立对方.
我们就可以根据上面的跨立方式得出要判断(((p_0-q_0) imes(q_1-q_0))*((q_1-q_0) imes(p_1-q_0))>0 (p_0p_1跨立q_0q_1))和(((q_0-p_0) imes(p_1-q_0))*((p_1-q_0) imes(q_1-p_0))>0 (q_0q_1跨立p_0p_1))
而当前面这一坨等于0的时候呢? 我们可以得知共线, 而由于已经通过了快速排斥实验, 所以端点一定在另一条线段上. 所以我们只要让(((p_0-q_0) imes(q_1-q_0))*((q_1-q_0) imes(p_1-q_0))geqslant0)和(((q_0-p_0) imes(p_1-q_0))*((p_1-q_0) imes(q_1-p_0))geqslant0)就可以了 (其实好啰嗦啊)

inline bool segCutSeg(const vec &p0,const vec &p1,const vec &q0,const vec &q1){
    if(min(p0.x,p1.x)>max(q0.x,q1.x)||min(q0.x,q0.y)>max(p0.x,p1.x)
        ||min(p0.y,p1.y)>max(q0.y,q1.y)||min(q0.y,q1.y)>max(p0.y,p1.y)) return false; //快速排斥实验
    return ((p0-q0)*(q1-q0))*((q1-q0)*(p1-q0))>=0&&((q0-p0)*(p1-p0))*((p1-p0)*(q1-p0))>=0; //跨立实验
}

8.判断矩形是否在矩形中
只需要简单的比较边界就行了.

inline bool rectInRect(const vec &p0,const vec &p1,const vec &q0,const vec &q1){
    return min(p0.x,p1.x)<min(q0.x,q1.x)&&max(p0.x,p1.x)>max(q0.x,q1.x)&&min(p0.y,p1.y)<min(q0.y,q1.y)&&max(p0.y,p1.y)>max(q0.y,q1.y);
}

9.判断圆是否在矩形中
圆在矩形中的充要条件是圆心在矩形内而且圆的半径小于到矩形四边距离的最小值.
话说圆要怎么表示? 一个点表示圆心一个double表示半径嘛= =

inline bool cirInRect(const vec &o,const double &r,const vec &p0,const vec &p1){
    return ptInRect(o,p0,p1)&&min(
        min(abs(min(p0.x,p1.x)-o.x),abs(max(p0.x,p1.x)-o.x)),
        min(abs(min(p0.y,p1.y)-o.y),abs(max(p0.y,p1.y)-o.y)))<r; //这么写好像清楚一点?
}

10.判断点是否在圆中
与圆心的距离小于半径啊~

inline bool ptInCir(const vec &p,const vec &o,const double &r){
    return dis(p,o)<r;
}

11.判断线段、折线、多边形是否在圆中
圆是个凸集, 所以我们依次判断每个端点是否在圆内即可..
代码就不给了OvO

12.判断圆(O_1)是否在圆(O_2)内
首先嘛肯定要(r_1<r_2)... 其次呢要有(|O_1O_2|<=r_2-r_1)(其实好像可以和上一条一起判.. 因为负数显然不会大于距离..) (等号成立时两圆内切(突然不知道算不算在圆内了..先算是吧...~~(大不了到时候抠掉嘛)~~))

inline bool cirInCir(const vec &o1,const double &r1,const vec &o2,const double &r2){
    return dis(o1,o2)<=r2-r1;
}

13.点(p)到直线(q0q1)的距离,利用叉积求出平行四边形的面积, 然后算出底边长. 一比得到高即可.

inline bool ptDisLine(const vec &p,const vec &q0,const vec &q1){
    return (p-q0)*(q1-q0)/len(q1-q0);
}

14.点(p)到线段(q0q1)的距离,先判断垂足是否在线段上.
是则同上, 不是则找距离两个端点中近的那个
如何判断呢?
垂足不在线段上
如图所示, 如果垂足在直线外, 那么夹角( heta)一定是钝角, 则(cos heta<0)
那么(vec{e_1}cdotvec{e_2}=|vec{e_1}||vec{e_2}|cos heta<0),
而且这样求出的( heta)的顶点就是较近的点..
两边一求就好了...

double ptDisSeg(const point &p,const point &a,const point &b){
    if(dcmp((p-a)^(b-a))<0) return dist(p,a);
    if(dcmp((p-b)^(a-b))<0) return dist(p,b);
    return fabs((p-a)*(b-a))/dist(b,a);
}

15.两条相交直线(线段)的交点. (所以要先判断是否相交...)
好像两种方法?
一种是利用解析几何推导.
对于直线(l:Ax+By+C=0)和直线上的两点((x_0,y_0),(x_1,y_1))
斜率(k=-frac{A}{B}=frac{y_0-y_1}{x_0-x_1}),
那我们不妨设(A=y_0-y_1,B=x_1-x_0), 代入方程得到(C=x_0y_1-x_1y_0)
然后回到问题, 设两条直线分别是

[ l_1:a_1x+b_1y+c_1=0,\ l_2:a_1x+b_1y+c_2=0. ]

 这样(x,y)就是联立得到方程组的解.
 根据相关知识, 我们可以得到

[ left{egin{matrix} x=frac{b_2c_1-b_1c_2}{b_1a_2-b_2a_1}\ y=frac{a_1c_2-a_2c_1}{b_1a_2-b_2a_1} end{matrix} ight. ]

 而这些量都可以用叉积的形式表示出来..
 所以就可以这么写:

inline vec lineCutLineNode(const vec &p0,const vec &p1,const vec &q0,const vec &q1){
    double a1,b1,c1,a2,b2,c2,d;
    a1=p1.y-p0.y; b1=p0.x-p1.x; c1=p0*p1;
    a2=q1.y-q0.y; b2=q0.x-q1.x; c2=q0*q1;
    d=a1*b2-a2*b1;
    return vec((b2*c1-b1*c2)/d,(a1*c2-a2*c1)/d);
}

这种方法略微麻烦, 但是精度比较好.
另一种方法是利用叉积的比值.
利用参数方程的知识(数学选修2-1 今天刚学哒XD), 我们可以用一个定点和一个方向向量确定一条直线.
令两条直线

[l1=P+tvec v, \ l2=Q+tvec w ]

然后因为是交点, 所以同时满足两条直线的方程, 令交点在两条直线的参数分别为(lambda)和(mu), 则有

[ P+lambdavec v=Q+muvec w ]

两边分别叉乘(vec w)

[ (P+lambdavec v) imesvec w=(Q+muvec w) imes w ]

去括号,

[P imesvec w+lambdavec v imesvec w=Q imesvec w+muvec w imesvec w ]

又(vec w)一定与自己共线, 所以(vec w imesvec w=0), 这样就没有mu了, 移项, 合并同类项,

[(Q-P) imesvec w=lambdavec v imesvec w ]

系数除到另一边去, 得

[lambda=frac{(Q-P) imesvec w}{vec v imesvec w} ]

那么这样就都是已知量了, 代入进去即可.
最后要求的交点就是(P+lambdavec v)啦~

留一道例题:poj1269

完整的代码在github里面哟~