[PKU 1127 1410 USTC 1121]判断线段相交 判断点在多边形内

{

Bloxorz的程序突然出现了bug

调试得郁闷

不过这也算是好事吧

先做了几道计算几何

学了2个基本的计算几何方法

}

 

计算几何有一个基本的问题

就是判断线段相交

这里介绍2种方法

一种是“外积”法 另一种是求交点

先讲"外积"法

这里的叉积好像有点不伦不类 说是外积吧 求出来却是一个纯量

说不是吧 又好像有外积求模的影子 似乎是外积的简化吧

不扯了 我就用*号代替×号好了

(本来计算机领域也是用*的 又能不混淆 两全其美)

这种运算针对2个向量(xp,yp) (xq,yq)

p*q=(xp,yp)*(xq,yq)=xp·yq-xq·yp 得到一个纯量A

这个P很能反映2个向量的关系

我们归纳出以下几个性质

　　p*q=-q*p 即*运算有方向性 q乘到p是不同与p乘到q的

　　A<0 则向量p在q的逆时针方向

　　A=0 则向量p与q共线

　　A>0 则向量p在q的顺时针方向

这个性质让我们可以很好的判断线段和线段的位置关系

从而可以鉴别线段是否有跨立

我们认识到如果每条线段的2个端点分别在另一条线段的两侧

显然2条线段就相交了

我们称一条线段的2个端点分别在另一条线段的两侧称跨立

这就好比一个人两脚岔开 跨在一条白线两侧一样 很形象

那么 如何检测跨立呢?

我们用求*积解决

　　譬如判断l1对l2是否跨立

　　我们可以利用已知线段的端点 构造3个向量p q r

　　显然我们只需要判断p q是否分别在r两测即可

　　即判断(p*r)*(r*q)是否是正的 如果为正就可以判定跨立

　　如果为负就可以判定不跨立

　　一旦两次都判断跨立 则线段相交

讲到这里 我们还遗留有问题

如果线段的公共点是端点怎么办? 如果线段共线怎么办?

其实这就是(p*r)*(r*q)=0的情况 我们进行跨立试验的不足之处

我们还要引入快速排斥试验

　　正如其名 快速排斥好理解的多

　　只要判定分别以两线段为对角线的矩形是否相交

　　如果不相交必然线段也不相交

　　如果相交 则继续进行跨立试验

　　如果求出值为零 由于已经通过快速排斥试验 必然是有公共点的

　　也就可以直接判断了

判断矩形相交只需满足下图的条件即可

　　

这样我们完成了"外积"法判断的全过程

给出代码

 

var    xnp,xmq,xnq,xmp,ynp,ymq,ynq,ymp,
    xp1,yp1,xp2,yp2,xq1,yq1,xq2,yq2:longint;
    i,n:longint;
 function ifxy:boolean;
begin
if xp1<xp2 then begin xnp:=xp1; xmp:=xp2 end else begin xnp:=xp2; xmp:=xp1; end;
if yp1<yp2 then begin ynp:=yp1; ymp:=yp2 end else begin ynp:=yp2; ymp:=yp1; end;
if xq1<xq2 then begin xnq:=xq1; xmq:=xq2 end else begin xnq:=xq2; xmq:=xq1; end;
if yq1<yq2 then begin ynq:=yq1; ymq:=yq2 end else begin ynq:=yq2; ymq:=yq1; end;
if (xnp>xmq)or(xnq>xmp)or(ynp>ymq)or(ynq>ymp) then exit(false);
ifxy:=(((xp1-xq1)*(yq2-yq1)-(xq2-xq1)*(yp1-yq1))*((xq2-xq1)*(yp2-yq1)-(xp2-xq1)*(yq2-yq1))>=0)
and(((xq1-xp1)*(yp2-yp1)-(xp2-xp1)*(yq1-yp1))*((xp2-xp1)*(yq2-yp1)-(xq2-xp1)*(yp2-yp1))>=0);
end;
begin
assign(input,'seg.in'); reset(input);
assign(output,'seg.out'); rewrite(output);
readln(n);
for i:=1 to n do
    begin
    readln(xp1,yp1,xp2,yp2,xq1,yq1,xq2,yq2);
    writeln(ifxy);
    end;
close(input); close(output);
end.

还需要补充的一点是

 

我这里为了节约代码量 把判断正负的工作缩减为乘法

不过有时候会卡乘法 用乘法int64会爆界 要改成比较大小

 

接下来是交点法

名副其实 就是求交点 很好理解

得用一点解析几何的结论

比如 过2点 (x1,y1) (x2,y2)的直线

可以表示为 (y2-y1)x+(x1-x2)y+(y1-y2)x1+(x2-x1)y1=0

两条直线联立求2元一次方程组解即可

注意讨论有无解 是不是无数多解

还需注意除法的精度问题

 

var    x1,x2,x3,x4,y1,y2,y3,y4,i,n,
    a1,a2,b1,b2,c1,c2:longint;
    x0,y0:real;
 function max(a,b:longint):longint;
 begin
 if a>b then max:=a else max:=b;
 end;
 function min(a,b:longint):longint;
 begin
if a<b then min:=a else min:=b;
end;
function d:boolean;
begin
a1:=y2-y1; b1:=x1-x2; c1:=(y1-y2)*x1+(x2-x1)*y1;
a2:=y4-y3; b2:=x3-x4; c2:=(y3-y4)*x3+(x4-x3)*y3;
if (a1*b2=a2*b1)
    then if (a1*c2=a2*c1)
        then exit(true)
        else exit(false)
    else begin
    x0:=-(b2*c1-b1*c2)/(a1*b2-a2*b1);
    y0:=-(a2*c1-a1*c2)/(a2*b1-a1*b2);
    if (min(x1,x2)<=x0)and(max(x1,x2)>=x0)and(min(x3,x4)<=x0)and(max(x3,x4)>=x0)
    and(min(y1,y2)<=y0)and(max(y1,y2)>=y0)and(min(y3,y4)<=y0)and(max(y3,y4)>=y0)
        then exit(true) else exit(false);
    end;
end;
begin
assign(input,'seg.in'); reset(input);
assign(output,'seg0.out'); rewrite(output);
readln(n);
for i:=1 to n do
    begin
    readln(x1,y1,x2,y2,x3,y3,x4,y4);
    writeln(d);
    end;
close(input); close(output);
end.

p.s.

 

有人说用了除法会慢 我测过数据了 没这个事

200w个询问 2种方法都是16s左右 耗时相当

不过第一种会爆界 第二种精度会不够 也算是各有千秋

要选择去用

我用这个基本方法解决了2个问题 都很裸

一个是pku 1127 需要+并查集 or Floyd 不难

另一个是pku 1410 注意矩形包含内部 不仅仅是边界

pku1127的代码

 

const    maxn=13;
    maxd=maxlongint shr 1;
 var    xnp,xmp,xnq,xmq,ynp,ymp,ynq,ymq,i,n,x,y,k,j:longint;
    x1,x2,y1,y2:array[1..maxn]of longint;
    f:array[1..maxn,1..maxn]of longint;
 function ifxy(i,j:longint):boolean;
 begin
 if x1[i]<x2[i] then begin xnp:=x1[i]; xmp:=x2[i]; end else begin xnp:=x2[i]; xmp:=x1[i]; end;
 if x1[j]<x2[j] then begin xnq:=x1[j]; xmq:=x2[j]; end else begin xnq:=x2[j]; xmq:=x1[j]; end;
 if y1[i]<y2[i] then begin ynp:=y1[i]; ymp:=y2[i]; end else begin ynp:=y2[i]; ymp:=y1[i]; end;
 if y1[j]<y2[j] then begin ynq:=y1[j]; ymq:=y2[j]; end else begin ynq:=y2[j]; ymq:=y1[j]; end;
 if (xnp>xmq)or(xnq>xmp)or(ynp>ymq)or(ynq>ymp) then exit(false);
ifxy:=(((x1[i]-x1[j])*(y2[j]-y1[j])-(x2[j]-x1[j])*(y1[i]-y1[j]))*((x2[j]-x1[j])*(y2[i]-y1[j])-(x2[i]-x1[j])*(y2[j]-y1[j]))>=0)
 and(((x1[j]-x1[i])*(y2[i]-y1[i])-(x2[i]-x1[i])*(y1[j]-y1[i]))*((x2[i]-x1[i])*(y2[j]-y1[i])-(x2[j]-x1[i])*(y2[i]-y1[i]))>=0)
 end;
 begin
assign(input,'jack.in'); reset(input);
assign(output,'jack.out'); rewrite(output);
readln(n);
 while n<>0 do
    begin
    for i:=1 to n do readln(x1[i],y1[i],x2[i],y2[i]);
    for i:=1 to n do for j:=1 to n do if  ifxy(i,j) then f[i,j]:=1 else f[i,j]:=maxd;
    for i:=1 to n do f[i,i]:=0;
    for k:=1 to n do
        for i:=1 to n do
            for j:=1 to n do
                if f[i,j]>f[i,k]+f[k,j] then f[i,j]:=f[i,k]+f[k,j];
    readln(x,y);
    while x<>0 do
        begin
        if f[x,y]<maxlongint shr 1
            then writeln('CONNECTED')
            else writeln('NOT CONNECTED');
        readln(x,y);
        end;
    readln(n);
    end;
close(input); close(output);
 end.
 

pku1410的代码

 

 

var    xnp,xmp,ynp,ymp,xnq,xmq,ynq,ymq,i,n,x1,x2,y1,y2,xp1,xp2,yp1,yp2:longint;
 function ifxy(xq1,yq1,xq2,yq2:longint):boolean;
 begin
 if xp1<xp2 then begin xnp:=xp1; xmp:=xp2 end else begin xnp:=xp2; xmp:=xp1; end;
if yp1<yp2 then begin ynp:=yp1; ymp:=yp2 end else begin ynp:=yp2; ymp:=yp1; end;
if xq1<xq2 then begin xnq:=xq1; xmq:=xq2 end else begin xnq:=xq2; xmq:=xq1; end;
if yq1<yq2 then begin ynq:=yq1; ymq:=yq2 end else begin ynq:=yq2; ymq:=yq1; end;
if (xnp>xmq)or(xnq>xmp)or(ynp>ymq)or(ynq>ymp) then exit(false);
ifxy:=(((xp1-xq1)*(yq2-yq1)-(xq2-xq1)*(yp1-yq1))*((xq2-xq1)*(yp2-yq1)-(xp2-xq1)*(yq2-yq1))>=0)
and(((xq1-xp1)*(yp2-yp1)-(xp2-xp1)*(yq1-yp1))*((xp2-xp1)*(yq2-yp1)-(xq2-xp1)*(yp2-yp1))>=0);
end;
begin
assign(input,'inter.in'); reset(input);
assign(output,'inter.out'); rewrite(output);
readln(n);
for i:=1 to n do
    begin
    readln(xp1,yp1,xp2,yp2,x1,y1,x2,y2);
    if x1>x2 then begin xnp:=x1; x1:=x2; x2:=xnp; end;
    if y1>y2 then begin xnp:=y1; y1:=y2; y2:=xnp; end;
    if     (x1<xp1)and(xp1<x2)and(x1<xp2)and(xp2<x2)and(y1<yp1)and(yp1<y2)and(y1<yp2)and(yp2<y2)or
        (ifxy(x1,y1,x1,y2))or(ifxy(x1,y1,x2,y1))or(ifxy(x1,y2,x2,y2))or(ifxy(x2,y1,x2,y2))
        then writeln('T') else writeln('F');
    end;
close(input); close(output);
end.

 

 

 

 

还有一个方法是判断点在多边形内部

一种方法是射线法 需要用到上面的判断线段相交

并且就是我说的爆界爆得利害的地方...

我们由点向外射一条射线 由于多边形是有界的

必然一部分在外边 我们沿着射线向内走

碰到一次边界说明我们从多边形外到了多边形内 或 从多边形内到了多边形外

显然这样改变的次数为奇数时 我们到达那个点时是在多边形内的即点在形内

这只是一个粗略的想法

还有很多细节需要注意

　　这里列出来

　　取射线可以取 垂直于坐标轴的射线 方便处理

　　射线其实就是另一端点很远的线段

　　　　这里由于线段的端点很大 很有可能在判断相交时爆界要注意

　　如果点在多边形边上直接特判

　　碰到与射线重合的边不予考虑

　　碰到与在射线上的点 如果是那条边上横坐标较大的点就计数 否则忽略

这样的需要扫描所有的边 复杂度是O(n)

还有很多方法譬如面积法 有兴趣可以去baidu

给出代码

 

const    maxn=50;
    inf=100000000;
var    n,m,i,ans:longint;
    x,y:array[1..maxn+1]of longint;
    xnp,xmq,xnq,xmp,ynp,ymq,ynq,ymp,x0,y0:int64;
function ifxy(xp1,yp1,xp2,yp2,xq1,yq1,xq2,yq2:int64):boolean;
begin
if xp1<xp2 then begin xnp:=xp1; xmp:=xp2 end else begin xnp:=xp2; xmp:=xp1; end;
if yp1<yp2 then begin ynp:=yp1; ymp:=yp2 end else begin ynp:=yp2; ymp:=yp1; end;
if xq1<xq2 then begin xnq:=xq1; xmq:=xq2 end else begin xnq:=xq2; xmq:=xq1; end;
if yq1<yq2 then begin ynq:=yq1; ymq:=yq2 end else begin ynq:=yq2; ymq:=yq1; end;
if (xnp>xmq)or(xnq>xmp)or(ynp>ymq)or(ynq>ymp) then exit(false);
ifxy:=(((xp1-xq1)*(yq2-yq1)-(xq2-xq1)*(yp1-yq1)>=0)and((xq2-xq1)*(yp2-yq1)-(xp2-xq1)*(yq2-yq1)>=0)
or((xp1-xq1)*(yq2-yq1)-(xq2-xq1)*(yp1-yq1)<=0)and((xq2-xq1)*(yp2-yq1)-(xp2-xq1)*(yq2-yq1)<=0))
and(((xq1-xp1)*(yp2-yp1)-(xp2-xp1)*(yq1-yp1)>=0)and((xp2-xp1)*(yq2-yp1)-(xq2-xp1)*(yp2-yp1)>=0)
or((xq1-xp1)*(yp2-yp1)-(xp2-xp1)*(yq1-yp1)<=0)and((xp2-xp1)*(yq2-yp1)-(xq2-xp1)*(yp2-yp1)<=0));
end;
function ifp:longint;
var    i,j:longint;
begin
j:=0;
for i:=1 to n do
    begin
    if ifxy(x0,y0,x0,y0,x[i],y[i],x[i+1],y[i+1]) then exit(true);
    if x[i]=x[i+1] then continue;
    if ifxy(x[i],y[i],x[i],y[i],x0,y0,x0,inf) then begin if x[i]>x[i+1] then inc(j); continue; end;
    if ifxy(x[i+1],y[i+1],x[i+1],y[i+1],x0,y0,x0,inf) then begin if x[i+1]>x[i] then inc(j); continue; end;
    if ifxy(x0,y0,x0,inf,x[i],y[i],x[i+1],y[i+1]) then inc(j);
    end;
if j and 1=1 then ifp:=true else ifp:=false;
end;
begin
assign(input,'check.in'); reset(input);
assign(output,'check.out'); rewrite(output);
while not eof do
    begin
    readln(n,m);
    for i:=1 to n do readln(x[i],y[i]);
    x[n+1]:=x[1]; y[n+1]:=y[1];
    for i:=1 to m do
        begin
        readln(x0,y0);
        writeln(ifb);
        end;
    end;
close(input); close(output);
end.

注意这里就对判断线段相交作了改动

 

我用这个方法解决了USTC 1121

题目也很裸 多边形内的点权为2 在外为1 求权和

我爆界爆了一个下午 哭...

 

const    maxn=50;
    inf=100000000;
var    n,m,i,ans:longint;
    x,y:array[1..maxn+1]of longint;
    xnp,xmq,xnq,xmp,ynp,ymq,ynq,ymp,x0,y0,x1,y1:int64;
function ifxy(xp1,yp1,xp2,yp2,xq1,yq1,xq2,yq2:int64):boolean;
begin
if xp1<xp2 then begin xnp:=xp1; xmp:=xp2 end else begin xnp:=xp2; xmp:=xp1; end;
if yp1<yp2 then begin ynp:=yp1; ymp:=yp2 end else begin ynp:=yp2; ymp:=yp1; end;
if xq1<xq2 then begin xnq:=xq1; xmq:=xq2 end else begin xnq:=xq2; xmq:=xq1; end;
if yq1<yq2 then begin ynq:=yq1; ymq:=yq2 end else begin ynq:=yq2; ymq:=yq1; end;
if (xnp>xmq)or(xnq>xmp)or(ynp>ymq)or(ynq>ymp) then exit(false);
ifxy:=(((xp1-xq1)*(yq2-yq1)-(xq2-xq1)*(yp1-yq1)>=0)and((xq2-xq1)*(yp2-yq1)-(xp2-xq1)*(yq2-yq1)>=0)
or((xp1-xq1)*(yq2-yq1)-(xq2-xq1)*(yp1-yq1)<=0)and((xq2-xq1)*(yp2-yq1)-(xp2-xq1)*(yq2-yq1)<=0))
and(((xq1-xp1)*(yp2-yp1)-(xp2-xp1)*(yq1-yp1)>=0)and((xp2-xp1)*(yq2-yp1)-(xq2-xp1)*(yp2-yp1)>=0)
or((xq1-xp1)*(yp2-yp1)-(xp2-xp1)*(yq1-yp1)<=0)and((xp2-xp1)*(yq2-yp1)-(xq2-xp1)*(yp2-yp1)<=0));
end;
function ifp:longint;
var    i,j:longint;
begin
j:=0;
for i:=1 to n do
    begin
    if ifxy(x0,y0,x0,y0,x[i],y[i],x[i+1],y[i+1]) then exit(1);
    if x[i]=x[i+1] then continue;
    if ifxy(x[i],y[i],x[i],y[i],x0,y0,x0,inf) then begin if x[i]>x[i+1] then inc(j); continue; end;
    if ifxy(x[i+1],y[i+1],x[i+1],y[i+1],x0,y0,x0,inf) then begin if x[i+1]>x[i] then inc(j); continue; end;
    if ifxy(x0,y0,x0,inf,x[i],y[i],x[i+1],y[i+1])
        then inc(j);
    end;
ifp:=j and 1;
end;
begin
assign(input,'island.in'); reset(input);
assign(output,'island.out'); rewrite(output);
readln(n,m);
while n<>0 do
    begin
    for i:=1 to n do
    readln(x[i],y[i]);
    x[n+1]:=x[1]; y[n+1]:=y[1];
    ans:=0;
    for i:=1 to m do
        begin
        readln(x0,y0);
        inc(ans); ans:=ans+ifp;
        end;
    writeln(ans);
    readln(n,m);
    end;
close(input); close(output);
end.

计算几何代码都一坨一坨的

 

要有耐心 要注意常数 精度 爆界...