交大计算几何模板.docx
- 文档编号:8271687
- 上传时间:2023-01-30
- 格式:DOCX
- 页数:51
- 大小:37.11KB
交大计算几何模板.docx
《交大计算几何模板.docx》由会员分享,可在线阅读,更多相关《交大计算几何模板.docx(51页珍藏版)》请在冰豆网上搜索。
交大计算几何模板
计算几何_main
#include
#include
#include
#include
#include
#include
#include
usingnamespacestd;
#definerep(i,n)for(i=0;i<(n);i++)
#defineforu(i,a,b)for(i=(a);i<=(b);i++)
#defineford(i,a,b)for(i=(a);i>=(b);i--)
doubleeps=1e-8;
structline{doublea,b,c;};
intcmp(doublex){
if(x>eps)return1;
if(x<-eps)return-1;
return0;
}
classpoint{
public:
doublex,y;
point(){}
point(doublex,doubley):
x(x),y(y){}
voidinput(){scanf("%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;returna;}
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);
returnh;
}
//================点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&p0,constpoint&p1,constpoint&p2,point&cp)
{
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){
pointh;
h.y=-(a.c*b.a-b.c*a.a)/(a.b*b.a-b.b*a.a);//=====makesureaandbaren'tparallel
if(abs(a.a) elseh.x=(-a.c-h.y*a.b)/a.a; returnh; } //==========线段平移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<0&&d2<=0&&d1>0)wn--; } returnwn! =0; } //========================================判断线段是否有在多边形内部========= boolcompareab(constpoint&a,constpoint&b){ if(a.x 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); a[n]=a[0]; 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; } } } closed++;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; } returnfalse; } //多边形的重心 voidPolygonCentroids(pointp[],intn,point&cp){ //if面积为0需要特判 doublesum=0,s=0;cp.x=0;cp.y=0; for(inti=1;i 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 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&cp){ 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&cp){ 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: (242,89),(212,185),(71,128),OUTPUT: (208.8229,171.0697) voidOrthocenter(pointA,pointB,pointC,point&cp){ 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 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&cp1,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 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>eps)return2;elsereturn1; } //====两园面积交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 } pointoperator*(double&a,constpoint&b){ returnpoint(a*b.x,a*b.y); } doubledot(constpoint&a,constpoint&b){ returna.x*b.x+a.y*b.y; } doubledet(constpoint&a,constpoint&b){ returna.x*b.y-a.y*b.x; } //===========用有向面积,划分成一个三角形和圆的面积的交 doublearea2(pointpa,pointpb){ if(pa.len() if(pb.len() 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 } elseif(b>r){ theta=PI-B-asin(sin(B)/r*a); S=.5*a*r*sin(theta)+(C-theta)/2*r*r; } elseS
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 交大 计算 几何 模板