NOI 2009 描邊

問題簡述

平面上有一些筆畫,筆畫兩端爲半圓,中部爲矩形,筆畫可以互相重疊。任務是求平面上筆畫的總面積。

問題建模

這是一個計算幾何問題,我們可以把一個筆畫分解看爲一個矩形和兩個圓相互重疊(如下圖所示)。

image image

由於題中給定的是筆畫的兩個端點的座標和半徑,所以首先要求出矩形的四個端點座標。如下圖所示,設矩形兩端中點分別爲A,B,矩形A所在邊上一個頂點爲C。只需求出clip_image006的長度相等的平面法向量clip_image008,然後clip_image008[1]除以clip_image010求出單位向量,再乘以R就能求出clip_image012。用A,B座標分別加減clip_image012[1],即可求出矩形四個頂點的座標。

image

解法 豎線掃描

算法描述
用一條豎線去切割平面,會與上面的矩形或圓相交而且有一些相交的區間,把這些區間求並以後並乘以豎線的寬度,即可近似認爲是豎線覆蓋到的圖形面積。用豎線掃描完整個平面,其面積之和就是平面上圖形面積。掃描的豎線越細,結果就越精確。

image

我們要解決豎線與線段相交,以及豎線與圓相交的問題。

如圖所示,掃描線x=x0與線段AB相交,我們要求出交點D的縱座標。分別過B作平行於x=x0的直線,過A作垂直於x=x0的直線,相交於點C,AC交x=x0與點E。顯然有

clip_image018

所以

clip_image020

於是用E的座標加上clip_image022可算出D座標。

image

如圖是掃描線與圓A相交點C,作AB垂直於掃描線與點B,連接AC。|AC|=R,|AB|爲掃描線橫座標與點A橫座標之差絕對值,根據勾股定理即可求出|BC|,進而求出掃描線與圓的兩交點座標。

在程序實現時,我們需要求出每個“筆畫”的左右邊界橫座標,按照左邊界座標排序所有“筆畫”,並用一個雙向鏈表實現隊列維護當前掃描線可以掃描到的“筆畫”。掃描線向右移動時,判斷是否有新的“筆畫”左界小於等於掃描線橫座標,並加入隊列。同時,當發現一個“筆畫”右界小於掃描線橫座標,則把它從隊列中刪除。

實際測試
時間複雜度爲O(MC),M爲筆畫數,C爲掃描線數。在實際的測試中,前9個測試點掃描豎線寬度設定爲1*10^(-6),均得到了10分,程序運行時間如下表。

測試點 1 2 3 4 5 6 7 8 9
時間() 10 10 34 25 50 27 10 17 335

而第10個測試點由於規模過大,運行400秒僅能得到7分。當我把掃描豎線寬度設定爲5*10^(-8),,並運行了11天零4個小時後,終於獲得了10分。

參考程序

/* 
 * Problem: NOI2009 path
 * Author: Guo Jiabao
 * Time: 2009.10.9 13:50
 * State: Solved
 * Memo: 豎線掃描
*/
#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <cstring>
#include <ctime>
using namespace std;
const int MAXN = 2501;
const double INF = 1e100;

struct point
{
    double x,y;
    point (double tx,double ty) : x(tx),y(ty) {}
    point () : x(0) , y(0) {}
};

struct shape
{
    point c[2],v[4];
    double left,right;
};

struct itvs
{
    double a;
    bool s;
}isl[MAXN];

struct delist
{
    struct item
    {
        item *c[2];
        shape *e;
        item(item *last,shape *te)
        {
            c[0] = last;
            c[1] = NULL;
            e = te;
        }
    }*front,*rear;
    delist()
    {
        front = rear = 0;
    }
    void insert(shape *e)
    {
        if (rear)
            rear = rear -> c[1] = new item(rear,e);
        else
            front = rear = new item(NULL,e);
    }
    void remove(item *i)
    {
        if (i -> c[0])
            i -> c[0] -> c[1] = i -> c[1];
        else
            front = i -> c[1];
        if (i -> c[1])
            i -> c[1] -> c[0] = i -> c[0];
        else
            rear = i -> c[0];
        delete i;
    }
};

int M,N,ic;
double R,Total,Left,Right,miny,maxy;
point P[MAXN];
shape A[MAXN];
delist DL;

void operator *= (point &a,double k)
{
    a.x *=k;
    a.y *=k;
}

point operator + (point a,point b)
{
    return point(a.x + b.x,a.y + b.y);
}

point operator - (point a,point b)
{
    return point(a.x - b.x,a.y - b.y);
}

inline double MIN(double a,double b)
{
    return a<b?a:b;
}

inline double MAX(double a,double b)
{
    return a>b?a:b;
}

inline double ABS(double a)
{
    return a<0?-a:a;
}

inline double dist(point a,point b)
{
    return sqrt((a.x-b.x)*(a.x-b.x) + (a.y-b.y)*(a.y-b.y));
}

void init()
{
    int i,a,b;
    freopen("path.in","r",stdin);
    freopen("path.out","w",stdout);
    scanf("%d",&N);
    for (i=1;i<=N;i++)
        scanf("%lf%lf",&P[i].x,&P[i].y);
    scanf("%d",&M);
    for (i=1;i<=M;i++)
    {
        scanf("%d%d",&a,&b);
        A[i].c[0] = P[a];
        A[i].c[1] = P[b];
    }
    scanf("%lf",&R);
}

void count_point()
{
    int i;
    point v,v1,v2;
    Left = INF; Right = -INF;
    for (i=1;i<=M;i++)
    {
        v = A[i].c[1] - A[i].c[0];
        v1.x = v.y; v1.y = -v.x;
        v1 *= R / dist(A[i].c[0],A[i].c[1]);
        A[i].v[0] = A[i].c[0] + v1;
        A[i].v[1] = A[i].c[0] - v1;
        A[i].v[2] = A[i].c[1] + v1;
        A[i].v[3] = A[i].c[1] - v1;

        A[i].left = MIN(A[i].c[0].x-R,A[i].c[1].x-R);
        if (A[i].left < Left)
            Left = A[i].left;

        A[i].right = MAX(A[i].c[0].x+R,A[i].c[1].x+R);
        if (A[i].right > Right)
            Right = A[i].right;
    }
}

inline int cmp(const void *a,const void *b)
{
    return ((shape *)a) -> left < ((shape *)b) -> left ? -1 : 1;
}

inline int cmp1(const void *a,const void *b)
{
    if ( ((itvs *)a) -> a < ((itvs *)b) -> a )
        return -1;
    else if ( ((itvs *)a) -> a > ((itvs *)b) -> a )
        return 1;
    if (((itvs *)a) -> s)
        return -1;
    return 1;
}

void insert_ivl(double s,double t)
{
    if (s >= t)
        return;
    isl[++ic].a = s;
    isl[ic].s = true;
    isl[++ic].a = t;
    isl[ic].s = false;
}

double ivl_merge()
{
    int i,lv=0;
    double s,t=0;
    qsort(isl+1,ic,sizeof(isl[0]),cmp1);
    for (i=1;i<=ic;i++)
    {
        if (isl[i].s)
        {
            if (++lv == 1)
                s = isl[i].a;
        }
        else
        {
            if (--lv == 0)
                t += isl[i].a - s;
        }
    }
    return t;
}

inline bool isin(double x,double x1,double x2)
{
    if (x1 < x2)
        return x1 <= x && x <= x2;
    else
        return x2 <= x && x <= x1;
}

void cross_segment(double x,point a,point b)
{
    if (b.x == a.x || !isin(x,a.x,b.x))
        return;
    double y = (b.y - a.y) * (x - a.x) / (b.x - a.x);
    y += a.y;
    if (y < miny)
        miny = y;
    if (y > maxy)
        maxy = y;
}

void scan()
{
    const double dx = 1e-6;
    double x;
    double per,lper=0;
    int p=0;
    for (x=Left;x<=Right;x+=dx)
    {
        //////
        per = (x - Left)/(Right - Left) * 100;
        if (per - lper >=1)
        {
            lper = per;
            fprintf(stderr,"%.0lf percent\n",per);
        }
        //////
        while (p+1 <=M && A[p+1].left <= x)
            DL.insert(&A[++p]);
        ic = 0;
        for (delist :: item *it = DL.front;it;it = it ->c[1])
        {
            shape *S;
            for (;it;)
            {
                S = it -> e;
                if (x <= S->right)
                    break;
                delist :: item *itn = it -> c[1];
                DL.remove(it);
                it = itn;
            }
            if (!it)
                break;

            //Cross circle
            for (int i=0;i<2;i++)
            {
                if (S->c[i].x-R < x && x < S->c[i].x+R)
                {
                    double tx,ty;
                    tx = ABS(x - S->c[i].x);
                    ty = sqrt(R*R - tx*tx);
                    insert_ivl(S->c[i].y - ty , S->c[i].y + ty);
                }
            }

            miny = INF; maxy = -INF;
            cross_segment(x,S->v[0],S->v[1]);
            cross_segment(x,S->v[2],S->v[3]);
            cross_segment(x,S->v[0],S->v[2]);
            cross_segment(x,S->v[1],S->v[3]);
            insert_ivl(miny,maxy);
        }
        double t = ivl_merge();
        Total += t * dx;
    }
}

void solve()
{
    count_point();
    qsort(A+1,M,sizeof(A[0]),cmp);
    scan();
    printf("%.12lf\n",Total);
}

int main()
{
    int c_start,c_end;
    c_start = clock();
    init();
    solve();
    c_end = clock();
    cerr <<"Running Time : " << double(c_end - c_start) / CLOCKS_PER_SEC << endl;
    return 0;
}

相關日誌