AC » Archive » define型计算几何模板

define型计算几何模板

二月 23rd, 2011 › 图片, 心血来潮, 模板 › ben1222 ›

本想在寒假期间写点什么的…无奈域名出了点问题…
导致这一个半月里面网站一直处于无法访问的状态…

想起以前积累的计算几何模板…想差不多可以整理下发出来…
不过充斥着define…以及define的奇怪用法…可以说完全是C语言风格的代码…
总之这是一个基于数组和define而不是类和template的模板…会有人看吗…

如果每段代码都注解出一步步怎么得到的话…工程浩大…
结果心血来潮…给几乎每段代码画了个图片…
蓝色是已知,红色是所求…
画了几个才发现…有些还真难用图片表达啊…

那么就开始吧…

先是关于点,线,圆的形式的说明

说明:
point,p,vector,v,点,向量:
  array[2]  [0]:x,dx  [1]:y,dy
line,vl,直线:
  array[3]  [0]:a  [1]:b  [2]:c
  对应直线方程ax+by+c=0
circle,o,圆:
  array[3]  [0]:x  [1]:y  [2]:r

常量:圆周率和微小值

#define PI 3.14159265358979
#define EPS 1e-10

然后是向量运算

#define IM(v1,v2) (v1[0]*v2[0]+v1[1]*v2[1])
#define IM2(ps1,pe1,ps2,pe2) ((pe1[0]-ps1[0])*(pe2[0]-ps2[0])+(pe1[1]-ps1[1])*(pe2[1]-ps2[1]))
/* 外积大于0,则v2在v1左手方向 */
#define OM(v1,v2) (v1[0]*v2[1]-v1[1]*v2[0])
#define OM2(ps1,pe1,ps2,pe2) ((pe1[0]-ps1[0])*(pe2[1]-ps2[1])-(pe1[1]-ps1[1])*(pe2[0]-ps2[0]))
#define SU(v1,v2,v3) v3[0]=v2[0]-v1[0];v3[1]=v2[1]-v1[1];

距离:

/*两点间距离平方*/
#define DIS2(p1,p2,ret) {double __v[2];SU(p1,p2,__v) ret=IM(__v,__v);}
/*两点间距离*/
double dis(long*p1,long*p2){
	double v[2];
	SU(p1,p2,v);
	return sqrt(IM(v,v));
}

关于圆:
三点求圆心:


/* 三点求圆心 */
#define SQS(n,d) n=(d)*(d)
#define ARCC(x1,y1,x2,y2,x3,y3,x0,y0) {\
	double SQS(__y1,y1),SQS(__y2,y2),SQS(__y3,y3),SQS(__x1,x1),SQS(__x2,x2),SQS(__x3,x3); \
	x0=((y3-y1)*(__y2-__y1+__x2-__x1)+(y2-y1)*(__y1-__y3+__x1-__x3))*0.5/ \
		(double)((x2-x1)*(y3-y1)-(x3-x1)*(y2-y1)); \
	y0=((x3-x1)*(__x2-__x1+__y2-__y1)+(x2-x1)*(__x1-__x3+__y1-__y3))*0.5/ \
		(double)((y2-y1)*(x3-x1)-(y3-y1)*(x2-x1));}
#define ARCCP(p1,p2,p3,pc) ARCC(p1[0],p1[1],p2[0],p2[1],p3[0],p3[1],pc[0],pc[1])

两点和直径求圆心:


/*
给定圆上两点和圆的直径求2个候选的圆心位置,不判断两点间距离是否超过直径
使用:math.h
*/
#define PDCC(x1,y1,x2,y2,d,rx1,ry1,rx2,ry2) {\
	double SQS(_x,x2-x1),SQS(_y,y2-y1);\
	double _u=sqrt(d*d/(_x+_y)-1);\
	rx1=((x1+x2)+(y2-y1)*_u)/2;\
	ry1=((y1+y2)-(x2-x1)*_u)/2;\
	rx2=((x1+x2)-(y2-y1)*_u)/2;\
	ry2=((y1+y2)+(x2-x1)*_u)/2;}
#define PDCCP(p1,p2,d,r1,r2) PDCC(p1[0],p1[1],p2[0],p2[1],d,r1[0],r1[1],r2[0],r2[1])

圆外一点到圆的切点:


/*
求圆外一点到圆的2条切线的2个切点
po是圆心坐标,pp是圆外的一点,r是圆的半径,ret1和ret2存放2个切点坐标
不检测点是否在圆内
使用:向量运算,math.h
*/
#define TPC(p1,op,p2) (po[p1]+r*(r*d1[p1] op d1[p2]*sd)/s1)
void tanpoint(double*po,double*pp,double r,double*ret1,double*ret2){
	double d1[2],s1,sd;
	SU(po,pp,d1)
	s1=IM(d1,d1);
	sd=s1-r*r;
	if(sd<0)sd=-sd;
	sd=sqrt(sd);
	ret1[0]=TPC(0,+,1);
	ret1[1]=TPC(1,-,0);
	ret2[0]=TPC(0,-,1);
	ret2[1]=TPC(1,+,0);
}

圆相交面积:


/*
圆相交面积
使用:向量运算,距离,常量,math.h
*/
double cosine(double b1, double b2, double b3){
	double b = b2 * b2 + b3 * b3 - b1 * b1;
	return b / b2 / b3 / 2;
}
double circle_ins(double*o1,double*o2){
	double dc=dis(o1,o2),r2=(o1[2]>o2[2]?o1[2]:o2[2]),r1=(o1[2]<o2[2]?o1[2]:o2[2]);
	if(r2-r1>=dc)return PI*r1*r1;
	else if(dc>=r2+r1)return 0;
	else{
		double a1,a2,ang1,ang2,ret=0,s1=PI*r1*r1,s2=PI*r2*r2;
		a1=cosine(r2,r1,dc);
		a2=cosine(r1,r2,dc);
		ang1=acos(a1);
		ang2=acos(a2);
		if(a1>0)ret=(r1*r1*ang1-r1*r1*sin(ang1*2)/2);
		else ret=(r1*r1*ang1+r1*r1*sin((PI-ang1)*2)/2);
		if(a2>0)ret+=(r2*r2*ang2-r2*r2*sin(ang2*2)/2);
		else ret+=(r2*r2*ang2+r2*r2*sin((PI-ang2)*2)/2);
		return ret;
	}
}

关于直线:
直线解析式:


/*
给定两点,求其所在直线的解析式
ret[0]=a,ret[1]=b,ret[2]=c
ax+by+c=0
使用:向量运算
*/
#define PT_LINE(p1,p2,ret) {ret[1]=p2[0]-p1[0];ret[0]=p1[1]-p2[1];ret[2]=-IM(p1,ret);}

点关于直线的对称点:


/*
点关于直线的对称点
vl:直线解析式,pt:点,ret:返回对称点
使用:向量运算
*/
#define PL_SYM(vl,pt,ret) {\
	double __d=IM(vl,vl); \
	ret[0]=((vl[1]*vl[1]-vl[0]*vl[0])*pt[0]-2*vl[0]*(vl[1]*pt[1]+vl[2]))/__d; \
	ret[1]=((vl[0]*vl[0]-vl[1]*vl[1])*pt[1]-2*vl[1]*(vl[0]*pt[0]+vl[2]))/__d;}

直线交点:


/*
判断两条直线[ps1,pe1],[ps2,pe2]的交点情况:
	0:平行 1:相交 2:重叠
如果有交点,返回值在ret中
使用:向量运算
*/
long lineins(long*ps1,long*pe1,long*ps2,long*pe2,double*ret){
	long v1[2],v2[2];
	v1[0]=pe1[0]-ps1[0];v1[1]=pe1[1]-ps1[1];
	v2[0]=pe2[0]-ps2[0];v2[1]=pe2[1]-ps2[1];
	long mul=OM(v1,v2);
	if(mul==0){
		v2[0]=pe2[0]-ps1[0];
		v2[1]=pe2[1]-ps1[1];
		return (OM(v1,v2)==0)?2:0;
	}else{
		double t=((ps2[0]-ps1[0])*v2[1]+(ps1[1]-ps2[1])*v2[0])/(double)mul;
		ret[0]=ps1[0]+v1[0]*t;
		ret[1]=ps1[1]+v1[1]*t;
		return 1;
	}
}

点到直线的距离:


/*
点pt到直线[pl1,pl2]的距离平方
使用:向量运算
*/
double dis2pl(double*pt,double*pl1,double*pl2){
	double v1[2],v2[2];
	SU(pt,pl1,v1)SU(pl1,pl2,v2)
	double om=OM(v1,v2);
	return om*om/IM(v2,v2);
}

直线和圆的交点:


/*
直线[pl1,pl2]与圆[po,r]的交点,返回值为交点个数
ax+by+c=0;
(x-ox)^2+(y-oy)^2=r^2
x=( b(b*ox-a*oy)-a*c土b*sqrt(r^2-dis2pl()^2) )/( a^2+b^2 )
y=( a(a*oy-b*ox)-b*c干a*sqrt(r^2-dis2pl()^2) )/( a^2+b^2 )
使用:向量运算,常量,直线解析式,math.h
*/
long lcins(double*po,double r,double*pl1,double*pl2,double*ret1,double*ret2){
	double vl[3],tmp1,tmp2;
	PT_LINE(pl1,pl2,vl)
	tmp1=IM(vl,po)+vl[2];
	double ab=IM(vl,vl),mdis2=ab*r*r-tmp1*tmp1;
	if(mdis2<=-EPS)return 0;
	else{
		tmp1=(vl[1]*OM(po,vl)-vl[0]*vl[2]);
		tmp2=(vl[0]*OM(vl,po)-vl[1]*vl[2]);
		if(mdis2<EPS){
			ret1[0]=tmp1/ab;ret1[1]=tmp2/ab;
			return 1;
		}else{
			double sq=sqrt(mdis2);
			ret1[0]=(tmp1+vl[1]*sq)/ab;ret1[1]=(tmp2-vl[0]*sq)/ab;
			ret2[0]=(tmp1-vl[1]*sq)/ab;ret2[1]=(tmp2+vl[0]*sq)/ab;
			return 2;
		}
	}
}

直线和椭圆的交点:


/*
直线[pl1,pl2]与椭圆的交点,返回值为交点个数
ax+by+c=0
(x/d)^2+(y/e)^2=1
x=( -c*ad*d土be*d*sqrt(ad^2+be^2-c^2) )/( ad^2+be^2 )
y=( -c*be*e干ad*e*sqrt(ad^2+be^2-c^2) )/( ad^2+be^2 )
使用:向量运算,常量,直线解析式,math.h
*/
long leins(double d,double e,double*pl1,double*pl2,double*ret1,double*ret2){
	double vl[3];
	PT_LINE(pl1,pl2,vl)
	double ab=vl[0]*vl[0]*d*d+vl[1]*vl[1]*e*e,delta=ab-vl[2]*vl[2];
	if(delta<=-EPS)return 0;
	else if(delta<EPS){
		ret1[0]=-vl[2]*vl[0]*d*d/ab;
		ret1[1]=-vl[2]*vl[1]*e*e/ab;
		return 1;
	}else{
		double sq=sqrt(delta);
		ret1[0]=(-vl[2]*vl[0]*d*d+vl[1]*e*d*sq)/ab;
		ret1[1]=(-vl[2]*vl[1]*e*e-vl[0]*e*d*sq)/ab;
		ret2[0]=(-vl[2]*vl[0]*d*d-vl[1]*e*d*sq)/ab;
		ret2[1]=(-vl[2]*vl[1]*e*e+vl[0]*e*d*sq)/ab;
		return 2;
	}
}

关于线段:

/*
整数点间连一条线,问穿过多少个单位正方形:
	=a+b-gcd(a,b)
扩展到三维:
	=a+b+c-gcd(a,b)-gcd(a,c)-gcd(b,c)+gcd(a,b,c)
*/

判断线段是否相交:


/*
判断两线段是否相交,跨立实验中用>则不包括端点,用>=则包括端点
使用:向量运算
*/
#define MAX(a,b) ((a)>(b)?(a):(b))
#define MIN(a,b) ((a)<(b)?(a):(b))
#define MN(a,b,t) (MAX(ps##a[t],pe##a[t])>=MIN(ps##b[t],pe##b[t]))
long segins(long*ps1,long*pe1,long*ps2,long*pe2){
	return MN(1,2,0)&&MN(2,1,0)&&MN(1,2,1)&&MN(2,1,1) && //快速排斥实验
		(OM2(ps1,ps2,ps1,pe1)*OM2(ps1,pe1,ps1,pe2)>=0) &&
		(OM2(ps2,ps1,ps2,pe2)*OM2(ps2,pe2,ps2,pe1)>=0); //跨立实验
}

点到线段的垂足:


/*
点pt到线段[ps,pe]所在直线的垂足,返回true表示垂足在线段上,垂足坐标存在ret中
使用:向量运算
*/
long plfoot(double*pt,double*ps,double*pe,double*ret){
	double vt[2],t,ret[2];
	SU(ps,pe,vt)
	t=(-(ps[0]-pt[0])*(pe[0]-ps[0])-(ps[1]-pt[1])*(pe[1]-ps[1]))/IM(vt,vt);
	ret[0]=ps[0]+t*(pe[0]-ps[0]);
	ret[1]=ps[1]+t*(pe[1]-ps[1]);
	return !((ret[0]<ps[0]&&ret[0]<pe[0])||(ret[0]>ps[0]&&ret[0]>pe[0])||
		(ret[1]<ps[1]&&ret[1]<pe[1])||(ret[1]>ps[1]&&ret[1]>pe[1]));
}

点到线段的距离:


/*
点pt到线段[pl1,pl2]的距离的平方
使用:向量运算,常量,距离,点到直线距离
*/
double dis2ps(double*pt,double*pl1,double*pl2){
	double ret;
	if(IM2(pl1,pt,pl1,pl2)<-EPS)DIS2(pt,pl1,ret)
	else if(IM2(pl2,pt,pl2,pl1)<-EPS)DIS2(pt,pl2,ret)
	else ret=dis2pl(pt,pl1,pl2);
	return ret;
}

角平分线上一点:


/*
求2个线段[pc,pl]和[pc,pr]组成的夹角的角平分线上的一个点
使用:向量运算,math.h
*/
void hangle(double*pl,double*pc,double*pr,double*ret){
	double vl[2],vr[2];
	SU(pc,pl,vl)
	SU(pc,pr,vr)
	double dl=sqrt(IM(vl,vl)),dr=sqrt(IM(vr,vr));
	vl[0]/=dl;vl[1]/=dl;
	vr[0]/=dr;vr[1]/=dr;
	ret[0]=(vl[0]+vr[0])/2.0+pc[0];
	ret[1]=(vl[1]+vr[1])/2.0+pc[1];
}

关于三角形:
三角形面积:


/*
三角形面积之海伦公式:
p=(a+b+c)/2
s=sqrt((p-a)*(p-b)*(p-c)*p)
*/

三角形内切圆半径:


/*
三角形内切圆半径:
直角三角形:r=(a+b-c)/2,其中a,b是直角边长,c是斜边长
一般三角形:r=2S/(a+b+c),其中S是三角形面积
*/

三角形内周长一定时的最大面积:


/*
三角形内周长一定的最大面积
triple=x+y+z;
half=triple*0.5;
area=sqrt(half*(half-x)*(half-y)*(half-z));
R=area*2.0/triple;
if(x+y+z<=role)
	max=area;
else if(2.0*PI*R>=role)
	max=role*role/(4.0*PI);
else
{
	r=(triple-role)/(triple/R-2.0*PI);
	max=area+PI*r*r-(r*r*area/(R*R));
}

*/

三角形重心:

/* 三角形重心 */
#define CGT(p1,p2,p3,pc) {pc[0]=(p1[0]+p2[0]+p3[0])/3;pc[1]=(p1[1]+p2[1]+p3[1])/3;}

判断点是否在三角形内:


/*
外积判断点是否在三角形内
在三角形内的返回2,边上的1,外面的0
使用:向量运算
*/
long pint(long*p1,long*p2,long*p3,long*pm){
	long v1[2],v2[2],t,l;
	SU(p1,p2,v1)SU(p1,pm,v2)
	if((t=OM(v1,v2))==0)return 1;
	t=t>0?1:-1;
	SU(p2,p3,v1)SU(p2,pm,v2)
	l=t*OM(v1,v2);
	if(l==0)return 1;else if(l<0)return 0;
	SU(p3,p1,v1)SU(p3,pm,v2)
	l=t*OM(v1,v2);
	if(l==0)return 1;else if(l<0)return 0;
	return 2;
}

三角形内一点到3边的距离的平方和最小:(这个图片我实在想不出怎么画…)

/*
三角形内一点到3边的距离的平方和最小
void (double pts[3][2])
	double x,y,vl[3][3],d2[3];
	long i;
	for(i=0;i<3;i++){
		PT_LINE(pts[i],pts[(i+1)%3],vl[i]);
		d2[i]=IM(vl[i],vl[i]);
	}
	double X=0,Y=0,C=0,U=0,B=0;
	for(i=0;i<3;i++){
		X+=vl[i][0]*vl[i][1]/d2[i];
		Y+=vl[i][1]*vl[i][1]/d2[i];
		C+=vl[i][1]*vl[i][2]/d2[i];
	}
	for(i=0;i<3;i++){
		double t=vl[i][0]*Y-vl[i][1]*X;
		U+=t*(vl[i][2]*Y-vl[i][1]*C)/d2[i];
		B+=t*t/d2[i];
	}
	x=-U/B;
	y=-(X*x+C)/Y;
*/

判断矩形能否斜放在另一个矩形内:


/*
判断矩形能否斜放在另一个矩形内
minh=(2xya+d*(x*x-y*y))/(x*x+y*y)  (x>y)
d=sqrt(x*x+y*y-a*a)
return minh<=b

不判断面积大小和快速排斥
使用:math.h
*/
#define I __int64
long rcinrc(long w1,long h1,long w2,long h2){
	if((w1<h1?w1:h1)<(w2<h2?w2:h2))return 0;
	I sq=(I)w2*(I)w2+(I)h2*(I)h2,l=2*(I)w2*(I)h2*(I)w1,a2=(I)w1*(I)w1;
	if(w2==h2||sq==a2)return l<=sq*h1;
	else{
		if(w2<h2){long t=w2;w2=h2;h2=t;}
		if(sq<a2)return 0;
		I il=((I)w2*(I)w2-(I)h2*(I)h2);
		double dl=l+sqrt((double)(sq-a2))*il,dr=sq*h1;
		return dl<=dr;
	}
}

三角形和圆的公共部分面积:


/*
点pt在线段[pl1,pl2]的位置
IM2(pl1,pt,pl1,pl2)<-EPS  点在线段外,pl1一侧
IM2(pl2,pt,pl2,pl1)<-EPS  点在线段外,pl2一侧
PTINSEG:点在线段所在直线的前提上,判断点是否在线段内,线段加长一段微小量
PTINSEGS:同上,线段缩小一段微小量
使用:向量运算,常量
*/
#define PTINSEG(pt,pl1,pl2) ((IM2(pl1,pt,pl1,pl2)>=-EPS)&&(IM2(pl2,pt,pl2,pl1)>=-EPS))
#define PTINSEGS(pt,pl1,pl2) ((IM2(pl1,pt,pl1,pl2)>EPS)&&(IM2(pl2,pt,pl2,pl1)>EPS))

/*
2点和圆心组成的三角形与这个圆的公共部分面积
r是圆半径,原点是圆心
ps和pe是三角形另外两个点,逆时针为正
使用:向量运算,常量,直线和圆的交点,math.h
*/
#define VRAD(v,ret) {if(v[0]==0){if(v[1]>0)ret=PI/2;else ret=PI*3/2;}\
	else{double _t=atan(v[1]/v[0]);if(v[0]<0)_t=PI+_t;ret=(_t<0?_t+PI*2:_t);}}
#define MARC(rad) {if(rad>PI)rad-=PI*2;if(rad<-PI)rad+=PI*2;}
double area_c2p(double r,double*ps,double*pe){
	double d1=IM(ps,ps),d2=IM(pe,pe),r2=r*r;
	if(d1<=r2&&d2<=r2)return OM(ps,pe)/2;
	double rad1,rad2,po[2],ins1[2],ins2[2],drad;
	VRAD(ps,rad1)VRAD(pe,rad2)
	drad=rad2-rad1;MARC(drad)
	po[0]=po[1]=0;
	long nins=lcins(po,r,ps,pe,ins1,ins2);
	if(nins<=1)return drad/2*r*r;
	if((IM2(ps,ins1,ps,pe)<-EPS&&IM2(ps,ins2,ps,pe)<-EPS)||\
		(IM2(pe,ins1,pe,ps)<-EPS&&IM2(pe,ins2,pe,ps)<-EPS))return drad/2*r*r;
	else if(PTINSEG(ins1,ps,pe)&&PTINSEG(ins2,ps,pe)){
		double radi1,radi2,drad2,ret;
		if(PTINSEGS(ins1,ps,ins2)){VRAD(ins1,radi1)VRAD(ins2,radi2)ret=OM(ins1,ins2)/2;}
		else{VRAD(ins2,radi1)VRAD(ins1,radi2)ret=OM(ins2,ins1)/2;}
		drad=radi1-rad1;MARC(drad)
		drad2=rad2-radi2;MARC(drad2)
		return ret+(drad+drad2)/2*r*r;
	}else if(PTINSEG(ps,ins1,ins2)&&PTINSEG(pe,ins1,ins2))return OM(ps,pe)/2;
	else{
		double*pins,*pout,radi;
		if(IM2(ps,ins1,ps,pe)<-EPS||IM2(pe,ins1,pe,ps)<-EPS){pins=ins2;pout=ins1;}
		else{pins=ins1;pout=ins2;}
		VRAD(pins,radi)
		if(IM2(pe,pout,pe,ps)<-EPS){
			drad=radi-rad1;MARC(drad)
			return OM(pins,pe)/2+drad/2*r*r;
		}else{
			drad=rad2-radi;MARC(drad)
			return OM(ps,pins)/2+drad/2*r*r;
		}
	}
}

凸包:


/*排序用*/
#define POINT struct TPoint
POINT{long xy[2];}pstart;
long operator<(POINT a,POINT b){
	long t=OM2(pstart.xy,a.xy,pstart.xy,b.xy);
	if(t==0){
		long d1,d2;
		DIS2(pstart.xy,a.xy,d1)
		DIS2(pstart.xy,b.xy,d2)
		return d1<d2;
	}else return t>0;
}
/*
O(nlogn)的凸包,返回值是凸多边形的顶点数,存在pts中,逆时针方向
使用:向量运算,距离,std::sort(algorithm)
*/
long convex(POINT*pts,long num,long*stack){
	long i,minp,minx=0x7fffffff,miny;
	for(i=0;i<num;i++){
		if(pts[i].xy[0]<minx){minx=pts[i].xy[0];miny=pts[i].xy[1];minp=i;}
		else if(pts[i].xy[0]==minx&&pts[i].xy[1]<miny)miny=pts[minp=i].xy[1];
	}
	if(minx==0x7fffffff)return 0;
	POINT _t=pts[minp];
	pts[minp]=pts[0];
	pts[0]=_t;
	pstart=pts[0];
	sort(pts+1,pts+num);
	long top=3;
	for(i=0;i<3;i++)stack[i+1]=i;
	for(i=3;i<num;i++){
		while(OM2(pts[stack[top-1]].xy,pts[stack[top]].xy,pts[stack[top-1]].xy,pts[i].xy)<0)top--;
		stack[++top]=i;
	}
	for(i=1;i<=top;i++)pts[i-1]=pts[stack[i]];
	return top;
}