交大计算几何模板Word文件下载.docx
《交大计算几何模板Word文件下载.docx》由会员分享，可在线阅读，更多相关《交大计算几何模板Word文件下载.docx（51页珍藏版）》请在冰豆网上搜索。
%lf%lf"
&
x,&
y）;
pointoperator-（pointa）{a.x=x-a.x;
a.y=y-a.y;
returna;
pointoperator+（pointa）{a.x=x+a.x;
a.y=y+a.y;
pointoperator/（doublea）{returnpoint（x/a,y/a）;
booloperator==（constpoint&
b）{return!
cmp（x-b.x）&
&
!
cmp（y-b.y）;
};
doublearea（pointa,pointb,pointc）{
return（b.x-a.x）*（c.y-a.y）-（b.y-a.y）*（c.x-a.x）;
doubledot（pointa,pointb,pointc）{
return（b-a）^（c-a）;
doubledis（pointa）{returnsqrt（a.x*a.x+a.y*a.y）;
doubledis（pointa,pointb）{returndis（b-a）;
//================两点求线
linepoint_make_line（pointa,pointb）{
lineh;
h.a=b.y-a.y;
h.b=-（b.x-a.x）;
h.c=-a.x*b.y+a.y*b.x;
returnh;
//===========旋转角度p的向量
pointrotate_point（pointa,doublep）{
pointh;
h.x=a.x*cos（p）-a.y*sin（p）;
h.y=a.x*sin（p）+a.y*cos（p）;
//================点P到线段st的距离======================
doubledis_point_segment（pointp,points,pointt）{
if（（（p-t）^（s-t））>
0&
（（（p-s）^（t-s））>
0））returnfabs（（p-s）*（t-s））/dis（s-t）;
elsereturnmin（dis（p-s）,dis（p-t））;
//========一个点关于直线作镜像
voidPointProjLine（constpoint&
p0,constpoint&
p1,constpoint&
p2,point&
cp）{
doublet=dot（p1,p2,p0）/dot（p1,p2,p2）;
cp.x=p1.x+t*（p2.x-p1.x）;
cp.y=p1.y+t*（p2.y-p1.y）;
//===!
!
或者===
doublePointToLine（constpoint&
p0,constpoint&
p1,constpoint&
p2,point&
cp）
{
doubled=dis（p1,p2）;
doubles=area（p1,p2,p0）/d;
cp.x=p0.x+s*（p2.y-p1.y）/d;
cp.y=p0.y-s*（p2.x-p1.x）/d;
returns;
voidReflectPoint（constpoint&
pointp3;
PointToLine（p0,p1,p2,p3）;
//PointProjLine（p0,p1,p2,p3）;
都是求影射点
cp=p3+（p3-p0）;
//=========判点是否在线段上
boolPointOnSegment（pointp,points,pointt）{
returncmp（area（p,s,t））==0&
cmp（dot（p,s,t））<
=0;
//================两线交点
pointline_make_point（linea,lineb）{
h.y=-（a.c*b.a-b.c*a.a）/（a.b*b.a-b.b*a.a）;
//=====makesureaandbaren'
tparallel
if（abs（a.a）<
eps）h.x=（-b.c-h.y*b.b）/b.a;
elseh.x=（-a.c-h.y*a.b）/a.a;
//==========线段平移D的长度
linemove_d（linea,constdoubled）{
return（line）{a.a,a.b,a.c+d*sqrt（a.a*a.a+a.b*a.b）};
//==========判平行
boolparallel（linea,lineb）{
if（cmp（a.b*b.a-b.b*a.a）==0）returntrue;
returnfalse;
//========================点与多边形线段与多边形=================
intPointInPolygon（pointcp,pointa[],intn）{
inti,k,d1,d2,wn=0;
a[n]=a[0];
rep（i,n）{
if（PointOnSegment（cp,a[i],a[i+1]））return2;
k=cmp（area（a[i],a[i+1],cp））;
d1=cmp（a[i+0].y-cp.y）;
d2=cmp（a[i+1].y-cp.y）;
if（k>
0&
d1<
=0&
d2>
0）wn++;
if（k<
d2<
d1>
0）wn--;
}
returnwn!
//========================================判断线段是否有在多边形内部=========
boolcompareab（constpoint&
a,constpoint&
b）{
if（a.x<
b.x||（a.x==b.x&
a.y<
b.y））returntrue;
elsereturnfalse;
pointstack[11000];
boolSegmentCrossPolygon（points,pointt,pointa[],intn）{
inti,j,k,m1,m2,closed;
linee1,e2;
pointcross;
if（PointInPolygon（s,a,n）==1||PointInPolygon（t,a,n）==1）returntrue;
closed=1;
stack[closed]=s;
e1=point_make_line（s,t）;
rep（i,n）{
k=i+1;
e2=point_make_line（a[i],a[k]）;
if（!
parallel（e1,e2））{
cross=line_make_point（e1,e2）;
if（PointOnSegment（cross,s,t）&
PointOnSegment（cross,a[i],a[k]））{
closed++;
stack[closed]=cross;
}
stack[closed]=t;
sort（stack+1,stack+closed+1,compareab）;
foru（i,1,closed-1）{
cross=（stack[i]+stack[i+1]）/2;
if（PointInPolygon（cross,a,n）==1）
returntrue;
//多边形的重心
voidPolygonCentroids（pointp[],intn,point&
cp）{
//if面积为0需要特判
doublesum=0,s=0;
cp.x=0;
cp.y=0;
for（inti=1;
i<
n-1;
i++,sum+=s）{
s=area（p[0],p[i],p[i+1]）;
cp.x+=s*（p[0].x+p[i].x+p[i+1].x）;
cp.y+=s*（p[0].y+p[i].y+p[i+1].y）;
cp.x/=sum*3;
cp.y/=sum*3;
pointgravity（point*p,intn）{
//if面积为0需要特判
doublearea=0;
pointcenter;
center.x=0;
center.y=0;
p[n]=p[0];
for（inti=0;
i<
n-1;
i++）{
area+=（p[i].x*p[i+1].y-p[i+1].x*p[i].y）/2;
center.x+=（p[i].x*p[i+1].y-p[i+1].x*p[i].y）*（p[i].x+p[i+1].x）;
center.y+=（p[i].x*p[i+1].y-p[i+1].x*p[i].y）*（p[i].y+p[i+1].y）;
}
area+=（p[n-1].x*p[0].y-p[0].x*p[n-1].y）/2;
center.x+=（p[n-1].x*p[0].y-p[0].x*p[n-1].y）*（p[n-1].x+p[0].x）;
center.y+=（p[n-1].x*p[0].y-p[0].x*p[n-1].y）*（p[n-1].y+p[0].y）;
center.x/=6*area;
center.y/=6*area;
returncenter;
//===============================圆
doubleangle（pointp0,pointp1,pointp2）{
doublecr=area（p0,p1,p2）;
doubledt=dot（p0,p1,p2）;
if（cmp（cr）==0）cr=0.0;
if（cmp（dt）==0）dt=0.0;
returnatan2（cr,dt）;
//-pi~pi
voidCircleCenter（pointp0,pointp1,pointp2,point&
doublea1=p1.x-p0.x,b1=p1.y-p0.y,c1=（sqr（a1）+sqr（b1））/2;
doublea2=p2.x-p0.x,b2=p2.y-p0.y,c2=（sqr（a2）+sqr（b2））/2;
doubled=a1*b2-a2*b1;
cp.x=p0.x+（c1*b2-c2*b1）/d;
cp.y=p0.y+（a1*c2-a2*c1）/d;
//三角形内心INPUT:
（242,89）,（212,185）,（71,128）,OUTPUT:
（189.5286,137.4987）
doubleIncenter（pointA,pointB,pointC,point&
doubles,p,r,a,b,c;
a=dis（B,C）,b=dis（C,A）,c=dis（A,B）;
p=（a+b+c）/2;
s=sqrt（p*（p-a）*（p-b）*（p-c））;
r=s/p;
cp.x=（a*A.x+b*B.x+c*C.x）/（a+b+c）;
cp.y=（a*A.y+b*B.y+c*C.y）/（a+b+c）;
returnr;
//三角形外心INPUT:
（208.8229,171.0697）
voidOrthocenter（pointA,pointB,pointC,point&
CircleCenter（A,B,C,cp）;
cp.x=A.x+B.x+C.x-2*cp.x;
cp.y=A.y+B.y+C.y-2*cp.y;
//园外一点p0,半径为r，直线ax+by+c=0的交点
intCircleLine（pointp0,doubler,doublea,doubleb,doublec,point&
cp1,point&
cp2）{
doubleaa=a*a,bb=b*b,s=aa+bb;
doubled=r*r*s-sqr（a*p0.x+b*p0.y+c）;
if（d+eps<
0）return0;
if（d<
eps）d=0;
elsed=sqrt（d）;
doubleab=a*b,bd=b*d,ad=a*d;
doublexx=bb*p0.x-ab*p0.y-a*c;
doubleyy=aa*p0.y-ab*p0.x-b*c;
cp2.x=（xx+bd）/s;
cp2.y=（yy-ad）/s;
cp1.x=（xx-bd）/s;
cp1.y=（yy+ad）/s;
if（d>
eps）return2;
elsereturn1;
//两园交线CommonAxisof|P-P1|=r1and|P-P2|=r2oftheax+by+c=0form
voidCommonAxis（pointp1,doubler1,pointp2,doubler2,double&
a,double&
b,double&
c）{
doublesx=p2.x+p1.x,mx=p2.x-p1.x;
doublesy=p2.y+p1.y,my=p2.y-p1.y;
a=2*mx;
b=2*my;
c=-sx*mx-sy*my-（r1+r2）*（r1-r2）;
//两园交点Crossingof|P-P1|=r1and|P-P2|=r2
//两个圆不能共圆心，请特判
intCircleCrossCircle（pointp1,doubler1,pointp2,doubler2,point&
cp2）{
doublemx=p2.x-p1.x,sx=p2.x+p1.x,mx2=mx*mx;
doublemy=p2.y-p1.y,sy=p2.y+p1.y,my2=my*my;
doublesq=mx2+my2,d=-（sq-sqr（r1-r2））*（sq-sqr（r1+r2））;
if（d+eps<
0）return0;
if（d<
eps）d=0;
elsed=sqrt（d）;
doublex=mx*（（r1+r2）*（r1-r2）+mx*sx）+sx*my2;
doubley=my*（（r1+r2）*（r1-r2）+my*sy）+sy*mx2;
doubledx=mx*d,dy=my*d;
sq*=2;
cp1.x=（x-dy）/sq;
cp1.y=（y+dx）/sq;
cp2.x=（x+dy）/sq;
cp2.y=（y-dx）/sq;
if（d>
//====两园面积交dist=是距离dis是距离的平方
doubletwoCircleAreaUnion（pointa,pointb,doubler1,doubler2）{
if（r1+r2<
=（a-b）.dist（））return0;
if（r1+（a-b）.dist（）<
=r2）returnpi*r1*r1;
if（r2+（a-b）.dist（）<
=r1）returnpi*r2*r2;
doublec1,c2;
c1=（r1*r1-r2*r2+（a-b）.dis（））/（a-b）.dist（）/r1/2.0;
c2=（r2*r2-r1*r1+（a-b）.dis（））/（a-b）.dist（）/r2/2.0;
doubles1,s2;
s1=acos（c1）;
s2=acos（c2）;
doubleans=0;
ans+=s1*r1*r1-r1*r1*sin（s1）*cos（s1）;
ans+=s2*r2*r2-r2*r2*sin（s2）*cos（s2）;
returnans;
多边形和圆相交的面积
structpoint{
doublex,y;
point（）{}
point（double_x,double_y）:
x（_x）,y（_y）{}
doublelen（）{returnsqrt（x*x+y*y）;
voidoutput（）{printf（"
%.15lf%.15lf\n"
x,y）;
}a,b,c,o;
constdoubleeps=1e-8;
constdoublePI=acos（-1.）;
doubler;
inlineintsign（doublex）{
if（x<
eps）return-1;
elsereturn（x>
eps）;
pointoperator*（double&
b）{
returnpoint（a*b.x,a*b.y）;
doubledot（constpoint&
returna.x*b.x+a.y*b.y;
doubledet（constpoint&
returna.x*b.y-a.y*b.x;
//===========用有向面积，划分成一个三角形和圆的面积的交
doublearea2（pointpa,pointpb）{
if（pa.len（）<
pb.len（））swap（pa,pb）;
if（pb.len（）<
eps）return0;
doublea,b,c,B,C,sinB,cosB,sinC,cosC,S,h,theta;
a=pb.len（）;
b=pa.len（）;
c=（pb-pa）.len（）;
//sinB=abs（det（pb,pb-pa））/a/c;
cosB=dot（pb,pb-pa）/a/c;
B=acos（cosB）;
//sinC=abs（det（pa,pb））/a/b;
cosC=dot（pa,pb）/a/b;
C=acos（cosC）;
//printf（"
area2（%.4lf,%.4lf,%.4lf）\n"
a,b,C/PI*180）;
if（a>
r）{
S=C/2*r*r;
h=a*b*sin（C）/c;
if（h<
r&
B<
PI/2）S-=（acos（h/r）*r*r-h*sqrt（r*r-h*h））;
elseif（b>
theta=PI-B-asin（sin（B）/r*a）;
S=.5*a*r*sin（theta）+（C-theta）/2*r*r;
elseS
升级会员 