首页 > 技术文章 > hihocoder 北大网络赛 E.

zhyfzy 2015-09-20 20:38 原文

给一个1000个点的多边形(从某个点依次按照外形给出每个节点),这个多边形不一定是凸多边形

再给一个圆,问这个多边形与圆相交区域的周长

 

我们将这个问题分成两个部分,第一部分是求线段在圆内的长度,第二部分是求圆弧的长度。

对于第一个部分,求线段在圆内的长度,

注意线段两端都在圆外边但是线段与圆有交点的情况。相切的情况可以不去考虑

 

第二个部分,我们在求第一部分时顺带求出所有线段与圆的交点,依次枚举每段圆弧,我们需要判断这段圆弧是否在多边形内

其实我们只需求出圆弧的中点,对于中点,判断是否在多边形内,射线法或者转角法均可,如果在多边形内就把其长度累加到答案中。

 

另外特殊判断一下整个圆都在多边形内部的情况

 

程序是赛后写出来的,还没有提交过_(:зゝ∠)_

#include<cstdio>
#include<cstring>
#include<cmath>
#include<iostream>
#include<algorithm>
#include<set>
#include<map>
#include<stack>
#include<vector>
#include<queue>
#include<string>
#include<sstream>
#define eps 1e-9
#define ALL(x) x.begin(),x.end()
#define INS(x) inserter(x,x.begin())
#define rep(i,j,k) for(int i=j;i<=k;i++)
#define MAXN 1005
#define MAXM 40005
#define INF 0x3fffffff
#define PB push_back
#define MP make_pair
#define X first
#define Y second
#define clr(x,y) memset(x,y,sizeof(x));
using namespace std;
typedef long long LL;
int i,j,k,n,m,x,y,T,big,cas,num,len;
bool flag;

const double pi=acos(-1.0);
int dcmp(double x) {
    if(fabs(x) < eps) return 0;
    else
    return x < 0 ? -1 : 1;
}
struct Vector
{
    double x, y;
    Vector (double x=0, double y=0) :x(x),y(y) {}
    Vector operator + (const Vector &B) const { return Vector (x+B.x,y+B.y); }
    Vector operator - (const Vector &B) const { return Vector(x - B.x, y - B.y); }
    Vector operator * (const double &p) const { return Vector(x*p, y*p); }
    Vector operator / (const double &p) const { return Vector(x/p, y/p); }
    double operator * (const Vector &B) const { return x*B.x + y*B.y;}//点积
    double operator ^ (const Vector &B) const { return x*B.y - y*B.x;}//叉积
    bool operator < (const Vector &b) const { return x < b.x || (x == b.x && y < b.y); }
    bool operator ==(const Vector &b) const { return dcmp(x-b.x) == 0 && dcmp(y-b.y) == 0; }
};
typedef Vector Point;
Point Read(){double x, y;scanf("%lf%lf", &x, &y);return Point(x, y);}
double Length(Vector A){ return sqrt(A*A); }//向量的模
double Angle(Vector A, Vector B){return acos(A*B / Length(A) / Length(B)); }//向量的夹角,返回值为弧度
double Area2(Point A, Point B, Point C){ return (B-A)^(C-A); }//向量AB叉乘AC的有向面积
Vector VRotate(Vector A, double rad){return Vector(A.x*cos(rad) - A.y*sin(rad), A.x*sin(rad) + A.y*cos(rad));}//向量A旋转rad弧度
Point  PRotate(Point A, Point B, double rad){return A + VRotate(B-A, rad);}//将B点绕A点旋转rad弧度
Vector Normal(Vector A){double l = Length(A);return Vector(-A.y/l, A.x/l);}//求向量A向左旋转90°的单位法向量,调用前确保A不是零向量

Point GetLineIntersection/*求直线交点,调用前要确保两条直线有唯一交点*/(Point P, Vector v, Point Q, Vector w){double t = (w^(P - Q)) / (v^w);return P + v*t;}//在精度要求极高的情况下,可以自定义分数类
double DistanceToLine/*P点到直线AB的距离*/(Point P, Point A, Point B){Vector v1 = B - A, v2 = P - A;return fabs(v1^v2) / Length(v1);}//不加绝对值是有向距离
double DistanceToSegment/*点到线段的距离*/(Point P, Point A, Point B)
{
    if (A==B) return Length(P-A);
    Vector v1=B-A,v2=P-A,v3=P-B;
    if (dcmp(v1*v2)<0) return Length(v2);else
    if (dcmp(v1*v3)>0) return Length(v3);else
    return fabs(v1^v2)/Length(v1);
}

Point GetLineProjection/*点在直线上的射影*/(Point P, Point A, Point B)
{
    Vector v=B-A;
    return A+v*((v*(P-A))/(v*v));
}

bool OnSegment/*判断点是否在线段上(含端点)*/(Point P,Point a1,Point a2)
{
    Vector v1=a1-P,v2=a2-P;
    if (dcmp(v1^v2)==0 && min(a1.x,a2.x)<=P.x  && P.x<=max(a1.x,a2.x)  && min(a1.y,a2.y)<=P.y && P.y<=max(a1.y,a2.y)) return true;
    return false;
}

bool SegmentInter/*线段相交判定*/(Point a1, Point a2, Point b1, Point b2)
{
    //if (OnSegment(a1,b1,b2) || OnSegment(a2,b1,b2) || OnSegment(b1,a1,a2) || OnSegment(b2,a1,a2)) return 1;
    //如果只判断线段规范相交(不算交点),上面那句可以删掉
    double c1=(a2-a1)^(b1-a1),c2=(a2-a1)^(b2-a1);
    double c3=(b2-b1)^(a1-b1),c4=(b2-b1)^(a2-b1);
    return dcmp(c1)*dcmp(c2)<0 && dcmp(c3)*dcmp(c4)<0;
}

bool InTri/*判断点是否在三角形内*/(Point P, Point a,Point b,Point c)
{
    if (dcmp(fabs((c-a)^(c-b))-fabs((P-a)^(P-b))-fabs((P-b)^(P-c))-fabs((P-a)^(P-c)))==0) return true;
    return false;
}

double PolygonArea/*求多边形面积,注意凸包P序号从0开始*/(Point *P ,int n)
{
    double ans = 0.0;
    for(int i=1;i<n-1;i++)
        ans+=(P[i]-P[0])^(P[i+1]-P[0]);
    return ans/2;
}
bool CrossOfSegAndLine/*判断线段是否与直线相交*/(Point a1,Point a2,Point b1,Vector b2)
{
    if (OnSegment(b1,a1,a2) || OnSegment(b1+b2,a1,a2)) return true;
    return dcmp(b2^(a1-b1))*dcmp(b2^(a2-b1))<0;
}

Point stand(Point u)
{
    double len=Length(u);
    return u/len;
}

int intersCL(Point cen,double R,Point A,Point B,Point &I1,Point &I2)//直线与圆的交点
{
    double dis=DistanceToLine(cen,A,B);
    int tt=dcmp(dis-R);
    if (dcmp(dis-R)>0) return 0;
    if (dcmp(dis-R)==0)
    {
        I1=GetLineProjection(cen,A,B);
        return 1;
    }

    double x=sqrt(R*R-dis*dis);
    Point p=GetLineProjection(cen,A,B);
    Vector v=stand(A-B)*x;
    I1=p+v;
    I2=p-v;
    return 2;
}

Point p[10005],c;
double R,ans;

bool cmp(Point a,Point b)
{
    return atan2(a.y,a.x)<atan2(b.y,b.x);
}

//转角发判定点P是否在多边形内部
int isPointInPolygon(Point P, Point* Poly, int n)
{
    int wn=0;
    for(int i = 0; i < n; ++i)
    {
        if(OnSegment(P, Poly[i], Poly[(i+1)%n]))    return -1;    //在边界上
        int k = dcmp((Poly[(i+1)%n] - Poly[i])^( P - Poly[i]));
        int d1 = dcmp(Poly[i].y - P.y);
        int d2 = dcmp(Poly[(i+1)%n].y - P.y);
        if(k > 0 && d1 <= 0 && d2 > 0)    wn++;
        if(k < 0 && d2 <= 0 && d1 > 0)    wn--;
    }
    if(wn != 0)    return 1;    //内部
    return 0;                //外部
}

Point middle (Point A,Point B)
{
    return Point((A.x+B.x)/2,(A.y+B.y)/2);
}
vector <Point> pt;

Point print(Point A)
{
    printf("%lf %lf\n",A.x,A.y);
}

int main()
{
    while (scanf("%d",&n),n)
    {
        for (i=0;i<n;i++) p[i]=Read();
        c=Read();scanf("%lf",&R);
        int tot=0;
        for (i=0;i<n;i++)
        {
            Vector v1=p[i]-p[(i-1+n)%n];
            Vector v2=p[(i+1)%n]-p[i];
            tot+=dcmp(v1^v2);
        }
        if (tot<0)
        {
            for (i=0;i<n/2;i++)
            {
                swap(p[i],p[n-i-1]);
            }
        }
        ans=0;
        pt.clear();
        for (i=0;i<n;i++)
        {
            Point A=p[i],B=p[(i+1)%n];
            Point I1,I2;
            int num=intersCL(c,R,A,B,I1,I2);
            if (num<=1) continue;


            bool T1=OnSegment(I1,A,B);
            bool T2=OnSegment(I2,A,B);
            bool T3=OnSegment(A,I1,I2);
            bool T4=OnSegment(B,I1,I2);

            if (A==I1) pt.PB(A-c);
            if (A==I2) pt.PB(A-c);
            if (B==I1) pt.PB(B-c);
            if (B==I2) pt.PB(B-c);

            if (T1 && T2)
            {
                ans+=Length(I2-I1);
                pt.PB(I1-c);pt.PB(I2-c);
            }else
            if (T3 && T4)
            {
                ans+=Length(A-B);
            }else
            if (!T1 && !T2)
            {
                continue;
            }else
            if (T1)
            {
                if (T4) ans+=Length(I1-B);
                else ans+=Length(I1-A);

                pt.PB(I1-c);
            }else
            {
                if (T4) ans+=Length(I2-B);
                else ans+=Length(I2-A);

                pt.PB(I2-c);
            }

        }
        sort(pt.begin(),pt.end(),cmp);
        unique(pt.begin(),pt.end(),cmp);
        
        int sz=pt.size();

        if (sz==0 && ans==0)
        {
            printf("%d\n",(int)round(pi*2*R));
            continue;
        }

        for (i=0;i<sz;i++)
        {
            Point p1=pt[i];
            Point p2;
            if (i==sz-1)
            {
                p2=pt[0];
            }else
                p2=pt[i+1];
            double a1=atan2(p1.y,p1.x);
            double a2=atan2(p2.y,p2.x);
            if (i==sz-1) a2+=2*pi;

            double a3=(a1+a2)/2;


            Point pm=VRotate(Point(R,0.0),a3);

            pm=stand(pm)*R+c;

            if (isPointInPolygon(pm,p,n)>0)
            {
                ans+=(a2-a1)*R;
            }
        }
        printf("%d\n",(int)round(ans));
    }
    return 0;
}

 

推荐阅读