汇总篇:计算几何汇总
原文地址:https://www.cnblogs.com/lxglbk/archive/2012/08/12/2634192.html
代码部分为个人原创
圆与多变形相交的面积
因为任意简单多边形都可以划分成若干三角形,我们可以把这个简单多边形划分成三角形后,求三角形与圆的面积交,然后在把所有三角形的解合并。
计算一个圆与一个三角形的面积交(其中一个三角形顶点是圆心,如上图),我采用的方法是分4种情况。
1.另外2个顶点在圆内(上),这个非常好算直接求三角形的有向面积即可。
2.另外两个顶点有1个再圆内(上),另外1个再圆外,求得直线与圆一个交点后求一个三角形面积+上一个扇形面积。
需要函数:线段与圆的交点, 扇形面积
3.2个顶点在圆外,且2个顶点所在边与圆相交,先求圆外2顶点所在直线与圆交点,然后定比分点公式求另外2条直线与圆交点,然后求一个三角形+2个扇形面积即可。
4.2个顶点都在圆外且2顶点所在边与圆不相交,这个情况求2个交点后算出那个扇形面积就行了。
代码
const int eps = 1e-2;
double myRound(double a){//因为小数有误差,所以判断相切时要精确到固定位数
return floor(a * 100 + 0.5) / 100; /*保留小数点后2位*/
}
class point{
public:
double x;
double y;
point(double x_=0,double y_=0):x(x_),y(y_){}
friend const point operator+(const point& p1,const point& p2){
return point(p1.x+p2.x,p1.y+p2.y);
};
friend const point operator-(const point& p1,const point& p2){
return point(p1.x-p2.x,p1.y-p2.y);
};
friend const point operator*(const point& p,const double& m){
return point(p.x*m,p.y*m);
};
friend const point operator*(const double& m,const point& p){
return point(p.x*m,p.y*m);
};
friend const point operator/(const point& p,const double& m){
return point(p.x/m,p.y/m);
};
friend ostream& operator <<(ostream& out,point& a){
printf("(%lf,%lf)",a.x,a.y);
return out;
};
};
typedef point vect2;//重命名,向量也是用坐标表示
class line{
public:
point start;
point end;
line(point s=point(0,0),point e=point(0,0)):start(s),end(e){}
};
class circle{
public:
double r;
point O;
circle(point O_=point(0,0),double r_=0):O(O_),r(r_){}
};
class polygon{//这里使用浅拷贝,更快
public:
vector<point> p;
int n;
polygon(vector<point> q= vector<point>(),int n_=0){p=q;n=n_;}
};
//函数声明------------------------------------------------------------------
double dot(point O,point A,point B);
double cross(point O,point A,point B);
double dis(const point &p1,const point &p2);
double angle(point O,point A,point B);
double sectorArea(point O,point A,point B,double r);
double disOfPointToLine(point O,line l);
point shadowPointOfPointToLine(point O,line l);
bool pointIsOnLine(point O,line L);
bool interPointOfLineAndCircle(line l,point O,double r,vector<point> &ans);
double InterAreaOfCircleAndPolygon(polygon poly,circle c);
//--------------------------------------------------------------------------
double dot(point O,point A,point B){//点乘
double oa_x=A.x-O.x;
double oa_y=A.y-O.y;
double ob_x=B.x-O.x;
double ob_y=B.y-O.y;
return oa_x*ob_x+oa_y*ob_y;
}
double cross(point O,point A,point B){//叉乘
double oa_x=A.x-O.x;
double oa_y=A.y-O.y;
double ob_x=B.x-O.x;
double ob_y=B.y-O.y;
return oa_x*ob_y-oa_y*ob_x;
}
double dis(const point &p1,const point &p2){//求两点之间距离
double ans=(p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y);
return sqrt(ans);
}
double dis2(const point &p1,const point &p2){//求两点之间距离
return(p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y);
}
//两直线夹角
double angle(point O,point A,point B){
return acos(dot(O,A,B)/(dis(O,A)*dis(O,B)));
}
//扇形面积alpha*pi*r*r/2*pi
double sectorArea(point O,point A,point B,double r){
double alpha=angle(O,A,B);
return alpha*r*r/2;
}
//点到直线的距离
double disOfPointToLine(point O,line l){
double cos0=dot(l.start,O,l.end)/(dis(O,l.start)*dis(l.start,l.end));
return dis(O,l.start)*sin(acos(cos0));
}
//点在直线上的投影
point shadowPointOfPointToLine(point A,line l){//投影点出了问题
point B = l.start;
point C = l.end;
point D;
vect2 BC = C - B;
vect2 BA = B - A;
vect2 BD = dot(B,A,C)/dis2(B,C)*BC;
D = B + BD;
return D;
}
//点是否在线段上
bool pointIsOnLine(point O,line L){
point S=L.start;
point E=L.end;
if(cross(O,S,E)==0//在直线L上
and min(S.x,E.x)<=O.x&&O.x<=max(S.x,E.x)
and min(S.y,E.y)<=O.y&&O.y<=max(S.y,E.y))
return true;
else return false;
}
//线段到圆的交点,所需函数:点到直线的距离,点在直线上的投影,两点之间的距离,判断点是否在线段上
bool interPointOfLineAndCircle(line l,point O,double r,vector<point> &ans){
ans.clear();
point C,D,E;
double d= disOfPointToLine(O,l);
E = shadowPointOfPointToLine(O,l);
vect2 AB = l.end-l.start;
vect2 e=AB/dis(l.start,l.end);
double base2 =myRound(r*r-d*d);//取约数,不然本来是0的会得不到0
//相切base^2 =0 两个点相等,相离 base^2<0
if(base2<0)return false;
double base = sqrt(base2);
C = E - e*base;
D = E + e*base;
//判断C,D时候在线段AB上即可
if(pointIsOnLine(C,l))ans.push_back(C);
if(pointIsOnLine(D,l)&&!(D.x==C.x&&D.y==C.y))//避免点重复
ans.push_back(D);
if(ans.size()>0)return true;
else return false;
}
//多边形和圆相交的面积,所需函数:点到直线的距离,线段和圆的交点,扇形面积,两点之间的距离
double InterAreaOfCircleAndPolygon(polygon poly,circle c){
double sums=0;//有向面积
vector<point> pointsInCircle;
vector<point> pointsOutCircle;
vector<point> interPoints;
for(int i=0;i<poly.n;i++){
double s=0;
pointsInCircle.clear();
pointsOutCircle.clear();
line l(poly.p[i],poly.p[(i+1)%poly.n]);//取多边形的一条边
double d = disOfPointToLine(c.O,l);//圆心到这条边的距离
//判断边的两个端点是在圆内还是圆外
if(dis(l.start,c.O)<=d){
pointsInCircle.push_back(l.start);
}else pointsOutCircle.push_back(l.start);//靠l的start方向的点先入vector
if(dis(l.end,c.O)<=d){
pointsInCircle.push_back(l.end);
}else pointsOutCircle.push_back(l.end);
//分四种情况
if(pointsInCircle.size()==2){//1:两个点都在圆内
cout<<"三角形:"<<s<<endl;
s = cross(c.O,l.start,l.end)/2;
}else if(pointsInCircle.size()==1){//2:一个点在圆内,一个点在圆外
interPoints.clear();
//这里出错了没有求得交点
interPointOfLineAndCircle(l,c.O,c.r,interPoints);//只有一个交点
//定比分点求另一个交点,即圆外的顶点与圆心相连的线与圆的交点
vect2 e = (pointsOutCircle[0] - c.O)/dis(pointsOutCircle[0],c.O);
point interPoint2 = c.O + e*c.r;
double s1 = cross(c.O,pointsInCircle[0],interPoints[0]);
double s2 = sectorArea(c.O,interPoints[0],interPoint2,c.r);
cout<<"三角形:"<<s1<<"扇形:"<<s2<<endl;
s = s1 + s2;
}else if(pointsInCircle.size()==0){//两个顶点都在圆外
if(d<=c.r){//3:两顶点所在边与圆相交
interPoints.clear();
interPointOfLineAndCircle(l,c.O,c.r,interPoints);//这里有两个交点,在l的靠start方向的点先进入vector
//定比分点求另外两个交点,即圆外的顶点与圆心相连的线与圆的交点
vect2 e1 = (pointsOutCircle[0] - c.O)/dis(pointsOutCircle[0],c.O);
vect2 e2 = (pointsOutCircle[1] - c.O)/dis(pointsOutCircle[1],c.O);
interPoints.push_back(point(c.O + e1*c.r));
interPoints.push_back(point(c.O + e2*c.r));
double s1 = cross(c.O,interPoints[0],interPoints[1]);
double s2 = sectorArea(c.O,interPoints[0],interPoints[2],c.r);
double s3 = sectorArea(c.O,interPoints[1],interPoints[3],c.r);
printf("三角形:%lf 扇形:%lf %lf\n",s1,s2,s3);
s = s1 + s2 +s3;
}else{//4:两顶点所在边与圆不相交
interPoints.clear();
vect2 e1 = (pointsOutCircle[0] - c.O)/dis(pointsOutCircle[0],c.O);
vect2 e2 = (pointsOutCircle[1] - c.O)/dis(pointsOutCircle[1],c.O);
interPoints.push_back(point(c.O + e1*c.r));
interPoints.push_back(point(c.O + e2*c.r));
s = sectorArea(c.O,interPoints[0],interPoints[1],c.r);
cout<<"扇形:"<<s<<endl;
}
}
sums+=s;
cout<<"s="<<s<<endl;
}
return fabs(sums);
}