HDU 3961 Crisis【恶心的几何模拟】

链接:http://acm.hdu.edu.cn/showproblem.php?pid=3961
这题竟然有队伍场上过掉了。。大神们的脑子得有多清楚。。。我这题是做了一天才理清思路,过了样例。。这题的样例就已然强的一B了。。

题目大意即做法,题意很裸,直接把做法说出来了,所以就是模拟题咯。
先给你一个圆和圆外一点,然后给你一个多边形(不一定凸凹),将多边形绕圆心公转,同时自转,给你转动周期和时间。
让你求圆外点向圆发出光线,除去被多边形挡住的部分后,被照弧的弧长有多少。

就模拟呗:
用圆外点求到圆的两条公切线的切点,解一个二元二次方程和一个二元一次方程联立。
(x0,y0)圆心坐标 (xp,yp)圆外点坐标
(xp - x0)(x - x0) + (yp - y0)(y - y0) = r^2
(x - x0)^2 + (y - y0)^2 = r^2
得到切点cutl,cutr(根据cmp函数判断一下两者相对关系,为之后极角排序找点比较做准备)
仿射变换公式:(x,y) 转theta角到(x*cos(theta) - y*sin(theta),x*sin(theta)+y*cos(theta))
求多边形重心,记录一下重心到各顶点的向量,将重心公转到指定位置(应用二维仿射变换),再用向量还原多边形,再绕重心自转多边形。
此时对多边形上的进行极角排序,找出两端的点L,R。
判断找到的两点是否在(xp,yp),cutl,cutr 构成的三角形内。我用了个麻烦的方法,是先根据叉积的左右手性判断一部分情况的,后续说明。
再用(xp,yp)和L,R求直线方程,再用直线方程联立圆方程求解交点坐标。
ax+by+c=0
(x-x0)^2 + (y-y0)^2 = r^2
紧跟前步,用交点坐标到(xp,yp)的距离和L、R比较来判断多边形是否在圆后。
唉,这个好麻烦啊写的。。还是天王的方法好,直接用(xp,yp)、cutl、cutr、(x0,y0)构成四边形,判断交点是否在四边形内即可。
最后用点积除以模长再求acos的方法计算圆心角和弧长,再相减即可。
圆心角alpha = acos(vectorL * vectorR / (|vectorL|*|vectorR|))
注意用acos的时候要调整精度,因为除法可能会得出一个大于1的数字,这样acos就bug了。。。精度伤不起啊。
代码随手贴一下,写了一天半的代码,太矬了,弱爆了:
#include <iostream>
#include 
<cstdio>
#include 
<cmath>
#include 
<cstring>
#include 
<algorithm>
#define maxn 105
using namespace std;
const double pi = acos(-1.0);
struct point
{
    
double x,y;
    point(){}
    point(
double a,double b):x(a),y(b){}
    point 
operator - (const point &p)
    {
        
return point(x - p.x,y - p.y);
    }
    point 
operator + (const point &p)
    {
        
return point(x + p.x,y + p.y);
    }
    
double operator * (const point &p)
    {
        
return x * p.x + y * p.y;
    }
    
double operator ^ (const point &p)
    {
        
return x * p.y - y * p.x;
    }
    
double norm()
    {
        
return sqrt(x * x + y * y);
    }
}sun,earth,vec[maxn],p[maxn],cutl,cutr;
double multi(point a,point b,point o)
{
    
return (a - o) ^ (b - o);
}
double r;
int T1,T2,T;
int n;
const double eps = 1e-6;
int comp(double x)
{
    
if(fabs(x) < eps)
        
return 0;
    
else if(x < -eps)
        
return -1;
    
else
        
return 1;
}
point rotate(point a,point o,
double theta)
{
    point vec 
= a - o;
    point to 
= point(vec.x * cos(theta) - vec.y * sin(theta),vec.x * sin(theta) + vec.y * cos(theta));
    
return o + to;
}
double tarea(point p0,point p1,point p2)
{
    
double k = p0.x * p1.y + p1.x * p2.y + p2.x * p0.y - p1.x * p0.y - p2.x * p1.y - p0.x * p2.y;
    
return k / 2;
}
bool cmp(point p,point q)
{
    
double t1 = atan2(p.y - sun.y,p.x - sun.x),t2 = atan2(q.y - sun.y,q.x - sun.x);
    
return (comp(t1 - t2) > 0);
}
struct line
{
    
double a,b,c;
    line(){}
    line(point p,point q)
    {
        
if(comp(p.x - q.x) == 0)
        {
            a 
= 1.0;
            b 
= 0.0;
            c 
= -p.x;
        }
        
else if(comp(p.y - q.y) == 0)
        {
            a 
= 0.0;
            b 
= 1.0;
            c 
= -p.y;
        }
        
else
        {
            a 
= 1.0;
            b 
= -(p.x - q.x) / (p.y - q.y);
            c 
= -(p.x + b * p.y);
        }
    }
};
bool arc(line llll,point &P)
{
    
double a = llll.a,b = llll.b,c = llll.c;
    
double A = b * b + 1,B = 2.0 * b * c + 2.0 * b * earth.x - 2.0 * earth.y,C = c * c + 2.0 * c * earth.x + earth.x * earth.x + earth.y * earth.y - r * r;
    
double det = B * B - 4.0 * A * C;
    
if(comp(det) < 0)
        
return false;
    
else if(comp(det) == 0)
    {
        
double x,y;
        y 
= (-+ sqrt(det)) / (2.0 * A);
        x 
= -(b * y + c);
        point w(x,y);
        P 
= w;
        
return true;
    }
    
else
    {
        
double y1 = (-+ sqrt(det))/(2.0 * A),y2 = (-- sqrt(det)) / (2.0 * A);
        
double x1 = -(b * y1 + c),x2 = -(b * y2 + c);
        point w1(x1,y1),w2(x2,y2);
        
double dis1 = (w1 - sun).norm(),dis2 = (w2 - sun).norm();
        
if(comp(dis1 - dis2) > 0)
            P 
= w2;
        
else
            P 
= w1;
        
return true;
    }
}
void calc()
{
    
double a,b,c;
    a 
= (sun.y - earth.y) * (sun.y - earth.y) + (sun.x - earth.x) * (sun.x - earth.x);
    b 
= -2.0 * earth.y * ((sun.y - earth.y) * (sun.y - earth.y) + (sun.x - earth.x) * (sun.x - earth.x)) - 2.0 * r * r * (sun.y - earth.y);
    c 
= r * r * r * r + 2.0 * r * r * (sun.y - earth.y) * earth.y - r * r * (sun.x - earth.x) * (sun.x - earth.x) + earth.y * earth.y * ((sun.y - earth.y) * (sun.y - earth.y) + (sun.x - earth.x) * (sun.x - earth.x));
    
double det = b * b - 4.0 * a * c;
    
double y1 = (-+ sqrt(det))/(2.0 * a),y2 = (-- sqrt(det)) / (2.0 * a);
    
double x1 = (r * r - (sun.y - earth.y) * (y1 - earth.y)) / (sun.x - earth.x) + earth.x;
    
double x2 = (r * r - (sun.y - earth.y) * (y2 - earth.y)) / (sun.x - earth.x) + earth.x;
    cutl 
= point(x1,y1),cutr = point(x2,y2);
    
if(!cmp(cutl,cutr))
        swap(cutl,cutr);
}
double len(point a,point b)
{
    point vec1 
= point(a - earth),vec2 = point(b - earth);
    
double value = (vec1 * vec2) / (vec1.norm() * vec2.norm());
    
if(comp(value - 1.0>= 0)
        value 
= 1.0;
    
double theta = acos(value);
    
return r * theta;
}
int main()
{
    
while(scanf("%lf %lf %lf %lf %lf",&sun.x,&sun.y,&earth.x,&earth.y,&r) == 5)
    {
        
//calc arc
        calc();
        scanf(
"%d %d %d %d",&n,&T1,&T2,&T);
        
for(int i = 0;i < n;i++)
        {
            
double x,y;
            scanf(
"%lf %lf",&x,&y);
            p[i] 
= point(x,y);
        }
        
int t1 = T % T1,t2 = T % T2;
        
double gong = 2.0 * pi * (double)t1 / (double)T1,zi = 2.0 * pi * (double) t2 / (double) T2;
        point zhong,temp(
0,0);
        
double area = 0,sumx = 0,sumy = 0,sarea = 0;
        
for(int i = 1;i < n - 1;i++)
        {
            temp.x 
= p[0].x + p[i].x + p[i+1].x;
            temp.y 
= p[0].y + p[i].y + p[i+1].y;
            area 
= tarea(p[0],p[i],p[i+1]);
            sarea 
+= area;
            sumx 
+= temp.x * area;
            sumy 
+= temp.y * area;
        }
        zhong 
= point(sumx / sarea / 3,sumy / sarea / 3);
        
for(int i = 0;i < n;i++)
            vec[i] 
= p[i] - zhong;
        
//rotate the gravity point
        zhong = rotate(zhong,earth,gong);
        
//gongzhuan
        for(int i = 0;i < n;i++)
            p[i] 
= vec[i] + zhong;
        
//zizhuan
        for(int i = 0;i < n;i++)
            p[i] 
= rotate(p[i],zhong,zi);
        sort(p,p 
+ n,cmp);
        
if(!cmp(cutl,cutr))
            swap(cutl,cutr);
        point L 
= p[0],R = p[n - 1];
        
if(!cmp(L,R))
            swap(L,R);
        point xl,xr;
        line LE 
= line(sun,L),RI = line(sun,R);
        
if(comp(multi(L,cutl,sun)) <= 0)
            xl 
= cutl;
        
else if(comp(multi(L,cutr,sun)) >= 0)
            xl 
= cutr;
        
else
            arc(LE,xl);
        
if(comp(multi(R,cutl,sun)) <= 0)
            xr 
= cutl;
        
else if(comp(multi(R,cutr,sun)) >= 0)
            xr 
= cutr;
        
else
            arc(RI,xr);
        
double len1 = len(cutl,cutr),len2 = len(xl,xr);
        
if(comp((L - sun).norm() - (xl - sun).norm()) > 0 || comp((R - sun).norm() - (xr - sun).norm()) > 0)
            len2 
= 0;
        
double ans = len1 - len2;
        printf(
"%.2lf\n",ans);
    }
}

posted on 2011-08-24 21:26 BUPT-[aswmtjdsj] @ Penalty 阅读(240) 评论(0)  编辑 收藏 引用 所属分类: 计算几何HDU Solution Report


只有注册用户登录后才能发表评论。
【推荐】超50万行VC++源码: 大型组态工控、电力仿真CAD与GIS源码库
网站导航: 博客园   IT新闻   BlogJava   知识库   博问   管理


<2011年8月>
31123456
78910111213
14151617181920
21222324252627
28293031123
45678910

导航

统计

常用链接

留言簿(1)

随笔分类(150)

随笔档案(71)

搜索

积分与排名

最新评论

阅读排行榜

评论排行榜