HDU 4667 Building Fence 計算幾何 凸包+圓

1.三角形的所有端點

2.過所有三角形的端點對所有圓做切線，得到所有切點。

3.做任意兩圓的外公切線，得到所有切點。

對上述所有點求凸包，標記每個點是三角形上的點還是某個圓上的點。

求完凸包后，因為所有點都是按逆時針（或順時針）排好序的，如果相鄰兩點在同一圓上，那么求這段圓弧的距離，否則求這段直線的距離。最后得到所有周長。

 

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <algorithm>

using namespace std;

const double eps = 1e-9;
const double PI = acos(-1.0);
const int MAXN = 60;

struct Point
{
    double x, y;
    int id;         //點標號，標記是否在同一個圓上
    Point() { }
    Point( double x, double y ):x(x), y(y) { }
    Point( double x, double y, int id ):x(x), y(y), id(id) { }
    void readPoint()
    {
        scanf( "%lf%lf", &x, &y );
        return;
    }
};

struct Circle
{
    Point c;   //圓心坐標
    double r;  //半徑
    Circle() {}
    Circle( Point c, double r ): c(c), r(r) {}
    Point getPoint( double theta )   //根據極角返回圓上一點的坐標
    {
        return Point( c.x + cos(theta)*r, c.y + sin(theta)*r );
    }
    void readCircle()
    {
        scanf("%lf%lf%lf", &c.x, &c.y, &r );
        return;
    }
};

typedef Point Vector;

Vector operator+( Vector A, Vector B )       //向量加
{
    return Vector( A.x + B.x, A.y + B.y );
}

Vector operator-( Vector A, Vector B )       //向量減
{
    return Vector( A.x - B.x, A.y - B.y );
}

Vector operator*( Vector A, double p )      //向量數乘
{
    return Vector( A.x * p, A.y * p );
}

Vector operator/( Vector A, double p )      //向量數除
{
    return Vector( A.x / p, A.y / p );
}

int dcmp( double x )    //控制精度
{
    if ( fabs(x) < eps ) return 0;
    else return x < 0 ? -1 : 1;
}

bool operator<( const Point& A, const Point& B )   //兩點比較
{
    return dcmp( A.x - B.x) < 0 || ( dcmp(A.x - B.x ) == 0 && dcmp( A.y - B.y ) < 0 );
}

bool operator>( const Point& A, const Point& B )   //兩點比較
{
    return dcmp( A.x - B.x) > 0 || ( dcmp(A.x - B.x ) == 0 && dcmp( A.y - B.y ) > 0 );
}

bool operator==( const Point& a, const Point& b )   //兩點相等
{
    return dcmp( a.x - b.x ) == 0 && dcmp( a.y - b.y ) == 0;
}

double Cross( Vector A, Vector B )   //向量叉積
{
    return A.x * B.y - A.y * B.x;
}

double PointDis( Point a, Point b ) //兩點距離的平方
{
    return (a.x - b.x)*(a.x - b.x) + (a.y - b.y)*(a.y - b.y);
}

//求凸包，graham算法，O(nlogn)，返回凸包點的個數
int graham( Point *p, int n, Point *ch )
{
    if ( n <= 2 ) return 0;
    int top = 0;
    sort( p, p + n );

    ch[ top ] = p[0];
    ch[ ++top ] = p[1];
    ch[ ++top ] = p[2];

    top = 1;

    for ( int i = 2; i < n; ++i )
    {
        while ( top && dcmp( Cross( ch[top] - ch[top - 1], p[i] - ch[top - 1] ) ) <= 0 ) --top;
        ch[++top] = p[i];
    }
    int len = top;
    ch[++top] = p[n - 2];
    for ( int i = n - 3; i >= 0; --i )
    {
        while ( top > len && dcmp( Cross( ch[top] - ch[top - 1], p[i] - ch[top - 1] ) ) <= 0 ) --top;
        ch[++top] = p[i];
    }
    return top;
}

//過定點做圓的切線,得到切點，返回切點個數
//tps保存切點坐標
int getTangentPoints( Point p, Circle C, Point *tps )
{
    int cnt = 0;

    double dis = sqrt( PointDis( p, C.c ) );
    int aa = dcmp( dis - C.r );
    if ( aa < 0 ) return 0;  //點在圓內
    else if ( aa == 0 ) //點在圓上，該點就是切點
    {
        tps[cnt] = p;
        ++cnt;
        return cnt;
    }

    //點在圓外，有兩個切點
    double base = atan2( p.y - C.c.y, p.x - C.c.x );
    double ang = acos( C.r / dis );
    //printf( "base = %f ang=%f\n", base, ang );
    //printf( "base-ang=%f  base+ang=%f \n", base - ang, base + ang );

    tps[cnt] = C.getPoint( base - ang ), ++cnt;
    tps[cnt] = C.getPoint( base + ang ), ++cnt;

    return cnt;
}

//求兩圓外公切線切點，返回切線個數
//p是圓c2在圓c1上的切點
int makeCircle( Circle c1, Circle c2, Point *p )
{
    int cnt = 0;
    double d = sqrt( PointDis(c1.c, c2.c) ), dr = c1.r - c2.r;
    double b = acos(dr / d);
    double a = atan2( c2.c.y - c1.c.y, c2.c.x - c1.c.x );
    double a1 = a - b, a2 = a + b;
    p[cnt++] = Point(cos(a1) * c1.r, sin(a1) * c1.r) + c1.c;
    p[cnt++] = Point(cos(a2) * c1.r, sin(a2) * c1.r) + c1.c;
    return cnt;
}

double DisOnCircle( Point a, Point b, Circle c )  //求圓上一段弧長
{
    Point o = c.c;
    double A = sqrt( PointDis( o, a ) );
    double B = sqrt( PointDis( o, b ) );
    double C = sqrt( PointDis( a, b ) );
    double alpha = acos( ( A*A + B*B - C*C ) / ( 2.0*A*B ) );
    if ( dcmp( Cross( o-a, o-b ) ) < 0 ) return alpha * c.r;
    else return ( 2.0*PI - alpha ) * c.r;
}

/**********************以上模板**********************/

int cntC, cntT;  //圓的個數，三角形的個數
Circle yuan[MAXN];   //所有圓
Point PP[300100];    //所有點
Point tubao[300100]; //凸包
int totPP;           //點總數

void showP( Point *p, int nn )
{
    printf( "allP = %d\n", nn );
    for ( int i = 0; i < nn; ++i )
        printf("%f %f\n", p[i].x, p[i].y );
    puts("-------------------------");
    return;
}

int main()
{
    //freopen( "10022.in", "r", stdin );
    //freopen( "s.out", "w", stdout );
    while ( scanf( "%d%d", &cntC, &cntT ) == 2 )
    {
        totPP = 0;
        for ( int i = 0; i < cntC; ++i )
            yuan[i].readCircle();
        for ( int i = 0; i < cntT; ++i )
        {
            for ( int j = 0; j < 3; ++j )
            {
                PP[totPP].readPoint();
                PP[totPP].id = -(totPP+2);
                ++totPP;
            }
        }

        if ( cntC == 1 && cntT == 0 )
        {
            printf("%.6f\n", 2.0 * PI * yuan[0].r );
            continue;
        }

        int pretot = totPP;
        //求兩圓的外切點
        for ( int i = 0; i < cntC; ++i )
            for ( int j = i + 1; j < cntC; ++j )
            {
                Point PonA[4], PonB[4];
                makeCircle( yuan[i], yuan[j], PonA );
                int ans = makeCircle( yuan[j], yuan[i], PonB );
                for ( int k = 0; k < ans; ++k )
                {
                    PonA[k].id = i;
                    PonB[k].id = j;
                    PP[totPP++] = PonA[k];
                    PP[totPP++] = PonB[k];
                }
            }

        //求所有點與所有圓的切點
        for ( int i = 0; i < pretot; ++i )
        {
            for ( int j = 0; j < cntC; ++j )
            {
                Point qiedian[4];
                int ans = getTangentPoints( PP[i], yuan[j], qiedian );
                for ( int k = 0; k < ans; ++k )
                {
                    qiedian[k].id = j;
                    PP[totPP++] = qiedian[k];
                }
            }
        }

        //showP( PP, totPP );
        int cntBao = graham( PP, totPP, tubao );
        //puts("*********");
        //showP( tubao, cntBao );
        double girth = 0.0;
        tubao[cntBao] = tubao[0];

        for ( int i = 1; i <= cntBao; ++i )
        {
            if ( tubao[i].id == tubao[i - 1].id )  //如果兩點在同一個圓上
                girth += DisOnCircle( tubao[i], tubao[i - 1], yuan[ tubao[i].id ] );
            else
                girth += sqrt( PointDis( tubao[i], tubao[i - 1] ) );

        }

        printf( "%.5lf\n", girth );
    }
    return 0;
}