CII2187-一般任务安排问题的最小费用最大流解法及其他

PKU1031-求点对多边形所成张角

Matrix posted @ 2008年7月22日 00:00 in Problem Analysis with tags 点对多边形张角 计算几何 , 1863 阅读

by F. Yang     

    PKU1031题意描述:在坐标原点处有一个点光源S. 在距离S为r处的亮度为

I = k/r,

其中k为与光源有关的常数. 而光对在某处一高度为h长度为dl的小量竖直矩形的照度为

 dI = Ih \mid \cos \alpha \mid dl

其中I为该小量矩形处的亮度,\alpha为光源到dl处的方向与dl的法线方向的夹角.

    现有一高度为 h的简单多边形篱笆,用点列P_0 \dots P_{n-1}描述. 原点不在多边形的边上. 求光源对多边形各处的照度的总和.

     分析:dl的两端点为A, B. 且A, B在光源S到直线AB的垂足的同侧, 否则只需要将AB从垂足处分为两条线段讨论即可. 另外假设A离垂足较近. 设S到AB的垂线与SA的夹角为\theta, 而S到AB的小量张角为d\theta. SA的长度为r. 由小量分析的方法可得

 dl = rd\theta/\cos\theta.

    易知\theta即为上述的\alpha. 故得到在dl处的照度为

 dI = h \cdot\frac{k}{r}\cdot\cos\theta\cdot\frac{rd\theta}{\cos\theta} = hkd\theta.

    显然,照度与距离无关,只与张角有关. 故这个问题变成求一个点对一个简单多边形的张角.

     求点对多边形张角的问题可以用投影和区间并的思想解决:

(1) 令P_n = P_0;

(2 )维护一个闭区间序列ints[], 初始为空;

(3) 对每条条线段P_iP_{i+1}, 0 \leq i < n, 求出向量\overrightarrow{SP_i}\overrightarrow{SP_{i+1}}的辐角\alpha, \beta, 设\alpha \leq \beta. 若\beta - \alpha < \pi, 则将区间[\alpha, \beta]加入到序列ints[]中, 否则分别将区间[0, \alpha], [\beta, 2\pi]加入到序列ints[]中.

(4) 将区间序列ints[]按左端排序, 合并各个相交的区间;

(5) 最后剩下的不相交区间的长度总和就是所求夹角.

    下面为解决这个问题的代码, 函数接口为double angfun(pt p[], int n) , 由于本题光源刚好在原点,所以不当作参数, 对其他点也是容易处理的, 坐标平移一下即可.

  1. struct pt{
  2.         double x, y;
  3. };
  4.  
  5. struct interval{
  6.         double a, b;
  7. };
  8.  
  9. const double pi=acos(-1.0);
  10.  
  11. bool cmp(interval s1, interval s2)
  12. {
  13.         return s1.a<s2.a || s1.a==s2.a && s1.b<s2.b;
  14. }
  15.  
  16. double angle(pt p)
  17. {
  18.         double r, alpha;
  19.        
  20.         r=sqrt(p.x*p.x+p.y*p.y);
  21.         alpha=acos(p.x/r);
  22.         if(p.y<0) alpha=2*pi-alpha;
  23.         return alpha;
  24. }
  25.  
  26. double angfun(pt p[], int n)
  27. {
  28.         int i, m=0, t;
  29.         interval s[110];
  30.         double a, b, total;
  31.        
  32.         p[n]=p[0];
  33.        
  34.         for(i=0; i<n; i++){
  35.                 a=angle(p[i]);
  36.                 b=angle(p[i+1]);
  37.                 if(a>b) swap(a, b);
  38.                 if(b-a<pi-eps){
  39.                         s[m].a=a;
  40.                         s[m].b=b;
  41.                         m++;
  42.                 }
  43.                 else{
  44.                         s[m].a=0;
  45.                         s[m].b=a;
  46.                         m++;
  47.                         s[m].a=b;
  48.                         s[m].b=2*pi;
  49.                         m++;
  50.                 }
  51.         }
  52.        
  53.         sort(s, s+m, cmp);
  54.        
  55.         t=0;
  56.         for(i=1; i<m; i++)
  57.                 if(s[i].a<=s[t].b){
  58.                         if(s[t].b<s[i].b) s[t].b=s[i].b;
  59.                 }
  60.                 else s[++t]=s[i];
  61.                
  62.         total=0;
  63.         for(i=0; i<=t; i++) total+=s[i].b-s[i].a;
  64.        
  65.         return total;
  66. }
  67.  

    讨论: 假设允许点在多边形上, 这个方法也是可用的, 不过我们需要知道这个时候点是看作影响多边形外还是影响多边形内. 另外还需要判断每条边的左侧是内部还是右侧是内部. 这只要按叉积求一次面积即可判断. 具体在此从略. 本部分感谢bigheadghost提示.

    PS: PKU上该题做成多case会wrong到死, 要做成单case才能过, 晕死!

 

Avatar_small
bigheadghost 说:
2008年7月22日 03:37

拍了再看;p

Avatar_small
bigheadghost 说:
2008年7月22日 03:52

为什么最后会有不相交区间的?

难道意思是

要么剩下一个区间[a, b], 要么剩下两个[a, 2 pi] 和 [0, b]?

Avatar_small
bigheadghost 说:
2008年7月22日 03:52

and...

开始那段推导看不懂T_T

Head_small
Matrix 说:
2008年7月22日 05:27

 对!

为什么最后会有不相交区间的?

难道意思是

要么剩下一个区间[a, b], 要么剩下两个[a, 2 pi] 和 [0, b]?

Head_small
Matrix 说:
2008年7月22日 05:32

 

画图比较麻烦, 所以说得不太清楚, 呵呵...

小量分析里有一个方法: 顶角很小的等腰三角形近似于一个直角三角形, 而且两底角均看作直角. 底边长为rd\theta. 其中r为腰长, d\theta为顶角.

 

and...

开始那段推导看不懂T_T

 

Avatar_small
xmjt621 说:
2009年6月06日 23:59

@bigheadghost:
设光源到直线垂直距离h
l = h * tan(alpha)
=> d(l) = h * d(tan(alpha))
= h * (sec(alpha)) ^ 2 * d(alpha)
= h * 1 / (cos(alpha)) ^ 2 * d(alpha)
= r / cos(alpha) * d(alpha)

Avatar_small
walker01 说:
2010年4月10日 19:28

@xmjt621: Cool!Thanks a lot!


登录 *


loading captcha image...
(输入验证码)
or Ctrl+Enter