企业网站制作费用,淘宝评价采集wordpress,营销宝,桂林网站建设桂林首先基础知识#xff1a;
精度控制
计算几何经常牵扯到浮点数的运算#xff0c;所以就会产生精度误差#xff0c;因此我们需要设置一个eps#xff08;偏差值#xff09;#xff0c;一般取1e-7到1e-10之间#xff0c;并用下面的函数控制精度。#xff08;计算误差不得…首先基础知识
精度控制
计算几何经常牵扯到浮点数的运算所以就会产生精度误差因此我们需要设置一个eps偏差值一般取1e-7到1e-10之间并用下面的函数控制精度。计算误差不得小于这个设置的差值
const double eps 1e-8; //const 是c中用来定义常量的
int dcmp(double x){if(fabs(x)eps) return 0; //fabs是求浮点小数的绝对值 else return x0?-1:1;
}
向量
有大小有方向的量又称为矢量。
二维的向量常用一个对数(x,y)表示代码中常用一个结构体来实现向量
struct vector{double x,y;vector (double X 0,double Y 0){x X;y Y;}};
向量的模,即向量的长度
a(x,y)
|a|sqrt(x^2y^2)
double len(vector a){return sqrt(a.x*a.xa.y*a.y)
}
二维平面中的点同样可以用对数(x,y)来表示所以向量的存储方式同样可以用于点
typedef vector point;
需要注意的是1点加减向量为点2点减点为向量
极角
极坐标一共两个参数极角和极径 一个函数图像上某一点到原点的距离就是极径,极径与x轴的夹角就是极角.
对于向量a(x,y)可以用函数atan2(y,x)来计算他的极角
按照极角为关键字排序后的顺序为极角序
向量的四则运算
vector operator (vector a,vector b) {return vector (a.xb.x ,a.yb.y);
}
vector operator - (vector a,vector b){return vector (a.x-b.x ,a.y - b.y);
}
vector operator * (vector a,double p){return vector (a.x*p,a.y*p);
}
vector operator / (vector a,double p){return vector (a.x/p,a.y/p);
}点积
a·b的几何意义为a在b上的投影长度乘以b的模长
a·b|a||b|cosθ,其中θ为a,b之间的夹角
坐标表示
a(x1,y1) b(x2,y2)
a·bx1*x2y1*y2;
double dot(vector a,vector n){return a.x*b.xa.y*b.y;
}
点积的应用
1判断两个向量是否垂直 a⊥b a·b0
2求两个向量的夹角点积0为钝角点积0为锐角
double Angle(vector a,vector b){return acos(dot(a,b)/len(a)/len(b));
}
3求模长
double len(vector a)
{ return sqrt(dot(a,a));
}法向量与单位向量垂直的向量称为单位法向量
vector normal(vector a){double l len(a);return vector (-a.y/l,a.x/l);
}
二维叉积
两个向量的叉积是一个标量a×b的几何意义为他们所形成的平行四边形的有向面积
坐标表示a(x1,y1) b(x2,y2)
a×bx1y2-x2y1
double cross(vector a,vector b){return a.x*b.y - a.y*b.x;
}
直观理解假如b在a的左边则有向面积为正假如在右边则为负。假如b,a共线则叉积为0,。所以叉积可以用来判断平行。
向量的旋转
a(x,y)可以看成是x*(1,0)y1*(0,1)
分别旋转两个单位向量则变成x*(cosθsinθ)y1*(-sinθcosθ)
点旋转
x,y绕原点逆时针旋转 得 a(x cos a - y sin a , x sin a y cos a )
点、直线、线段的关系
1 、
点到直线的距离
利用叉积求面积然后除以平行四边形的底边长得到平行四边形的高即点到直线的距离
double distl(point p,point a,point b){vector v p-a;vector u b-a;return fabs(cross(u,v)/len(u));
}
2、
点到线段的距离
比点到直线的距离稍微复杂。因为是线段所以如果平行四边形的高在区域之外的话就不合理这时候需要计算点到距离较近的端点的距离
double dists(point p,point a, point b){if (ab)return len(p-a);vector v1 b-a, v2 p -a ,v3 p -b;if(dcmp(dot(v1,v2))0)return len(v2);else if(dcmp(dot(v1,v3))0)return len(v3);return fabs(corss(v1,v2))/len(v1);
}
3、
判断两个线段是否相交
跨立实验判断一条线段的两端是否在另一条线段的两侧两个端点与另一线段的叉积乘积为负。需要正反判断两侧。
bool segment(point a,point b,point c,point d){double c1 cross ( b-a,c-a),c2 cross(b-a,d-a);double d1 cross (d-c,a-c),d2 cross(d-c,b-c);return dcmp(c1)*dcmp(c2)0dcmp(d1)*dcmp(d2)0;
}
4、
求两条直线的交点
有两种方法比较常用的一种是用叉积的比值计算。但是这种方法的精度不是很高。
point glt(point a,point a1,point b,point b1)
{ vector va1-a; vector wb1-b; vector ua-b; double tcross(w,u)/cross(v,w); return av*t;
} 还有一种比较麻烦不常用但是精度相对较好
point line_intersection(point a,point a0,point b,point b0)
{ double a1,b1,c1,a2,b2,c2; a1 a.y - a0.y; b1 a0.x - a.x; c1 cross(a,a0); a2 b.y - b0.y; b2 b0.x - b.x; c2 cross(b,b0); double d a1 * b2 - a2 * b1; return point((b1 * c2 - b2 * c1) / d,(c1 * a2 - c2 * a1) / d);
} 多边形相关 判断点在多边形内部
射线法以该点为起点引一条射线与多边形的边界相交奇数次说明在多边形的内部。
int pointin(point p,point* a,int n)
{ int wn0,k,d1,d2; for (int i1;in;i) { if (dcmp(dists(p,a[i],a[(i1-1)%n1]))0) return -1;//判断点是否在多边形的边界上 kdcmp(cross(a[(i1-1)%n1]-a[i],p-a[i])); d1dcmp(a[i].y-p.y); d2dcmp(a[(i1-1)%n1].y-p.y); if (k0d10d20) wn; if (k0d20d10) wn--; } if (wn) return 1; else return 0;
} 求多边形的面积
再多边形内取一个点进行三角剖分用叉积求三角形的面积。因为叉积是有向面积所以任意多边形都使用。
注意最后取绝对值。
double PolygonArea(Point* p,int n)
{double area0;for(int i1;in-1;i)areaCross(p[i]-p[0],p[i1]-p[0]);return area/2;
}求多边形的重心 同样方法将多边形三角剖分
算出每个三角形的重心
套用质点组的重心公式即可
质点组重心公式 三个点A,B,C
x(ma*xamb*xbmc*xc)/(mambmc)
y(ma*yamb*ybmc*yc)/(mambmc) m表示权三角形的有向面积
没写过用到再说啦
凸包
凸包就是把给定点包围在内部面积最小的凸多边形
Graham扫描法
1.把点按照x坐标排序
2.将p1,p2放到凸包中
3.从p3开始当新的点在凸包的前进方向的左边时加入该点否则弹出最后加入的点直到新点在左边
4.将上述流程反过来做一遍求上凸壳
void convexhull(point *p,int n,point *ch)
{ sort(p1,pn1); if (n1){ ch[m]n; return ; } m1; for (int i1;in;i){ while (m2dcmp(cross(ch[m-1]-ch[m-2],p[i]-ch[m-2]))0) m--; ch[m]p[i]; } int km; for (int in-1;i1;i--){ while (mkdcmp(cross(ch[m-1]-ch[m-2],p[i]-ch[m-2]))0) m--; ch[m]p[i]; } m--;
} 半平面交
增量法初始时加入一个范围巨大的框
每次拿新的半平面切割原来的凸集保留在新加直线左边的点删除右边的有向直线与多边形相交产生的新点加入到新的多边形中。
时间复杂度O(n^2)
void cut(point a,point b)
{ point tmp[N]; int cnt0; for (int i0;im;i){ double ccross(b-a,p[i]-a); double dcross(b-a,p[(i1)%m]-a); if (c0) tmp[cnt]p[i]; if (c*d0) tmp[cnt]glt(a,b,p[i],p[(i1)%m]); } mcnt; for (int i0;icnt;i) p[i]tmp[i];
} 旋转卡壳
求凸包上最远点对的距离
被凸包上一对平行直线卡住的点对叫做对踵点答案一定在对踵点上。 按照逆时针顺序枚举凸包上一条边则凸包上到该边所在直线最远的点也单调逆时针旋转相当于用一条直线卡对面一个顶点
double rotating_calipers(point* ch,int n)//旋转卡壳被凸包上被一对平行直线卡住的点叫对踵点最远点对一定在凸包的一对对踵点上
{ if (n1) return 0; if (n2) return dot(ch[0]-ch[1],ch[0]-ch[1]);//如果只有两个点那么就是两点的直线距离 int now1,i; double ans0; ch[n]ch[0]; for (int i0;in;i) { while(dcmp(distl(ch[now],ch[i],ch[i1])-distl(ch[now1],ch[i],ch[i1]))0)//最远点随着平行线的旋转是单调的所以点不会来回移动 now(now1)%n; ansmax(ans,dot(ch[now]-ch[i],ch[now]-ch[i])); ansmax(ans,dot(ch[now]-ch[i1],ch[now]-ch[i1]));//找到与当前直线平行的最远点用来更新答案 } return ans;
} 求两凸包上的最近距离
double rotating_calipers(point* p,point* q,int n,int m)//利用旋转卡壳求两个凸包上最近点的距离一个凸包上的平行线逆时针旋转另一个凸包上的最远点也单调逆时针旋转所以这个算法要正反进行两遍
{ int x0; int y0,i; double ans1e10,t; for (int i0;in;i) x(p[i].yp[x].y)?i:x; for (int i0;im;i) y(p[i].yp[y].y)?i:y; p[n]p[0]; q[m]q[0]; for (int i0;in;i) { while((tdcmp(cross(p[x]-p[x1],q[y1]-q[y])))0)//判断凸壳上下一条边的旋转方向 y(y1)%m; if (t0)//平行的时候四个点之间更新答案 { ansmin(ans,dists(p[x],q[y1],q[y])); ansmin(ans,dists(p[x1],q[y1],q[y])); ansmin(ans,dists(q[y],p[x],p[x1])); ansmin(ans,dists(q[y1],p[x],p[x1])); } else ansmin(ans,dists(q[y],p[x],p[x1]));//否则只用最远点更新答案 x(x1)%n; } return ans;
} 其实旋转卡壳的应用有很多变形因为只要是凸包上满足单调性的一些问题大多都可以用旋转卡壳来解。
最小圆覆盖
用最小的圆覆盖平面中的所有点。
point center(point a,point b,point c){ point p(ab)/2; point q(ac)/2; vector vrotate(b-a,pi/2); vector urotate(c-a,pi/2); if (dcmp(cross(v,u))0) { if(dcmp(len(a-b)len(b-c)-len(a-c))0) return (ac)/2; if(dcmp(len(a-c)len(b-c)-len(a-b))0) return (ab)/2; if(dcmp(len(a-b)len(a-c)-len(b-c))0) return (bc)/2; } return glt(p,v,q,u);
}
double mincc(point *p,int n,point c)
{ random_shuffle(p,pn); cp[0]; double r0; int i,j,k; for (i1;in;i) if (dcmp(len(c-p[i])-r)0){ cp[i],r0; for (j0;ji;j) if (dcmp(len(c-p[j])-r)0){ c(p[i]p[j])/2; rlen(c-p[i]); for (k0;kj;k) if (dcmp(len(c-p[k])-r)0){ ccenter(p[i],p[j],p[k]); rlen(c-p[i]); } } } return r;
}
