Question

求 半平面交 的程序 Pascal/C/C++ 均可

如题
时间复杂度在o(N^2）或以下

Answer 1

这个是N^2的：

Program Pyf;
//5294087 piyifan 2451 Accepted 1600K 532MS Pascal 2165B 2009-06-13 20:21:20
const
zero = 1e-8;
inf = 10000;
type
pt = record
x , y : double;
end;
Tline = record
a , b : pt;
end;
var
list , L : array[0..20000] of pt;
line : array[0..20000] of Tline;
n : longint;

function det( a , b , c : pt ) : double;
begin
exit( ( b.x - a.x ) * ( c.y - a.y ) - ( b.y - a.y ) * ( c.x - a.x ) );
end; { det }

function dis(a , b : pt ) : double;
begin
exit( sqrt( sqr( a.x - b.x ) + sqr( a.y - b.y ) ) );
end; { dis }

procedure getjiao( a , b , c , d : pt; var p : pt );
var
u , v : double;
begin
u := abs( det( a , c , b ) ); v := abs( det( a , b , d ) );
p.x := ( c.x * v + d.x * u ) / ( u + v );
p.y := ( c.y * v + d.y * u ) / ( u + v );
end; { getjiao }

function dblcmp( a : double ) : longint;
begin
if abs( a ) < zero then exit( 0 );
if a > 0 then exit( 1 )
else exit( -1 );
end; { dblcmp }

procedure work;
var
i , j , tot , t : longint;
ans : double;
begin
list[1].x := 0; list[1].y := 0;
list[2].x := inf; list[2].y := 0;
list[3].x := inf; list[3].y := inf;
list[4].x := 0; list[4].y := inf;
list[5] := list[1];
tot := 4;
for i := 1 to n do
begin
t := 0;
for j := 1 to tot do
begin
if dblcmp( det( line[i].a , line[i].b , list[j] ) ) >= 0 then
begin
inc( t );
L[t] := list[j];
end;
if dblcmp( det( line[i].a , line[i].b , list[j] ) )
xor dblcmp( det( line[i].a , line[i].b , list[j + 1] ) ) = -2 then
begin
inc( t );
getjiao( line[i].a , line[i].b , list[j] , list[j + 1] , L[t] );
end;
end;
tot := t;
for j := 1 to tot do list[j] := L[j];
list[tot + 1] := list[1];
if tot < 3 then tot := 0;
end;
ans := 0;
for i := 2 to tot do
ans := ans + det( list[1] , list[i - 1] , list[i] );
writeln( ans / 2 : 0 : 1 );
end; { work }

procedure init;
var
i : longint;
begin
readln( n );
for i := 1 to n do
readln( line[i].a.x , line[i].a.y , line[i].b.x , line[i].b.y );
end; { init }

begin
init;
work;
end.

这个是NlogN的：
program xlmj;
{$inline on}
{$R-,Q-,S-}
const
maxn = 54;
size = 10;
zero = 1e-6;

type
pointT = record
x,y : extended;
end;
lineT = record
p,q : pointT;
arc,c : extended;
end;

var
line : array [0..maxn] of lineT;
qj : array [0..maxn] of pointT;
q : array [0..maxn] of longint;
ini : array [0..maxn] of pointT;
m,n,i : longint;
ans : extended;

operator -(var a,b : pointT )c: pointT; inline;
begin
c.x:=a.x-b.x; c.y:=a.y-b.y;
end;

function eq(var a,b : extended ): boolean; inline;
begin
eq:=abs(a-b)<zero;
end;

function cha(p1,p2 : pointT ): extended; inline;
begin
cha:=p1.x*p2.y-p2.x*p1.y;
end;

function carc(x,y : extended ): extended; inline;
begin
carc:=arctan(y/x);
if x<0 then carc:=carc+pi
else if y<0 then carc:=carc+pi*2;
end;

procedure setline(x0,y0,x1,y1 : extended; i : longint ); inline;
var
len,sina,cosa,x,y : extended;
begin
with line[i] do
begin
p.x:=x0; p.y:=y0; q.x:=x1; q.y:=y1;
if eq(x0,x1) then
begin
c:=x0;
if y1>y0 then arc:=pi/2 else arc:=pi/2*3;
end else
if eq(y0,y1) then
begin
c:=y0;
if x1>x0 then arc:=0 else arc:=pi;
end else
begin
arc:=carc(x1-x0,y1-y0);
c:=x0-y0*(x1-x0)/(y1-y0);
end;
end;
end;

procedure init; inline;
var
k,t : extended;
p1,p2,p : pointT;
ch : char;
begin
inc(m); readln(ini[m].x,ini[m].y,ch,ch);
if eq(ini[m].x,ini[m-1].x) and eq(ini[m].y,ini[m-1].y) then exit;
if ch='S' then begin ans:=0; exit; end;
if ch='H' then t:=-1 else t:=1;
if eq(ini[m].y,ini[m-1].y) then
begin
p1.x:=(ini[m].x+ini[m-1].x)/2; p2.x:=p1.x;
p1.y:=0; p2.y:=1;
end else
begin
if eq(ini[m].x,ini[m-1].x) then k:=0
else k:=(ini[m-1].x-ini[m].x)/(ini[m].y-ini[m-1].y);
p1.x:=(ini[m].x+ini[m-1].x)/2; p2.x:=p1.x+1;
p1.y:=(ini[m].y+ini[m-1].y)/2; p2.y:=p1.y+k;
end;
p.x:=ini[m].x; p.y:=ini[m].y;
if cha(p2-p1,p-p1)*t>=0 then begin p:=p1; p1:=p2; p2:=p; end;
inc(n); setline(p1.x,p1.y,p2.x,p2.y,n);
// writeln(p1.x:0:0,' ',p1.y:0:0,' ',p2.x:0:0,' ',p2.y:0:0);
end;

function pre(var a,b : lineT ): boolean; inline;
begin
if a.arc+zero<b.arc then exit(true);
if b.arc+zero<a.arc then exit(false);
if (0<a.arc) and (a.arc<=pi) then pre:=a.c<b.c else pre:=a.c>b.c;
end;

procedure qsort(i,j : longint ); inline;
var
s,t : longint;
m,a : lineT;
begin
if i>=j then exit;
s:=i; t:=j; m:=line[(s+t) shr 1];
repeat
while pre(line[s],m) do inc(s);
while pre(m,line[t]) do dec(t);
if s<=t then
begin
a:=line[s]; line[s]:=line[t]; line[t]:=a;
inc(s); dec(t);
end;
until s>t;
qsort(i,t); qsort(s,j);
end;

procedure getline(var x1,y1,x2,y2 : extended; var a,b,c : extended ); inline;
begin
a:=y1-y2; b:=x2-x1;
c:=x1*y2-x2*y1;
end;

function getjiao(var p1,p2,p3,p4,p5 : pointT ): longint; inline;
var
a1,b1,c1,a2,b2,c2 : extended;
begin
getline(p1.x,p1.y,p2.x,p2.y,a1,b1,c1);
getline(p3.x,p3.y,p4.x,p4.y,a2,b2,c2);
p5.x:=(b1*c2-b2*c1)/(a1*b2-a2*b1);
p5.y:=(c1*a2-c2*a1)/(a1*b2-a2*b1);
end;

procedure halfS_cross; inline;
var
s,t,j,x,y : longint;
begin
if abs(ans)<zero then begin writeln('0.00'); exit; end;
qsort(1,n); i:=n; n:=0;
for i:=1 to i do
if (i=1) or (not eq(line[i].arc,line[i-1].arc)) then
begin
inc(n); line[n]:=line[i];
end;

s:=1; t:=1; q[1]:=1; qj[1]:=line[1].p;
qj[1].x:=qj[1].x-(line[1].q.x-line[1].p.x)*size;
for i:=2 to n do qj[i].x:=size;

for i:=2 to n do
begin
if eq(qj[t+1].x,qj[s].x) and eq(qj[t+1].y,qj[s].y) and
(cha(qj[s]-line[i].p,qj[s]-line[i].q)>=zero) then continue;
while (s<=t) and (cha(qj[t]-line[i].p,qj[t]-line[i].q)<-zero) do dec(t);
if s>t then begin ans:=0; writeln('0.00'); exit; end;
getjiao(line[q[t]].p,line[q[t]].q,line[i].p,line[i].q,qj[t+1]);
inc(t); q[t]:=i;

if line[i].p.y>line[i].q.y then
begin
while (s<=t) and (cha(qj[s+1]-line[i].p,qj[s+1]-line[i].q)<-zero) do inc(s);
if s>t then begin ans:=0; writeln('0.00'); exit; end;
getjiao(line[q[s]].p,line[q[s]].q,line[i].p,line[i].q,qj[s]);
qj[t+1]:=qj[s];
end;
end;
qj[t+1]:=qj[s]; ans:=0;
for i:=s+1 to t-1 do
ans:=ans+cha(qj[i]-qj[s],qj[i+1]-qj[s]);
writeln(ans/2:0:2);
end;

begin
ans:=100; n:=4;
setline(0,0,size,0,1);
setline(size,0,size,size,2);
setline(size,size,0,size,3);
setline(0,size,0,0,4);
while not seekeof do
begin
init;
halfS_cross;
end;
end.

Answer 2

第一次写半平面交 居然写了6个多小时..

而且由于这题的特殊性(是把一个凸多边形的边往里面移出来的半平面),该代码很可能还有问题 不过先贴上来

还有就是 该程序没有添加外围边框 所以要小心

////////////////////////////

该程序又过了1279 不过据说这题数据很弱..在prework里面已经添加了外围边框 不过这会改变n所以如果后面用到初始点数 小心点数被改了

/////////////////

该程序过了n题 看来应该是对的了 注意输入点必须顺时针噢 顺便加上计算面积的模块

////////////////

poj 1755 精度要求非常高 判断面积为0要做到10^(-16) 所以上限maxl也要加大

并且为了统一 更改半平面的极角f为其正向法向量极角 注意到ax+by+c>=0的正向法向量恰好是(a,b)

#include <iostream>
#include <cmath>
using namespace std;
const int maxn=210;
const double pi=3.1415926535897932384626433832795;
struct point//点
{
double x,y;
};
struct line//半平面直线 表示为a*x+b*y+c>=0
{
double a,b,c;
double f;
double func(point s)
{
return a*s.x+b*s.y+c;
}
};
point acro(line p,line q)//无视平行情况计算两条直线交点
{
point s;
// if (p.a*q.b-p.b*q.a==0) return;
s.y=(p.a*q.c-p.c*q.a)/(q.a*p.b-p.a*q.b);
s.x=(p.b*q.c-p.c*q.b)/(p.a*q.b-p.b*q.a);
return s;
}
line que[maxn];
int head,tail;
point a[maxn];
line b[maxn];
int n;
int needtopop(line s,line a,line b)//是否需要退栈
{
if (s.func(acro(a,b))<0) return 1;
return 0;
}
int addone(line s)//增加一个点
{
int hhead=head,ttail=tail;
while (head<tail&&needtopop(s,que[tail],que[tail-1])) tail--;
while (head<tail&&needtopop(s,que[head],que[head+1])) head++;
tail++;
que[tail]=s;
return 1;
}
line linemove(line s,double d)//半平面正向平移d
{
line c=s;
c.c-=d;
return c;
}
int comp(const void *a,const void *b)//排序
{
line *p=(line*)a,*q=(line*)b;
if (p->f>q->f) return 1;
if (p->f<q->f) return -1;
if (p->c>q->c) return 1;
if (p->c<q->c) return -1;
return 0;
}
int banpingmianjiao(double d)//主函数
{
line c[maxn];
int i,j;
for (i=0;i<n;i++) c[i]=linemove(b[i],d);//c储存要做的半平面
qsort(c,n,sizeof(line),comp);//极角排序
j=0;
for (i=1;i<n;i++)
if (c[i-1].f==c[i].f) j++;
else c[i-j]=c[i];//干掉极角相同的
c[n-j]=c[0];
int m=n-j;
if (m<=1) return 1;
head=0;tail=1;
que[0]=c[0];que[1]=c[1];

for (i=2;i<m;i++)
if (addone(c[i])==0) //加线进去
return 0;
while (head<tail&&needtopop(que[tail],que[head],que[head+1])) head++;//删掉头多余线
while (head<tail&&needtopop(que[head],que[tail],que[tail-1])) tail--;//删掉尾巴多余线
if (tail-head<2) return 0;// 必须有3条边 否则挂
return 1;
}

line poitoline(point p,point q)//将点q-〉p专成半平面 直线和法向量成右手螺旋
{
line c;
c.a=q.y-p.y;
c.b=p.x-q.x;
c.c=p.y*q.x-p.x*q.y;
double one=sqrt(c.a*c.a+c.b*c.b);
c.a/=one;c.b/=one;c.c/=one;
c.f=atan2(c.b,c.a);
return c;
}
void prework() //预处理将凸包（顺时针）转出半平面来
{
int i;
a[n]=a[0];
for (i=0;i<n;i++)
b[i]=poitoline(a[i],a[i+1]);

const double maxl=999999999;
point wei[5];
wei[0].x=-maxl;wei[0].y=-maxl;
wei[1].x=-maxl;wei[1].y=maxl;
wei[2].x=maxl;wei[2].y=maxl;
wei[3].x=maxl;wei[3].y=-maxl;
wei[4]=wei[0];
for (i=0;i<4;i++) b[n+i]=poitoline(wei[i],wei[i+1]);
n+=4;
}

double countarea(point* s,int m)
{
int i;
double tot=0;
s[m]=s[0];
for (i=0;i<m;i++) tot+=s[i].x*s[i+1].y-s[i].y*s[i+1].x;
if (tot<0) return -tot/2;
return tot/2;
}
void work()//二分
{
double l,r;
l=0;r=10000000;
while (l+0.000001<r)
{
double mid=(l+r)/2;
if (banpingmianjiao(mid)) l=mid;else r=mid;
}
printf("%.6lf\n",l);
}
int main()
{
// freopen("in.in","r",stdin);
while (scanf("%d",&n),n)
{
int i;
for (i=n-1;i>=0;i--) scanf("%lf%lf",&a[i].x,&a[i].y);
prework();
work();
}
return 0;
}

Answer 3

自己去看论文，我记得是朱泽园05或06年的，复杂度是NlogN

Answer 4

解释什么半平面交？