POJ 3675 Telescope
大意:求解圆和多边形的交。
分析:
任意一个凸N多边形均可分解成N-2个三角形。因此,这就是讨论分解后的三角形和圆的交的问题。
它有这些情况:
(1):
(2):
(3):
(4):
(5):
5)又可分为p1在外,p2在里;p1在里,p2在外。
在写代码的过程中遇到了一些无法解决的问题,特别是5.2情况,参考了
http://gzhu-101majia.iteye.com/blog/1128345
http://gzhu-101majia.iteye.com/blog/1128345
不能理解为什么自己的做法不能通过。
求解p1p2直线和圆的交点(假设不与横轴垂直,可参照(5)图):
分解后的三角形和圆的交按有向面积累加:
#include <iostream>
#include <stdio.h>
#include <cmath>
using namespace std;
const double eps=1e-7,PI=acos(-1.0);
struct point
{
double x,y;
}p[55];
double r;
double pp_dis(point p1,point p2) //两点距离
{
return sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y));
}
double xmulti(point p0,point p1,point p2){
return (p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y);
}
double angle(point p1,point o,point p2){
double a=pp_dis(o,p1);
double b=pp_dis(o,p2);
double c=pp_dis(p1,p2);
double cos=(a*a+b*b-c*c)/(2*a*b);
return fabs(acos(cos));
}
double dir(point o,point p1,point p2){
double area=xmulti(o,p1,p2);
if(area<-eps) return 1.0;
return -1.0;
}
point intersec(double x,double y) //计算直线与圆的交点
{
double k;
point temp;
if(x!=0)
{
k=y/x;
temp.x=fabs(r)/sqrt(1+k*k);
if(x<0) temp.x=-temp.x;
temp.y=k*temp.x;
}
else
{
temp.x=0;
if(y>0) temp.y=r;
else temp.y=-r;
}
return temp;
}
double get_area(point o,point p1,point p2) //三角剖分
{
double f=dir(o,p1,p2); //判断三角形面积加还是减
double s;
double a=pp_dis(o,p1);
double b=pp_dis(o,p2);
double c=pp_dis(p1,p2);
double d=fabs(xmulti(o,p1,p2))/c;
// 1
if(a<=r && b<=r)
{
double area=xmulti(o,p1,p2);
return fabs(area)/2.0*f;
}
// 2
else if(a>=r && b>=r && d>=r)
{
double sita1=angle(p1,o,p2);
double s=fabs(sita1*r*r/2.0); //扇形s=θ*r*r/2
return s*f;
}
// 3
else if(a>=r && b>=r && d<=r && (angle(o,p1,p2)-PI/2>eps || angle(o,p2,p1)-PI/2>eps))
{
double sita=angle(p1,o,p2);
s=fabs(sita*r*r/2.0);
return s*f;
}
// 4
else if(a>=r && b>=r && d<=r && angle(o,p1,p2)<=PI/2 && angle(o,p2,p1)<=PI/2)
{
point p3,p4;
if(fabs(p1.x-p2.x)>eps){
double k=(p2.y-p1.y)/(p2.x-p1.x);
double A=1+k*k;
double B=2*p1.y*k-2*k*k*p1.x;
double C=k*k*p1.x*p1.x+p1.y*p1.y-2*p1.y*k*p1.x-r*r;
double x1=(-B-sqrt(B*B-4*A*C))/(2*A);
double x2=(-B+sqrt(B*B-4*A*C))/(2*A);
if(fabs(x1-p1.x)<fabs(x2-p1.x)){
p3.x=x1;
p4.x=x2;
}
else {