计算几何——多边形和圆相交的面积

 汇总篇：计算几何汇总 

---

 

原文地址：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);		
}