LightOJ 1130 Intersection between Circle and Rectangle

最新推荐文章于 2020-04-05 19:05:54 发布

原创最新推荐文章于 2020-04-05 19:05:54 发布 · 1.3k 阅读

0 ·

CC 4.0 BY-SA版权

文章标签：

#LightOJ #精度 #面积

计算几何同时被 2 个专栏收录

56 篇文章

订阅专栏

LightOJ

54 篇文章

订阅专栏

本文介绍了一种计算给定圆与矩形重叠部分面积的方法，通过构造图形并将其分解为多个三角形和扇形来计算面积。特别关注了计算精度问题，并提供了实现代码。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

题目链接 http://lightoj.com/volume_showproblem.php?problem=1130

题意：求所给出的圆与矩形的重叠部分的面积

思路：构造出所求图形，求出矩形与圆的交点，将阴影部分分成若干三角形和扇形，分块计算面积。题目对计算精度要求很高

这个题因为精度问题纠结了好久……最后参考了一下别人的思路。

参考博文 http://blog.youkuaiyun.com/pvpishard/article/details/7220876

关于参考博文中提到的奇葩问题，我测试的时候全是WA……应该是LightOJ的测试数据更新了吧。

关于参考博文中提到“最精简的求1/sqrt()”的函数貌似精度只有小数点后3位，应该是不满足本题要求（小数点后6位）。

这个函数我也试着了解了一下，在较低精度要求下真的快不少！

参考链接 http://game.chinaitlab.com/devdoc/723647.html

http://www.cnblogs.com/nihgwu/archive/2006/08/04/467597.html

这个题应该还有更好的思路，LightOJ上目前最小MEMORY占用居然只有228……欢迎留言，一起讨论

#include <stdio.h>
#include <math.h>

const double PI=acos(-1.0);
const double STD=1e-9;      //浮点数大小判断标准

struct Point
{
	double x,y;
	void disp ()
	{
		printf("x=%lf y=%lf\n",x,y);
	}
	Point friend Vec (Point a,Point b)   //构造向量
	{
		Point res;
		res.x=a.x-b.x;
		res.y=a.y-b.y;
		return res;
	}
}pt[20];

struct Line
{
	Point s,e;
	void set (double sx,double sy,double ex,double ey)
	{
		s.x=sx;s.y=sy;e.x=ex;e.y=ey;
	}
}line[4];

struct Circle
{
	Point c;
	double r;
}cir;

int point_num;       //记录需要计算的点的数量
bool out[20];        //标记点与圆的位置关系,在圆内和圆上为false,在圆外为true

void Find (Line l,int mode)  //寻找矩形与圆的交点,结果保存在pt数组中
{
	Point temp;
	double t,d;
	switch (mode)
	{
	case 0:                 //寻找"左至右"边与圆的交点
		if (fabs(cir.c.y-l.s.y) > cir.r)   //与圆不相交
			return;
		d=cir.c.y-l.s.y;
		t=sqrt(fabs(cir.r*cir.r-d*d));
		temp.y=l.s.y;
		temp.x=cir.c.x-t;
		if (temp.x>l.s.x && temp.x<l.e.x)
			pt[point_num++]=temp;
		temp.x=cir.c.x+t;
		if (temp.x>l.s.x && temp.x<l.e.x)
			pt[point_num++]=temp;
		break;
	case 1:                 //寻找"上至下"边与圆的交点
		if (fabs(cir.c.x-l.s.x) > cir.r)   //与圆不相交
			return;
		d=cir.c.x-l.s.x;
		t=sqrt(fabs(cir.r*cir.r-d*d));
		temp.x=l.s.x;
		temp.y=cir.c.y+t;
		if (temp.y<l.s.y && temp.y>l.e.y)
			pt[point_num++]=temp;
		temp.y=cir.c.y-t;
		if (temp.y<l.s.y && temp.y>l.e.y)
			pt[point_num++]=temp;
		break;
	case 2:                   //寻找"右至左"边与圆的交点
		if (fabs(cir.c.y-l.s.y) > cir.r)     //与圆不相交
			return;
		d=cir.c.y-l.s.y;
		t=sqrt(fabs(cir.r*cir.r-d*d));
		temp.y=l.s.y;
		temp.x=cir.c.x+t;
		if (temp.x<l.s.x && temp.x>l.e.x)
			pt[point_num++]=temp;
		temp.x=cir.c.x-t;
		if (temp.x<l.s.x&&temp.x>l.e.x)
			pt[point_num++]=temp;
		break;
	case 3:                //寻找"下至上"边与圆的交点
		if (fabs(cir.c.x-l.s.x) > cir.r)    //与圆不相交
			return;
		d=cir.c.x-l.s.x;
		t=sqrt(fabs(cir.r*cir.r-d*d));
		temp.x=l.s.x;
		if (l.s.y<l.e.y)    //寻找"下至上"边与圆的交点
		temp.y=cir.c.y-t;
		if (temp.y>l.s.y && temp.y<l.e.y)
			pt[point_num++]=temp;
		temp.y=cir.c.y+t;
		if (temp.y>l.s.y && temp.y<l.e.y)
			pt[point_num++]=temp;
	}	
}
/*            
void Find (Line l)   //网上别人代码的修改+注释版
{
	Point temp;
	double t,d;
	if (fabs(l.s.x-l.e.x) < STD) //起始横坐标相同,即竖直线
	{
		if (fabs(cir.c.x-l.s.x) > cir.r) //起点在圆外
			return;
		d=cir.c.x-l.s.x;
		t=sqrt(fabs(cir.r*cir.r-d*d));
		temp.x=l.s.x;
		if (l.s.y<l.e.y) //寻找"下至上"边与圆的交点
		{
			temp.y=cir.c.y-t;
			if (temp.y>l.s.y && temp.y<l.e.y)
				pt[point_num++]=temp;
			temp.y=cir.c.y+t;
			if (temp.y>l.s.y && temp.y<l.e.y)
				pt[point_num++]=temp;
		}
		else  //寻找"上至下"边与圆的交点
		{
			temp.y=cir.c.y+t;
			if (temp.y<l.s.y && temp.y>l.e.y)
				pt[point_num++]=temp;
			temp.y=cir.c.y-t;
			if (temp.y<l.s.y && temp.y>l.e.y)
				pt[point_num++]=temp;
		}
	}
	else    //起始横坐标不同,即水平线
	{
		if (fabs(cir.c.y-l.s.y) > cir.r)     //起点在圆外
			return;
		d=cir.c.y-l.s.y;
		t=sqrt(fabs(cir.r*cir.r-d*d));
		temp.y=l.s.y;
		if (l.s.x<l.e.x)          //寻找"左至右"边与圆的交点
		{
			temp.x=cir.c.x-t;
			if (temp.x>l.s.x && temp.x<l.e.x)
				pt[point_num++]=temp;
			temp.x=cir.c.x+t;
			if (temp.x>l.s.x && temp.x<l.e.x)
				pt[point_num++]=temp;
		}
		else       //寻找"右至左"边与圆的交点
		{
			temp.x=cir.c.x+t;
			if (temp.x<l.s.x && temp.x>l.e.x)
				pt[point_num++]=temp;
			temp.x=cir.c.x-t;
			if (temp.x<l.s.x&&temp.x>l.e.x)
				pt[point_num++]=temp;
		}
	}
}*/

void Input ()         //输入及构造部分
{
	scanf("%lf%lf%lf",&cir.c.x,&cir.c.y,&cir.r);
	double lx,ly,rx,ry;
	scanf("%lf%lf%lf%lf",&lx,&ly,&rx,&ry);
	
	line[0].set(lx,ry,rx,ry);        //上边-左至右
	line[1].set(rx,ry,rx,ly);        //右边-上至下
	line[2].set(rx,ly,lx,ly);        //下边-右至左
	line[3].set(lx,ly,lx,ry);        //左边-下至上

	//求出需要考虑的点,即矩形顶点以及矩形与圆的交点
	point_num=0;
	pt[point_num++]=line[0].s;
	Find (line[0],0);
	pt[point_num++]=line[1].s;
	Find (line[1],1);
	pt[point_num++]=line[2].s;
	Find (line[2],2);
	pt[point_num++]=line[3].s;
	Find (line[3],3);
	
//	for (int i=0;i<point_num;i++) printf("Point %d ",i),pt[i].disp();
}

double get_dis (Point a,Point b)        //计算两点间距离
{
	return sqrt((a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y));
}

double get_dis2 (Point a,Point b)       //计算两点间距离的平方
{
	return (a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y);
}

void Judge ()
{
	double t;
	for (int i=0;i<point_num;i++)
	{
		t=get_dis(pt[i],cir.c);
		if (t>cir.r+STD)
			out[i]=true;
		else out[i]=false;
	}
}

double get_cj (Point a,Point b)        //计算叉积
{
	return a.x*b.y-a.y*b.x;
}

double get_dj(Point a,Point b)        //计算点积
{
	return a.x*b.x+a.y*b.y;
}


double get_tri (Point a,Point b)     //计算三角形面积
{
	return get_cj (Vec(a,cir.c),Vec(b,cir.c))/2;
}

double get_sec (Point a,Point b)     //计算扇形面积
{
	double da=get_dis2 (a,cir.c);
	double db=get_dis2 (b,cir.c);
	double cj=get_cj (Vec(a,cir.c),Vec(b,cir.c));
	double dj=get_dj (Vec(a,cir.c),Vec(b,cir.c));
	double dsin=cj/sqrt(da*db);       //da,db是平方过的,如果先开方再乘会损失精度WA……
	double alpha=asin(dsin);

	//asin函数求出的角是0到90度的,需要根据实际情况还原到原来的大小
	if (dj < -STD)           //当夹角大于90度时
		if (alpha>0)
			alpha=PI-alpha;
		else
			alpha=-PI-alpha;

	return alpha*cir.r*cir.r/2;
}

double Deal ()
{
	if (cir.r < STD)
		return 0;
	Judge ();

	double ans=0;
	for (int i=0;i<point_num;i++)
	{
		if (out[i] == false && out[(i+1)%point_num] == false)
			//如果两个点都不在圆外
			ans+=get_tri(pt[i],pt[(i+1)%point_num]);
		else  //若至少有一点在圆外
			ans+=get_sec(pt[i],pt[(i+1)%point_num]);
	}
	return fabs(ans);       //用叉积可能负面积
}

int main ()
{
	int T;
	scanf("%d",&T);
	for (int cas=1;cas<=T;cas++)
	{
		Input ();
		printf("Case %d: %.8lf\n",cas,Deal());
	}
	return 0;
}

/*
2
4 4 4
1 1 3 3
4 4 4
1 1 7 7

Out
Case 1: 3.93987658
Case 2: 35.75950633
*/