计算几何学函数

最新推荐文章于 2023-03-11 09:49:11 发布

winepot

最新推荐文章于 2023-03-11 09:49:11 发布

阅读量626

点赞数 1

分类专栏：技术文章文章标签： distance float struct ini c

本文链接：https://blog.csdn.net/winepot/article/details/2606390

版权

技术文章专栏收录该内容

9 篇文章 0 订阅

订阅专栏

/*********************************************/ * * * 计算几何学库函数 * * * * * /*********************************************/ #include <stdlib.h> #define infinity 1e20 #define EP 1e-10 struct TPoint{ float x,y; }; struct TLineSeg{ TPoint a,b; }; //求平面上两点之间的距离 float distance(TPoint p1,TPoint p2) { return(sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y))); } /******************************************************** 返回(P1-P0)*(P2-P0)的叉积。若结果为正，则<P0,P1>在<P0,P2>的顺时针方向；若为0则<P0,P1><P0,P2>共线；若为负则<P0,P1>在<P0,P2>的在逆时针方向; 可以根据这个函数确定两条线段在交点处的转向, 比如确定p0p1和p1p2在p1处是左转还是右转，只要求 (p2-p0)*(p1-p0)，若<0则左转，>0则右转，=0则共线 *********************************************************/ float multiply(TPoint p1,TPoint p2,TPoint p0) { return((p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y)); } //确定两条线段是否相交 int intersect(TLineSeg u,TLineSeg v) { return( (max(u.a.x,u.b.x)>=min(v.a.x,v.b.x))&& (max(v.a.x,v.b.x)>=min(u.a.x,u.b.x))&& (max(u.a.y,u.b.y)>=min(v.a.y,v.b.y))&& (max(v.a.y,v.b.y)>=min(u.a.y,u.b.y))&& (multiply(v.a,u.b,u.a)*multiply(u.b,v.b,u.a)>=0)&& (multiply(u.a,v.b,v.a)*multiply(v.b,u.b,v.a)>=0)); } //判断点p是否在线段l上 int online(TLineSeg l,TPoint p) { return( (multiply(l.b,p,l.a)==0)&&( ((p.x-l.a.x)*(p.x-l.b.x)<0 )||( (p.y-l.a.y)*(p.y-l.b.y)<0 )) ); } //判断两个点是否相等 int Euqal_Point(TPoint p1,TPoint p2) { return((abs(p1.x-p2.x)<EP)&&(abs(p1.y-p2.y)<EP)); } //一种线段相交判断函数，当且仅当u,v相交并且交点不是u,v的端点时函数为true; int intersect_A(TLineSeg u,TLineSeg v) { return((intersect(u,v))&& (!Euqal_Point(u.a,v.a))&& (!Euqal_Point(u.a,v.b))&& (!Euqal_Point(u.b,v.a))&& (!Euqal_Point(u.b,v.b))); } /*=============================================== 判断点q是否在多边形Polygon内，其中多边形是任意的凸或凹多边形， Polygon中存放多边形的逆时针顶点序列 ================================================*/ int InsidePolygon(int vcount,TPoint Polygon[],TPoint q) { int c=0,i,n; TLineSeg l1,l2; l1.a=q; l1.b=q; l1.b.x=infinity; n=vcount; for (i=0;i<vcount;i++) { l2.a=Polygon[i]; l2.b=Polygon[(i+1)%n]; if ( (intersect_A(l1,l2))|| ( (online(l1,Polygon[(i+1)%n]))&& ( (!online(l1,Polygon[(i+2)%n]))&& (multiply(Polygon[i],Polygon[(i+1)%n],l1.a)*multiply(Polygon[(i+1)%n],Polygon[(i+2)%n],l1.a)>0) || (online(l1,Polygon[(i+2)%n]))&& (multiply(Polygon[i],Polygon[(i+2)%n],l1.a)*multiply(Polygon[(i+2)%n],Polygon[(i+3)%n],l1.a)>0) ) ) ) c++; } return(c%2!=0); } /********************************************/ * * * 计算多边形的面积 * * * * 要求按照逆时针方向输入多边形顶点 * * 可以是凸多边形或凹多边形 * * * /********************************************/ float area_of_polygon(int vcount,float x[],float y[]) { int i; float s; if (vcount<3) return 0; s=y[0]*(x[vcount-1]-x[1]); for (i=1;i<vcount;i++) s+=y[i]*(x[(i-1)]-x[(i+1)%vcount]); return s/2; } /*==================================================================================== 寻找凸包 graham 扫描法 PointSet为输入的点集； ch为输出的凸包上的点集，按照逆时针方向排列; n为PointSet中的点的数目 len为输出的凸包上的点的个数；设下一个扫描的点PointSet[i]=P2,当前栈顶的两个点ch[top]=P1,ch[top-1]=P0, 如果P1P2相对于P0P1在点P1向左旋转(共线也不行)，则P0,P1一定是凸包上的点；否则P1一定不是凸包上的点，应该将其出栈。比如下面左图中点1就一定是凸包上的点，右图中点1就一定不是凸包上的点，因为可以连接点0,2构成凸包的边 2 | | _____2 1 1 / / ____0/ ____0/ ====================================================================================*/ void Graham_scan(TPoint PointSet[],TPoint ch[],int n,int &len) { int i,j,k=0,top=2; TPoint tmp; //选取PointSet中y坐标最小的点PointSet[k]，如果这样的点右多个，则取最左边的一个 for(i=1;i<n;i++) if ((PointSet[i].y<PointSet[k].y)||((PointSet[i].y==PointSet[k].y)&&(PointSet[i].x<PointSet[k].x))) k=i; tmp=PointSet[0]; PointSet[0]=PointSet[k]; PointSet[k]=tmp; //现在PointSet中y坐标最小的点在PointSet[0] for (i=1;i<n-1;i++) //对顶点按照相对PointSet[0]的极角从小到大进行排序，极角相同的按照距离PointSet[0]从近到远进行排序 { k=i; for (j=i+1;j<n;j++) if ( (multiply(PointSet[j],PointSet[k],PointSet[0])>0) || ( (multiply(PointSet[j],PointSet[k],PointSet[0])==0) &&(distance(PointSet[0],PointSet[j])<distance(PointSet[0],PointSet[k])) ) ) k=j; tmp=PointSet[i]; PointSet[i]=PointSet[k]; PointSet[k]=tmp; } ch[0]=PointSet[0]; ch[1]=PointSet[1]; ch[2]=PointSet[2]; for (i=3;i<n;i++) { while (multiply(PointSet[i],ch[top],ch[top-1])>=0) top--; ch[++top]=PointSet[i]; } len=top+1; }