ZOJ 2675 Little Mammoth（圆和矩形的交——三角剖分）

题目：ZOJ   2675  Little Mammoth 
http://acm.zju.edu.cn/onlinejudge/showProblem.do?problemCode=2675
大意：求解给定圆和矩形的交。

分析：三角剖分的应用，继上一篇博文说，这次使用那种容易理解的方法来做，不用那个吓人的模板，嗯嗯，正常工作。感人啊，给人继续做题的勇气。。。

#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[5];
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 {
                 p3.x=x2;
                 p4.x=x1;
             }
             p3.y=p1.y+k*(p3.x-p1.x);
             p4.y=p1.y+k*(p4.x-p1.x);
         }
         else {
             p3.x=p1.x;
             p4.x=p1.x;
             double h=sqrt(r*r-d*d);
             if(p1.y>eps){
                 p3.y=h;
                 p4.y=-h;
             }
             else {
                 p3.y=-h;
                 p4.y=h;
             }
         }

         double ag1=angle(p1,o,p3);
         double ag2=angle(p4,o,p2);
         double area=0;
         area+=r*r*(ag1+ag2)/2;
         area+=fabs(xmulti(p3,o,p4))/2.0;
         return area*f;
    }

    // 5.1
    else if(a>=r && b<=r)
    {
         double area=0;
         point p3;
         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.x*k*k+2*p1.y*k;
             double C=p1.y*p1.y+k*k*p1.x*p1.x-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( (x1>p1.x&&x1<p2.x) || (x1>p2.x&&x1<p1.x)) p3.x=x1;
             else p3.x=x2;
             p3.y=p1.y+k*(p3.x-p1.x);
         }
         else {
             p3.x=p1.x;
             double y1=sqrt(r*r-p3.x*p3.x);
             double y2=-y1;
             if( (y1>p1.y&&y1<p2.y) || (y1>p2.y&&y1<p1.y)) p3.y=y1;
             else p3.y=y2;
         }
         area+=fabs(xmulti(p3,o,p2))/2.0;
         double ag=angle(p3,o,p1);
         area+=r*r*ag/2;
         return area*f;
    }
    //  5.2
    else if(a<=r && b>=r)    
    {
        double area=0;
        point p3;
        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.x*k*k+2*p1.y*k;
             double C=p1.y*p1.y+k*k*p1.x*p1.x-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( (x1>p1.x&&x1<p2.x) || (x1>p2.x&&x1<p1.x)) p3.x=x1;
             else p3.x=x2;
             p3.y=p1.y+k*(p3.x-p1.x);
        }
        else
        {
             p3.x=p1.x;
             double y1=sqrt(r*r-p3.x*p3.x);
             double y2=-y1;
             if( (y1>p1.y&&y1<p2.y) || (y1>p2.y&&y1<p1.y)) p3.y=y1;
             else p3.y=y2;
        }
        //double sita1=angle(p2,o,p3);
        point pp2=intersec(p2.x,p2.y);
        double c=pp_dis(pp2,p3);
        double cos=(r*r+r*r-c*c)/(2*r*r);
        double sita1=acos(cos);  
        double s1=fabs(sita1*r*r/2.0); 
        double s3=fabs(p3.x*p1.y-p3.y*p1.x)/2.0;
        area=s1+s3;
        return area*f;
    }
    else return 0;
}

int main()
{
    //freopen("cin.txt","r",stdin);
    point o;
    while(~scanf("%lf%lf",&o.x,&o.y)){
        scanf("%lf",&r);
        double x1,y1,x2,y2;
        scanf("%lf%lf%lf%lf",&x1,&y1,&x2,&y2);
        p[0].x=x1-o.x;  p[0].y=y2-o.y;
        p[1].x=x2-o.x;  p[1].y=y2-o.y;
        p[2].x=x2-o.x;  p[2].y=y1-o.y;
        p[3].x=x1-o.x;  p[3].y=y1-o.y;
        p[4]=p[0];
        o.x=o.y=0;
        double ans=0;
        for(int i=0;i<4;i++)
        {
            ans+=get_area(o,p[i],p[i+1]);
        }
        printf("%.12lf\n",fabs(ans));
    }
    return 0;
}