三角剖分·圆和多边形的交

POJ 3675 Telescope

http://poj.org/problem?id=3675

大意：求解圆和多边形的交。

分析： 任意一个凸N多边形均可分解成N-2个三角形。因此，这就是讨论分解后的三角形和圆的交的问题。

它有这些情况:

（1）：

（2）：

（3）:

（4）:

(5):

5)又可分为p1在外，p2在里；p1在里，p2在外。

在写代码的过程中遇到了一些无法解决的问题，特别是5.2情况，参考了
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 {