HDU 5130 求圆和简单多边形公共部分面积

/*
模板：
求圆和简单多边形（可凸可凹，简单多边形是不自交的多边形）的公共面积
原理：
有序点集=多边形的点+圆和多边形的交点
然后以圆心为一点，每次按顺序（逆时针或者顺时针）取点集中两点，构成三角形，
比较三角形和对应弧形的绝对值面积大小，取绝对值较小的那块面积，按照逆时针为正面积，
顺时针为负面积求得所有三角形的有向面积和即可，最后取绝对值
*/
struct GetCirclePolyIntersectionArea{
Circle cir;
double Scir;

Point p[MAXN];
int tail;
GetCirclePolyIntersectionArea(){tail=0;}
GetCirclePolyIntersectionArea(Circle cir):cir(cir){
Scir = PI*cir.r*cir.r;
tail=0;
}
//tp[]是多边形的点集，n是点的个数。tp[]必须满足点是按顺时针或者逆时针排序的
void solve(Point tp[],int n){
for(int i=0;i<n;i++){
p[tail++]=tp[i];//p[]是囊括了圆和多边形交点的点集，也是按顺时针或者逆时针排序的
Line line = Line(tp[i],tp[(i+1)%n] - tp[i]);
double t1,t2;
vector<Point > sol;
sol.clear();
getLineCircleIntersection(line , cir , t1,t2,sol);

for(int j=0;j<sol.size();j++){
p[tail++]=sol[j];
}
}
double res=0;
for(int i=0;i<tail;i++){
Point O = cir.c;
double ang = Angle(p[(i+1)%tail]-O , p[i]-O);
if( dcmp( Cross( p[i]-O , p[(i+1)%tail]-O)) > 0 ) ang*=1;
else ang*=-1;
double Sshan = ang/(2*PI)*Scir;
double Strian = Area(O , p[i] ,p[(i+1)%tail] );
if(dcmp( abs(Sshan) - abs(Strian))<=0  ){
res += Sshan;
}else res += Strian;
}
printf("Case %d: %.10f\n",Cas++,abs(res));//.10f
}
};

/*完整代码*/
#include <iostream>
#include <cstdio>
#include <stack>
#include <cstring>
#include <queue>
#include <algorithm>
#include <cmath>
#include <map>
using namespace std;
#define N 300
typedef pair<int,int> PII;
int Cas=1;
const int INF=0x3f3f3f3f;
void Open()
{
#ifndef ONLINE_JUDGE
freopen("C:/OJ/in.txt","r",stdin);
//freopen("D:/my.txt","w",stdout);
#endif // ONLINE_JUDGE
}
const double eps = 1e-10;
const int MAXN=1111;
struct Point{
double x,y;
Point(double x=0 ,double y=0):x(x),y(y){}
void read(){
scanf("%lf%lf",&x,&y);
}
friend ostream &operator<< (ostream& output , Point& p){
output<<"("<<p.x<<","<<p.y<<")";
return output;
}

}tp[MAXN];
typedef Point Vector;

Vector operator - (Point A, Point B){
return Vector(A.x - B.x , A.y - B.y);
}
Vector operator + (Point A, Point B){
return Vector(A.x + B.x , A.y + B.y);
}
Vector operator * (Point A, double p){
return Vector(A.x * p , A.y * p);
}
Vector operator / (Point A, double p){
return Vector(A.x / p , A.y / p);
}

double Dot(Vector A, Vector B){
return A.x*B.x + A.y*B.y;
}
double Cross(Vector A, Vector B){
return A.x*B.y - A.y*B.x;
}
double Length(Vector A){
return sqrt(Dot(A,A));
}
double Area(Point A,Point B,Point C){
return 0.5*Cross(B-A,C-A);
}
double Angle(Vector A, Vector B){
return acos(Dot(A,B) / Length(A) / Length(B) );
}

int dcmp(double x){
if(x>eps) return 1;
else if(x<-eps) return -1;
else return 0;
}

struct Line{
Point P;
Vector v;
double ang;
Line(){}
Line(Point P, Vector v):P(P) ,v(v){ang = atan2(v.y , v.x);}
bool operator < (const Line& L) const{
return ang<L.ang;
}
Point point(double t){return P+v*t;}
};
struct Circle{
Point c;
double r;
Circle(){}
Circle(Point c, double r):c(c) , r(r){}
Point point(double a){
return Point(c.x + cos(a)*r ,c.y + sin(a)*r);
}
};

int getLineCircleIntersection(Line L ,Circle C, double& t1, double& t2 , vector<Point>& sol){
double a = L.v.x , b = L.P.x - C.c.x , c = L.v.y , d=L.P.y - C.c.y;
double e = a*a + c* c , f= 2*(a*b + c*d) , g=b*b + d*d - C.r*C.r;
double delta = f*f -4*e*g;
if(dcmp(delta ) < 0 ) return 0;
if(dcmp(delta) == 0 ){
t1 = t2 = -f/(2*e);
if(dcmp(t1)>0 && dcmp(t1-1)<0)
sol.push_back(L.point(t1));
return 1;
}
t1 = (-f - sqrt(delta)) / (2*e);
t2 = (-f + sqrt(delta)) / (2*e);

if(t1>t2) swap(t1,t2);

if(dcmp(t1)>0 && dcmp(t1-1)<0)
sol.push_back(L.point(t1));
if(dcmp(t2)>0 && dcmp(t2-1)<0)
sol.push_back(L.point(t2));
//    if(sol.size()==2){
//        cerr<<"2 sol:";
//        cerr<<t1<<" "<<t2<<endl;
//    }
return (int)sol.size();
}

Point p[MAXN+MAXN];
int tail=0;
double DOU(double x){
return x*x;
}
const double PI = acos(-1.0);

struct GetCirclePolyIntersectionArea{
Circle cir;
double Scir;

Point p[MAXN];
int tail;
GetCirclePolyIntersectionArea(){tail=0;}
GetCirclePolyIntersectionArea(Circle cir):cir(cir){
Scir = PI*cir.r*cir.r;
tail=0;
}
//tp[]是多边形的点集，n是点的个数。tp[]必须满足点是按顺时针或者逆时针排序的
void solve(Point tp[],int n){
for(int i=0;i<n;i++){
p[tail++]=tp[i];//p[]是囊括了圆和多边形交点的点集，也是按顺时针或者逆时针排序的
Line line = Line(tp[i],tp[(i+1)%n] - tp[i]);
double t1,t2;
vector<Point > sol;
sol.clear();
getLineCircleIntersection(line , cir , t1,t2,sol);

for(int j=0;j<sol.size();j++){
p[tail++]=sol[j];
}
}
double res=0;
for(int i=0;i<tail;i++){
Point O = cir.c;
double ang = Angle(p[(i+1)%tail]-O , p[i]-O);
if( dcmp( Cross( p[i]-O , p[(i+1)%tail]-O)) > 0 ) ang*=1;
else ang*=-1;
double Sshan = ang/(2*PI)*Scir;
double Strian = Area(O , p[i] ,p[(i+1)%tail] );
if(dcmp( abs(Sshan) - abs(Strian))<=0  ){
res += Sshan;
}else res += Strian;
}
printf("Case %d: %.10f\n",Cas++,abs(res));//.10f
}
};

int main()
{
//Open();
int n;
double K;

while(scanf("%d%lf",&n,&K)==2){
tail=0;
for(int i=0;i<n;i++){
tp[i].read();
}
double ax,ay,bx,by;
scanf("%lf%lf%lf%lf",&ax,&ay,&bx,&by);

double mu = K*K-1;
Point O=Point(-(bx - K * K * ax) / mu , -(by - K * K * ay) / mu);
double R = sqrt( DOU((bx - K * K * ax)/(mu)) + DOU((by - K * K * ay) / (mu)) - (K * K * (ax*ax + ay*ay)) / mu + (bx*bx + by*by) / mu );
Circle cir = Circle(O,R);

GetCirclePolyIntersectionArea gcp=GetCirclePolyIntersectionArea(cir);
gcp.solve(tp,n);
}
return 0;
}