【计算几何各种小模板总结贴】[不定期更新]

相当全的计算几何模板
http://paste.ubuntu.com/24458430/

精度控制

Ps：尽量不要用除法，三角函数，强制类型转换等操作，否则精度损失比较大
const double PI = acos(-1.0);
const double EPS = 1e-8;
任何double的比较运算都要用eps

向量运算

二维空间

向量旋转矩阵

我们想将向量 x,y−→−\overrightarrow {x, y} 以xx为轴点逆时针旋转，且旋转角为 α\alpha ，
[ x′y′]=[ cos(α)sin(α)−sin(α)cos(α)]×[ xy] \left[\begin{array}{cc}\ x' \\ y' \end{array}\right]=\left[\begin{array}{cc}\ \cos(\alpha) &&-\sin(\alpha) \\ \sin(\alpha) && \cos(\alpha) \end{array}\right] \times \left[\begin{array}{cc}\ x \\ y \end{array}\right]

叉积

double mul(Point p1,Point p2,Point p0){return((p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y));
}

两点求距离

double dis(Point p1,Point p2){return(sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y)));
}

三维空间

点,线,面

struct point{double x,y,z;
};
struct line{point a,b;
};
struct plane{point a,b,c;
};

计算 cross product product product product U x V

point xmult(point u,point v){point ret;ret.x=u.y*v.z-v.y*u.z;ret.y=u.z*v.x-u.x*v.z;ret.z=u.x*v.y-u.y*v.x;return ret;
}

计算 dot product product product product U . V

double dmult(point u,point v){return u.x*v.x+u.y*v.y+u.z*v.z;
}

矢量差 U - V

point subt(point u,point v){point ret;ret.x=u.x-v.x;ret.y=u.y-v.y;ret.z=u.z-v.z;return ret;
}

取平面法向量

point pvec(plane s){return xmult(subt(s.a,s.b),subt(s.b,s.c));
}
point pvec(point s1,point s2,point s3){return xmult(subt(s1,s2),subt(s2,s3));
}

两点距离,单参数取向量大小

double dis(point p1,point p2){return sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y)+(p1.z-p2.z)*(p1.z-p2.z));
}

向量大小

double vlen(point p){return sqrt(p.x*p.x+p.y*p.y+p.z*p.z);
}

判断四点是不是共面

bool judge(point a,point b,point c,point d){double tmp = dmult(prec(a,b,c),smult(d,a));return  ( abs(tmp) < eps );
}

点到平面距离

double ptoplane(point p,point s1,point s2,point s3){return fabs(dmult(pvec(s1,s2,s3),subt(p,s1)))/vlen(pvec(s1,s2,s3));
}

三角形面积

double Area_triangle(point a,point b,point c){double ab=dis(a,b),bc=dis(b,c),ac=dis(a,c);double p=(ab+bc+ac)/2;return sqrt(p* (p-ab) * (p-bc)*(p-ac));
}
double Area_triangle(point b,point p1,point p2){point a=xmult(smult(b,p1),smult(b,p2));return sqrt(a.x*a.x+a.y*a.y+a.z*a.z)/2.0;
}

极角排序

//极角排序

bool dy(double x,double y)  {   return x > y + eps;} // x > y
bool xy(double x,double y)  {   return x < y - eps;} // x < y
bool dyd(double x,double y) {   return x > y - eps;} // x >= y
bool xyd(double x,double y) {   return x < y + eps;} // x <= y
bool dd(double x,double y)  {   return fabs( x - y ) < eps;}  // x == ybool cmp(point a,point b)   // 第一次排序
{if( dd(a.y ,b.y ) )  return xy(a.x, b.x);return xy(a.y,b.y);
}

凸包

grabam 扫描法

O(nlog2n)O(nlog_{2}{n})复杂度寻求凸包

凸包运算（周长+面积+最小园覆盖）

const int maxN=50000;
Point p[maxN];
Point ch[maxN];
int n;
int len;int main()
{int t,l;//scanf("%d",&t);while(~scanf("%d",&n)){double distence=0.0;for(int i=0; i<n; i++){scanf("%lf%lf",&p[i].x,&p[i].y);}if(n==0) break;//这是一些值得特判的情况//if(n==0) if(n==1) if(n==2)if(n<=1){printf("0.50\n");continue;}if (n == 2){printf("%.2lf\n",dis(p[0],p[1])*0.50+0.50);continue;}Graham_scan(p,ch,n,len);for(int i=0; i<len; i++) //凸包 各个顶点的顺序cout<<ch[i].x<<" "<<ch[i].y<<endl;//求凸包的周长double perimeters=0.0;for(int i=1;i<len;i++){perimeters+=dis(ch[i],ch[i-1]);}perimeters+=dis(ch[0],ch[len-1]);printf("%.0lf\n",perimeters);//求凸包的面积double area=0.0;for(int i=2; i<len; i++){area+=mul(ch[i-1],ch[i],ch[0]);}printf("%d\n",(int )area/100);//求图包最小圆覆盖double maxr = -1;double a, b, c, r, s;for (int i=0; i<len; i++)  //枚举凸包上的点{for (int j=i+1; j<len; j++){for (int k=j+1; k<len; k++){a = dis(ch[i], ch[j]);b = dis(ch[i], ch[k]);c = dis(ch[j], ch[k]);if (a*a+b*b<c*c || a*a+c*c<b*b || b*b+c*c<a*a)r = max(max(a, b), c) / 2.0;//钝角时  半径为最长边的一半else{s = fabs(mul(ch[i], ch[j], ch[k])) / 2.0;r = a * b * c / (s * 4.0);//三角形外接圆公式}if (maxr < r) maxr = r;}}}printf ("%.2lf\n", maxr+0.50);//if(t) printf("\n");  //这是控制两组输出数据之间有一个空格的}return 0;
}

线与线

线段相交

bool IsIntersected(point s1,point e1,point s2,point e2)//两个线段相交
{return(max(s1.x,e1.x)>=min(s2.x,e2.x))&&(max(s2.x,e2.x)>=min(s1.x,e1.x))&&(max(s1.y,e1.y)>=min(s2.y,e2.y))&&(max(s2.y,e2.y)>=min(s1.y,e1.y))&&(multi(s1,s2,e1)*multi(s1,e1,e2)>0)&&(multi(s2,s1,e2)*multi(s2,e2,e1)>0);
}

求两线段交点

point intersection(point &A,point &B,point &C,point &D)
{/* AB与CD交点*/point p;double a1=A.y-B.y;double b1=B.x-A.x;double c1=A.x*B.y-B.x*A.y;double a2=C.y-D.y;double b2=D.x-C.x;double c2=C.x*D.y-D.x*C.y;p.x=(b1*c2-b2*c1)/(a1*b2-a2*b1);p.y=(a2*c1-a1*c2)/(a1*b2-a2*b1);return p;
}

两圆相交

两圆相交面积

struct Round {  double x, y;  double r;
};
double solve(Round a, Round b)
{double d = dis(a, b);if (d >= a.r + b.r)return 0;if (d <= fabs(a.r - b.r)){double r = a.r < b.r ? a.r : b.r;return PI * r * r;}double ang1 = acos((a.r * a.r + d * d - b.r * b.r) / 2. / a.r / d);double ang2 = acos((b.r * b.r + d * d - a.r * a.r) / 2. / b.r / d);double ret = ang1 * a.r * a.r + ang2 * b.r * b.r - d * a.r * sin(ang1);return ret;
}

四点共面

可以由4个点构成3个向量，3个向量共面的充要条件是向量为a，b，c，存在实数x,y,z不全为零，使得xa+yb+zc=0。转化为线性代数的3个向量线性相关的行列式为0。

double is_coplanar(point a, point b, point c, point d){point ab, ac, ad;ab.x=a.x-b.x, ab.y=a.y-b.y, ab.z=a.z-b.z;ac.x=a.x-c.x, ac.y=a.y-c.y, ac.z=a.z-c.z;ad.x=a.x-d.x, ad.y=a.y-d.y, ad.z=a.z-d.z;double gg=ab.x*ac.y*ad.z+ab.y*ac.z*ad.x+ac.x*ad.y*ab.z;gg-=ab.z*ac.y*ad.x+ab.x*ac.z*ad.y+ab.y*ac.x*ad.z;return gg;
}