当前位置:   article > 正文

计算几何基础

icpc计算几何基础

1.基础的基础


1.点与向量.

点与向量都用x,y两个坐标表示,所以我们在程序上暂时不做区分.

但是注意:

  (1).点+向量=点;

  (2).向量+向量=向量(向量-向量=向量+(相反向量));

  (3).点-点=向量;

  (4).点+点没有意义;

其中pt(point)表示点,vt(vector)表示向量.

1 struct pt{
2     double x,y;
3     point(double x=0,double y=0):x(x),y(y){}
4 };
5 typedef pt vt;
6 vt operator + (vt a,vt b){ return vt(a.x+b.x,a.y+b.y); }
7 vt operator - (pt a,pt b){ return vt(a.x-b.x,a.y-b.y); }
8 vt operator * (vt a,double x){ return vt(a.x*x,a.y*y); }
9 vt operator / (vt a,double x){ return vt(a.x/x,a.y/x); }
View Code

 

2.点的比较函数

在有些需要排序的时候会用到.

1 bool operator < (const pt &a,const pt &b){ return a.x<b.x||(a.x==b.x&&a.y<b.y); }
View Code

 

3.判断点相等的函数和避免误差的三态函数

浮点误差总是存在的,所以我们用三态函数dcmp来减少精度带来的问题.

1 const double eps=1e-8;
2 int dcmp(double x){ if(fabs(x)<eps) return 0; else return x<0?-1:1; }
3 bool operator == (const pt &a,const pt &b){  return dcmp(a.x-b.x)==0&&dcmp(a.y-b.y)==0; }
View Code

 

4.向量的极角

向量的极角表示从x轴正半轴旋转到该向量所需要的角度.

使用C标准库里面的atan2函数即可,单位是弧度.

1 atan2(y,x);
View Code

 

 

2.基本运算


1.向量的点积与叉积

cross(a,b)=a.xb.ya.yb.x

这个是怎么来的呢?我们比较两个向量所在直线的斜率,如果b向量所在直线的斜率大(两个向量的斜率都存在的情况),那么就有

b.yb.x>a.ya.x

然后去分母

a.xb.y>a.yb.x

我们注意到,此时b向量所在的直线位于a向量所在直线的右手螺旋方向(a在b的左手螺旋方向).如果是小于号即b向量所在直线的斜率小,那么b向量所在的直线位于a向量所在直线的左手螺旋方向(a在b的右手螺旋方向).如果是等号即斜率一样,那么两条直线平行即向量a,b共线.

但这是两直线斜率都存在的情况.我们来讨论一下斜率不存在的情况:

(1).a所在直线的斜率不存在即a.x=0.

若原不等式仍成立即

0>a.yb.x
,再讨论a.y的正负,可以证明b所在直线仍然位于a所在直线的右手螺旋方向.

(2).b所在直线的斜率不存在即b.x=0.

若原不等式仍成立即

a.xb.y>0
,再讨论b.y的正负,可以证明b所在直线仍然位于a所在直线的右手螺旋方向.

(3).a,b所在直线的斜率都不存在即a.x=b.x=0.

此时不等式取等,平行.

综上,原不等式对应的结论永远成立.所以我们可以用叉积来判断两条直线的左右位置关系.

另外,叉积的几何意义是以这两个向量为邻边的平行四边形的有向面积(可用余弦定理证明),也就是两倍的三角形面积.

为什么说是有向面积呢?因为根据定义我们很容易发现cross(a,b)=cross(b,a),所以算出来会带正负号!

1 double dot(vt a,vt b){ return a.x*b.x+a.y*b.y; }
2 double cross(vt a,vt b){ return a.x*b.y-a.y*b.x; }
3 double area2(pt a,pt b,pt c){ return cross(b-a,c-a); }
View Code

 

2.向量的模

1 double length(vt a){ return sqrt(dot(a,a)); }
View Code

 

3.两向量的夹角

通过点积来计算夹角的余弦值,再用一个反三角函数.

1 double angle(vt a,vt b){ return acos(dot(a,b)/length(a)/length(b)); }
View Code

 

4.向量的旋转.

将一个向量绕起点逆时针旋转α弧度.

怎么操作呢?

我们先把这个向量的起点平移到坐标原点.写出向量终点的极坐标:

x=rcosθ

y=rsinθ

旋转之后终点的极坐标为:

x=rcos(θ+α)

y=rsin(θ+α)

打开得到:

x=r(cosθcosαsinθsinα)=xcosαysinα

y=r(sinθcosα+cosθsinα)=xsinα+ycosα

即:

x=xcosαysinα

y=xsinα+ycosα

1 vt rotate(vt a,double rad){ return vt(a.x*cos(rad)-a.y*sin(rad),a.x*sin(rad)+a.y*cos(rad)) }
View Code

 

5.复数

白书上提到了,但是我没有搞懂用来干嘛...貌似也比较慢,FFT里都是手写的,不管了.

 

3.点和直线


1.线段相交

由于这是出现最频繁的问题,所以我们放在第一个说.

首先来定义一下规范相交:两线段相交且交点不是四个端点中的任何一个.

在解决判断线段是否相交的问题时,最容易想到的是解析几何的办法.但是和斜率相关的点斜式和斜截式都用不了,如果用一般式或两点式需要联立解方程,然后判断是否有解,如果有解还要判断是否在两条线段内部.这样很麻烦,而且误差不好控制.

对于不要求求出交点坐标的问题,我们采用一种神奇的方法.

首先一个引理:两线段ABCD规范相交ACAD分别在AB左右两侧且CACB分别在CD左右两侧

所以我们用叉积做就好了.

1 bool segment_proper_intersection(pt a,pt b,pt c,pt d){
2     return(dcmp(cross(b-a,c-a))^dcmp(cross(b-a,d-a))==-2&&
3            dcmp(cross(d-c,a-c))^dcmp(cross(d-c,b-c))==-2);
4 }
View Code

 

注意这里的dcmp(cross(ba,ca)cross(ba,da))==2dcmp(cross(ba,ca),cross(ba,da))<0.

因为小于0的情况是一个是1,另一个是-1.在计算的运算中,全部用补码(正数不变,负数取反+1).

这样1的补码是0000...001

   -1的补码是1111...111

异或之后结果是1111...110

开头是1,是个负数,-1之后变为了1111...101,然后取反(从反码变回原码),1000...010,也就是-2.

 

那么对于不规范相交呢?我们发现不规范相交都会有四个端点中至少一个端点在另一个线段上.所以在用叉积进行判断的时候,如果两个向量的叉积为0,就说明有一个点在另一个线段所在的直线上,如果这个点在那个线段上,就是不规范相交啦.

所以我们还需要写一个函数来判断直线上的一个点是否在直线的某一条线段上.

我们可以用坐标比较的方法.还有第二种方法是用点积.叉积可以用来比较左右位置,而点积可以用来比较前后位置.

设这个点为P,线段为AB.如果PAPB的点积小于等于0,就说明P在线段AB上(当PAB重合时取等).

 1 int between(pt x,pt a,pt b){
 2     if(fabs(a.x-b.x)>fabs(a.x-b.y)) return dcmp(x.x-a.x)*dcmp(x.x-b.x);
 3     else return dcmp(x.y-a.y)*dcmp(x.y-b.y);
 4 }
 5 int between(pt x,pt a,pt b){
 6     return dcmp(rot(x-a,x-b));
 7 }
 8 int segment_intersection(pt a,pt b,pt c,pt d){
 9     d1=dcmp(cross(b-a,c-a));
10     d2=dcmp(cross(b-a,d-a));
11     d3=dcmp(cross(d-c,a-c));
12     d4=dcmp(cross(d-c,b-c));
13     if(d1==0&&between(c,a,b)<=0||
14        d2==0&&between(d,a,b)<=0||
15        d3==0&&between(a,c,d)<=0||
16        d4==0&&between(b,c,d)<=0)
17        return 2;
18     return 0;
19 }
View Code

 

 

2.求直线交点

这回是真的要求交点了.

首先我们用参数方程表示直线(一个点+一个向量)

1 struct line{
2     pt p; vt v;
3     line(){}
4     line(pt a,pt b){ p=a; v=b-a; }
5 };
View Code

直线上所有点Q都满足Q=P+tv.

我们设两条直线分别为P+tvQ+tw.设交点在这两条直线上的参数分别为t1,t2,作别为(x,y).

列出方程

x=xP+t1xv=xQ+t2xw

y=yP+t1yv=yQ+t2yw

我们设u=PQ,则xu=xPxQ,yu=yPyQ

由1式得到:

t2=xu+t1xvxw

带入2式消元,解得:

t1=cross(w,u)cross(v,w)

这样就可以得到交点:

P+v1t1

1 pt get_line_intersection(pt p,vt v,pt q,vt w){
2     vt u=p-q;
3     double t=cross(w,u)/cross(v,w);
4     return p+v*t;
5 }
View Code

 

注意到这里的cross(v,w)不能为0,也就是说v,w不能共线,也就是说在调用函数之前一定要确保两直线不平行.

 

3.求线段交点

我们可以用叉积的几何意义来求线段的交点.图嘛...看黑书P357吧...我这里大概推一下.

首先通过三角形等高容易证明:

|DP||CP|=SADBSACB=AD×ABAC×AB

然后我们用向量来证明黑书下面那个公式:

OP=OD+|DP||DP|+|CP|DC

我们设S1=SADB,S2=SACB

xP=xD+S1S1+S2(xCxD)

xP=S1xC+S2xDS1+S2

同理可以算出纵坐标.

 1 int segment_intersection(pt a,pt b,pt c,pt d,pt &p){
 2     double s1,s2;
 3     int d1,d2,d3,d4;
 4     d1=dcmp(s1=cross(b-a,c-a));
 5     d2=dcmp(s2=cross(b-a,d-a));
 6     d3=dcmp(cross(d-c,a-c));
 7     d4=dcmp(cross(d-c,b-c));
 8     if((d1^d2)==-2&&(d3^d4)==-2){
 9         p.x=(c.x*s2-d.x*s1)/(s2-s1);
10         p.y=(c.y*s2-d.y*s1)/(s2-s1);
11         return 1;
12     }
13     if(d1==0&&between(c,a,b)<=0||
14        d2==0&&between(d,a,b)<=0||
15        d3==0&&between(a,c,d)<=0||
16        d4==0&&between(b,c,d)<=0)
17        return 2;
18     return 0;
19 }
View Code

 

好了我写不下去了,明天再继续吧...今天晚上回去看凸包,半平面交,不看完是小狗...

转载于:https://www.cnblogs.com/Sunnie69/p/5577243.html

声明:本文内容由网友自发贡献,不代表【wpsshop博客】立场,版权归原作者所有,本站不承担相应法律责任。如您发现有侵权的内容,请联系我们。转载请注明出处:https://www.wpsshop.cn/w/知新_RL/article/detail/745561
推荐阅读
相关标签
  

闽ICP备14008679号