半平面交

今天比较颓废。。不想看算法不想做题目。。。于是乎把半天的青春喂给老题半平面交了。

 

不过话说以前A poj1279的时候是各种打补丁的修改+提交+WA啊。。。这样怎么行。。。于是今天理清了下这题的思路。

 

那篇著名的论文我也是拜读过的，不过忘掉复杂度了。忘掉具体做法了。。。只记得一张图，貌似分了上下两半先后做完再处理。。。不知道记得对不对。。。

 

于是乎懒得再看了。。自己YY了一个算法。

 

总的思想就是把原来印象中的栈改成一个队列，每次加入一个半平面的时候同时更新队头队尾。

 

最后剩下就是约束力最强的若干半平面了、、

 

排序嘛仍然要分类的。。类似上下两个部分。不过把一类的优先级设定成高于另一个然后直接排序即可。这里我选择的关键字是斜率。。。这样必然牵涉到了斜率为0的问题。。。在我的做法里直接赋成极大值。反正已经对这些半平面分成两大类了、、

 

排序完了就是处理的顺序了~~如下图：

 

（这里用射线的方向代表半平面的约束。沿射线方向的右边是约束的区域）

若是分两部分做的话两条绿边是原本上下部分的起始半平面，红色的是约束力弱的半平面被淘汰了。

 

关键就是怎么处理队头了、、

 

首先判类别。如果是同一类是没有能力影响到队头半平面的、、

然后是斜率。容易想象斜率比队头斜率大的第二大类的半平面是没有能力影响队头半平面的、、

（以上是自己YY的、、望指教）

 

还要注意一点，就是判断队头半平面的约束力是否强于当前处理的半平面，是的话忽略，否则入队、

 

处理队尾的操作就是原来的基础操作了、

 

然后。。。

 

就差不多了。。。

 

总结一下：没有分两部分处理的做法在任意的时刻，队内半平面都是目前处理到的约束力最强的半平面。更新的时候分两步，更新队尾半平面（卷边），更新队头半平面（封口）。

 

以上就是我YY出来的半平面非主流做法。。。不充分的地方以后有机会会继续更新、今天太迟了。。。

 

因为想法还不够完善吧代码量仍然巨大。。。

 

不过经过poj三题的洗礼，相信这种算法的正确还是有一定保证的、、

poj 1279：http://blog.sina.com.cn/s/blog_7c4c33190101aent.html

poj 3335：http://blog.sina.com.cn/s/blog_7c4c33190101aeo3.html

poj 1474：http://blog.sina.com.cn/s/blog_7c4c33190101aeo9.html

 

就是之前的三篇博文。

 

见笑、

 

poj 1279：

 

program pku_1279;
const eps=1e-7;
type dotsty=record
              x,y:double;
            end;
     linesty=record
               x1,y1,x2,y2,k:double;
               sty:longint;
             end;
var line:array[1..1500] of linesty;
    dot:array[1..1500] of dotsty;
    q:array[1..1500] of longint;
    n,cases:longint;
//============================================================================
procedure swap(x,y:longint);
var tt:linesty;
begin
  tt:=line[x]; line[x]:=line[y]; line[y]:=tt;
end;
//============================================================================
procedure qsort(l,r:longint);
var k:double;
    i,j:longint;
begin
  k:=line[(l+r) shr 1].k; i:=l; j:=r;
  repeat
    while line[i].k>k do inc(i);
    while line[j].k
    if i<=j then
    begin
      swap(i,j);
      inc(i); dec(j);
    end;
  until i>j;
  if l
  if i
end;
//============================================================================
procedure init;
var i,j:longint;
begin
  readln(n);
  for i:=1 to n do
    readln(line[i].x1,line[i].y1);
  for i:=1 to n-1 do line[i].x2:=line[i+1].x1;
  for i:=1 to n-1 do line[i].y2:=line[i+1].y1;
  line[n].x2:=line[1].x1; line[n].y2:=line[1].y1;
  for i:=1 to n do
  begin
    if (line[i].x1
      ((line[i].x1=line[i].x2) and (line[i].y1
        line[i].sty:=1 else line[i].sty:=2;
    if line[i].x1<>line[i].x2 then
      line[i].k:=(line[i].y1-line[i].y2)/(line[i].x1-line[i].x2) else
        line[i].k:=maxlongint;
  end; i:=1; j:=n;
  repeat
   while line[i].sty=1 do inc(i);
   while line[j].sty=2 do dec(j);
   if i<=j then swap(i,j);
  until i>j; qsort(1,j); qsort(i,n);
end;
//============================================================================
function left(x,y:double; l:linesty):boolean;
begin
  if (l.x2-l.x1)*(y-l.y1)-(x-l.x1)*(l.y2-l.y1)>-eps then
    exit(true) else exit(false);
end;
//============================================================================
function cross(l1,l2:linesty):dotsty;
var s1,s2,ss,k,dx,dy:double;
begin
  s1:=(l2.x2-l1.x1)*(l1.y2-l1.y1)-(l1.x2-l1.x1)*(l2.y2-l1.y1);
  s2:=(l1.x2-l1.x1)*(l2.y1-l1.y1)-(l2.x1-l1.x1)*(l1.y2-l1.y1);
  ss:=s1+s2;
  if abs(ss)
  begin
    cross.x:=maxlongint; exit;
  end;
  k:=s1/ss;
  dx:=l2.x1-l2.x2; dy:=l2.y1-l2.y2;
  cross.x:=l2.x2+dx*k;
  cross.y:=l2.y2+dy*k;
end;
//============================================================================
function half_plane:double;
var tmp:dotsty;
    tt:double;
    i,be,en:longint;
    flag:boolean;
begin
  be:=1; en:=1; q[1]:=1;
  for i:=2 to n do
  begin
    if abs(line[q[en]].k-line[i].k)
      if left(line[q[en]].x1,line[q[en]].y1,line[i]) then
        dec(en) else continue;
    while be
    begin
      tmp:=cross(line[i],line[q[en]]);
      if tmp.x=maxlongint then exit(0);
      if left(tmp.x,tmp.y,line[q[en-1]]) then
        dec(en) else break;
    end;
    while be
    begin
      if line[i].sty=1 then break;
      if line[i].k>line[q[be]].k-eps then break;
      tmp:=cross(line[i],line[q[be]]);
      if tmp.x=maxlongint then exit(0);
      if left(tmp.x,tmp.y,line[q[be+1]]) then
        inc(be) else break;
    end;
    if be
    begin
      tmp:=cross(line[i],line[q[en]]);
      if tmp.x=maxlongint then exit(0);
      if left(tmp.x,tmp.y,line[q[be]]) then continue;
    end; inc(en); q[en]:=i;
  end; tt:=0;
  if en
  for i:=be to en-1 do
    dot[i]:=cross(line[q[i]],line[q[i+1]]);
  dot[en]:=cross(line[q[en]],line[q[be]]);
  for i:=be to en-1 do
    tt:=tt+dot[i].x*dot[i+1].y-dot[i+1].x*dot[i].y;
  tt:=tt+dot[en].x*dot[be].y-dot[be].x*dot[en].y;
  tt:=abs(tt/2); exit(tt);
end;
//============================================================================
procedure work;
var s:double;
begin
  init;
  s:=half_plane;
  writeln(s:0:2);
end;
//============================================================================
begin
  assign(input,'1.in'); reset(input);
  readln(cases);
  for cases:=1 to cases do work;
end.