计算几何：常见几何元素关系的判断以及计算

求三角形面积，海伦公式：S=sqrt(p*(p-a)*(p-b)*(p-c))；S=(b*c*sina)/2

 

在几何中，double二元组可以表示顶点，表示向量，表示线段，表示直线。一般直接把二元组看成向量来计算。向量最重要的两个计算，点积和叉积具有不同几何意义。

点积，u*v=x1*x2+y1*y2。通常用来计算两向量之间的夹角（角度=acos(u*v)），以及一个向量在另一个向量上的投影。叉积，u^v=x1*y2-x2*y1，结果是一个向量，遵守右手定则，结果的正负就代表了方向。叉积结果的绝对值，是平行四边形的面积；I=AB^BC时，I>0代表C在AB的左侧，I<0代表C在AB的右侧，I=0代表C在AB上。

需要注意的是要注意double在判断时的误差。由于double存储方式的原因，double对于不能精确表示的数只能用所能表示的数中最接近的数表示。在计算几何中，尤其要注意double的误差。选择合适的eps来限制精度

const double eps=1e-10;
a<0 --> a<-eps
a<=0 --> a<eps
a==0 --> fabs(a)<eps

接下来是一些常见的集合元素间位置关系判断的模块：

1.求直线交点，应当先对两条直线进行判断，若相交则求交点

P a,b,c,d;   //两条直线上的两点，P为自定义类型，代表向量（内含两个double类型的数，表示横纵坐标）
if((a-b)^(c-d)==0){
  //平行
}
else{
  //相交
  double k=((c-a)^(d-c))/(b-a)^(d-c));
  return a+(b-a)*k;
}

2.求线段交点。分为两种情况，若线段所在直线平行，则判断a线段的两个端点有没有可能在b线段上，以及b线段的两个端点有没有可能在a线段上；若线段所在直线相交，先求出直线交点，再判断交点在不在两条线段上。

P lineinter(P a,P b,P c,P d){
  ... //直线球交点
}
bool onseg(P a,P b,P r){  //判断r是否在以a、b为端点的线段上
  return (a-r)^(b-r)==0 && (a-r)*(b-r)<=0;
}
P seginter(P a,P b,P c,P d){
  if((a-b)^(d-c)==0){
    if(onseg(a,b,c)) ...
    else if(onseg(a,b,d)) ...
    else if(onseg(c,d,a)) ...
    else if(onseg(c,d,b)) ...
  }
  else{
    P r=lineinter(a,b,c,d);
    if(onseg(a,b,r)&&onseg(c,d,r)) return r;
  }
}

3.点到线段的最短距离

点到线段的距离要么是点到直线的距离，要么是点到两个端点距离的较小值。

double dis(P a,P b,P c){
  if(((c-a)*(b-a))^((c-b)*(a-b))>=0)  //判断点在直线上的投影是不是在线段上
    return fabs((c-a)^(c-b)/dis(a,b));
  return min(dis(a,c),dis(b,c));
}

4.直线与圆求交

三种关系：无交点、一个交点、两个交点。用圆心到直线的距离和半径进行判断

int circlelineinter(P a,P b,circle C,P &p1,P &p2){
  point o=C.c; double r=C.r;
  double d=fabs((a-o)^(b-o))/dis(a,b);  //圆心到直线的距离
  if(d>r) return 0;
  else if(fabs(d-r)<eps){
    double k=(o-a)*(b-a)/dis(a,b);
    p1=a+(b-a)*k/dis(a,b);
    return 1;
  }
  else{
    double k1=(o-a)*(b-a)/dis(a,b);
    p1=a+(b-a)*k1/dis(a,b);
    double k2=(o-b)*(a-b)/dis(a,b);
    p2=b+(a-b)*k2/dis(a,b);
    return 2;
  }
}

5.圆与圆相交

五种情况，内含，内切，相交，外切，相离。写程序时，认为内切外切是相交的极限情况，两圆之间要么没有交点，要么两个交点。首先是关于两圆之间关系的判断同时求交点

int circleinter(circle a,circle b,point &p1,point &p2){
  if(a.r>b.r) swap(a,b);
  double r=a.r,R=b.r,d=(a.c-b.c).len();
  //两圆重合
  if(fabs(d)<eps&&fabs(a.c.x-b.c.x)<eps&&fabs(a.c.y-b.c.y)<eps) return 0;
  if(R>r+d || R+r<d) return 0;   //内含或相离
  double cosang=(r*r+d*d-R*R)/(2*r*d);
  double sinang=sqrt(1.0-cosang*cosang);
  point o=(b.c-a.c)*(1.0/d);   //a.c->b.c方向上的单位向量
  o=a.c+o*(r*cosang);     //向量所指向的点是两圆相交线段的中点
  point p=point(a.c.y-b.c.y,b.c.x-a.c.x)*(1.0/d);
  p1=o+pt*(r*sinang);
  p2=o-pt*(r*sinang);
  if(fabs(p1.x-p2.x)<eps&&fabs(p1.y-p2.y)<eps) return 1;  //两交点重合
  return 2;
}

两圆求相交面积，用积分来做

double circleinter(circle a,circle b){
  if(a.r>b.r) swap(a,b);
  double r=a.r,R=b.r,d=(a.c-b.c).len();
  if(d<=R-r) return acos(-1.0)*r*r;
  if(r+R<d+eps) return 0;
  double x1=(r*r+d*d-R*R)/(2*d),x2=(R*R+d*d-r*r)/(2*d);
  double ans=x1*sqrt(r*r-x1*x1)+r*r*acos(x1/r);
  ans+=x2*sqrt(R*R-x2*x2)+R*R*acos(x2/r);
  return ans;
}

此外，务必要考虑到极限情况，这也是计算几何繁琐的原因，常见的极限情况：

1.直线与多边形相交，需要判断直线与点的相交情况

2.两直线平行、三点共线

3.判断实心物体是否相交，注意一个在另一个内部的情况

4.模拟物体运动，将初始的相互接触的状态认为是碰撞