计算几何目录 ㈠ 点的基本运算 1. 平面上两点之间距离 1 2. 判断两点是否重合 1 3. 矢量叉乘 1 4. 矢量点乘 2 5. 判断点是否在线段上 2 6. 求一点饶某点旋转后的坐标 2 7. 求矢量夹角 2 ㈡ 线段及直线的基本运算 1. 点与线段的关系 3 2. 求点到线段所在直线垂线的垂足 4 3. 点到线段的最近点 4 4. 点到线段所在直线的距离 4 5. 点到折线集的最近距离 4 6. 判断圆是否在多边形内 5 7. 求矢量夹角余弦 5 8. 求线段之间的夹角 5 9. 判断线段是否相交 6 10.判断线段是否相交但不交在端点处 6 11.求线段所在直线的方程 6 12.求直线的斜率 7 13.求直线的倾斜角 7 14.求点关于某直线的对称点 7 15.判断两条直线是否相交及求直线交点 7 16.判断线段是否相交，如果相交返回交点 7 ㈢ 多边形常用算法模块 1. 判断多边形是否简单多边形 8 2. 检查多边形顶点的凸凹性 9 3. 判断多边形是否凸多边形 9 4. 求多边形面积 9 5. 判断多边形顶点的排列方向，方法一 10 6. 判断多边形顶点的排列方向，方法二 10 7. 射线法判断点是否在多边形内 10 8. 判断点是否在凸多边形内 11 9. 寻找点集的graham算法 12 10.寻找点集凸包的卷包裹法 13 11.判断线段是否在多边形内 14 12.求简单多边形的重心 15 13.求凸多边形的重心 17 14.求肯定在给定多边形内的一个点 17 15.求从多边形外一点出发到该多边形的切线 18 16.判断多边形的核是否存在 19 ㈣ 圆的基本运算 1 .点是否在圆内 20 2 .求不共线的三点所确定的圆 21 ㈤ 矩形的基本运算 1.已知矩形三点坐标，求第4点坐标 22 ㈥ 常用算法的描述 22 ㈦ 补充 1．两圆关系： 24 2．判断圆是否在矩形内： 24 3．点到平面的距离： 25 4．点是否在直线同侧： 25 5．镜面反射线： 25 6．矩形包含： 26 7．两圆交点： 27 8．两圆公共面积： 28 9. 圆和直线关系： 29 10. 内切圆： 30 11. 求切点： 31 12. 线段的左右旋： 31 13．公式： 32 *//* 需要包含的头文件 */ #include 
     
       /* 常用的常量定义 */ const double INF  = 1E200    const double EP  = 1E-10 const int  MAXV = 300 const double PI  = 3.14159265 /* 基本几何结构 */ struct POINT {  double x;  double y;  POINT(double a=0, double b=0) { x=a; y=b;} //constructor }; struct LINESEG {  POINT s;  POINT e;  LINESEG(POINT a, POINT b) { s=a; e=b;}  LINESEG() { } }; struct LINE           // 直线的解析方程 a*x+b*y+c=0  为统一表示，约定 a >= 0 {    double a;    double b;    double c;    LINE(double d1=1, double d2=-1, double d3=0) {a=d1; b=d2; c=d3;} }; /********************** *                    *  *   点的基本运算     *  *                    *  **********************/ double dist(POINT p1,POINT p2)                // 返回两点之间欧氏距离 {  return( sqrt( (p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y) ) ); } bool equal_point(POINT p1,POINT p2)           // 判断两个点是否重合  {  return ( (abs(p1.x-p2.x)
      
       
        0：ep在矢量opsp的逆时针方向； r=0：opspep三点共线； r<0：ep在矢量opsp的顺时针方向 *******************************************************************************/ double multiply(POINT sp,POINT ep,POINT op) {  return((sp.x-op.x)*(ep.y-op.y)-(ep.x-op.x)*(sp.y-op.y)); } /* r=dotmultiply(p1,p2,op),得到矢量(p1-op)和(p2-op)的点积，如果两个矢量都非零矢量 r<0：两矢量夹角为钝角；r=0：两矢量夹角为直角；r>0：两矢量夹角为锐角 *******************************************************************************/ double dotmultiply(POINT p1,POINT p2,POINT p0) {  return ((p1.x-p0.x)*(p2.x-p0.x)+(p1.y-p0.y)*(p2.y-p0.y)); } /****************************************************************************** 判断点p是否在线段l上条件：(p在线段l所在的直线上) && (点p在以线段l为对角线的矩形内)*******************************************************************************/ bool online(LINESEG l,POINT p) {  return( (multiply(l.e,p,l.s)==0) &&( ( (p.x-l.s.x)*(p.x-l.e.x)<=0 )&&( (p.y-l.s.y)*(p.y-l.e.y)<=0 ) ) ); } // 返回点p以点o为圆心逆时针旋转alpha(单位：弧度)后所在的位置 POINT rotate(POINT o,double alpha,POINT p) {  POINT tp;  p.x-=o.x;  p.y-=o.y;  tp.x=p.x*cos(alpha)-p.y*sin(alpha)+o.x;  tp.y=p.y*cos(alpha)+p.x*sin(alpha)+o.y;  return tp; } /* 返回顶角在o点，起始边为os，终止边为oe的夹角(单位：弧度)  角度小于pi，返回正值  角度大于pi，返回负值  可以用于求线段之间的夹角 原理： r = dotmultiply(s,e,o) / (dist(o,s)*dist(o,e)) r'= multiply(s,e,o) r >= 1 angle = 0; r <= -1 angle = -PI -1
        
         <1 && r'>0 angle = arccos(r) -1
         
          <1 && r'<=0 angle = -arccos(r)*/ double angle(POINT o,POINT s,POINT e) {  double cosfi,fi,norm;  double dsx = s.x - o.x;  double dsy = s.y - o.y;  double dex = e.x - o.x;  double dey = e.y - o.y;  cosfi=dsx*dex+dsy*dey;  norm=(dsx*dsx+dsy*dsy)*(dex*dex+dey*dey);  cosfi /= sqrt( norm );  if (cosfi >=  1.0 ) return 0;  if (cosfi <= -1.0 ) return -3.1415926;  fi=acos(cosfi);  if (dsx*dey-dsy*dex>0) return fi;      // 说明矢量os 在矢量 oe的顺时针方向  return -fi; }   /*****************************\   *                             *   *      线段及直线的基本运算   *   *                             *   \*****************************/ /* 判断点与线段的关系,用途很广泛 本函数是根据下面的公式写的，P是点C到线段AB所在直线的垂足                 AC dot AB         r =     ---------                  ||AB||^2              (Cx-Ax)(Bx-Ax) + (Cy-Ay)(By-Ay)           = -------------------------------                           L^2     r has the following meaning:         r=0      P = A         r=1      P = B         r<0   P is on the backward extension of AB   r>1      P is on the forward extension of AB         0
          
           <1 P is interior to AB */ double relation(POINT p,LINESEG l) { LINESEG tl; tl.s=l.s; tl.e=p; return dotmultiply(tl.e,l.e,l.s)/(dist(l.s,l.e)*dist(l.s,l.e)); } // 求点C到线段AB所在直线的垂足 P POINT perpendicular(POINT p,LINESEG l) { double r=relation(p,l); POINT tp; tp.x=l.s.x+r*(l.e.x-l.s.x); tp.y=l.s.y+r*(l.e.y-l.s.y); return tp; } /* 求点p到线段l的最短距离,并返回线段上距该点最近的点np 注意：np是线段l上到点p最近的点，不一定是垂足 */ double ptolinesegdist(POINT p,LINESEG l,POINT &np) { double r=relation(p,l); if(r<0) { np=l.s; return dist(p,l.s); } if(r>1) { np=l.e; return dist(p,l.e); } np=perpendicular(p,l); return dist(p,np); } // 求点p到线段l所在直线的距离,请注意本函数与上个函数的区别 double ptoldist(POINT p,LINESEG l) { return abs(multiply(p,l.e,l.s))/dist(l.s,l.e); } /* 计算点到折线集的最近距离,并返回最近点. 注意：调用的是ptolineseg()函数 */ double ptopointset(int vcount,POINT pointset[],POINT p,POINT &q) { int i; double cd=double(INF),td; LINESEG l; POINT tq,cq; for(i=0;i
           
            = 0。//判断Q1Q2跨立P1P2的依据是：( Q1 - P1 ) × ( P2 - P1 ) * ( P2 - P1 ) × ( Q2 - P1 ) >= 0。bool intersect(LINESEG u,LINESEG v) { return( (max(u.s.x,u.e.x)>=min(v.s.x,v.e.x))&& //排斥实验 (max(v.s.x,v.e.x)>=min(u.s.x,u.e.x))&& (max(u.s.y,u.e.y)>=min(v.s.y,v.e.y))&& (max(v.s.y,v.e.y)>=min(u.s.y,u.e.y))&& (multiply(v.s,u.e,u.s)*multiply(u.e,v.e,u.s)>=0)&& //跨立实验 (multiply(u.s,v.e,v.s)*multiply(v.e,u.e,v.s)>=0)); } // (线段u和v相交)&&(交点不是双方的端点) 时返回true bool intersect_A(LINESEG u,LINESEG v) { return ((intersect(u,v))&& (!online(u,v.s))&& (!online(u,v.e))&& (!online(v,u.e))&& (!online(v,u.s))); } // 线段v所在直线与线段u相交时返回true；方法：判断线段u是否跨立线段v bool intersect_l(LINESEG u,LINESEG v) { return multiply(u.s,v.e,v.s)*multiply(v.e,u.e,v.s)>=0; } // 根据已知两点坐标，求过这两点的直线解析方程： a*x+b*y+c = 0 (a >= 0) LINE makeline(POINT p1,POINT p2) { LINE tl; int sign = 1; tl.a=p2.y-p1.y; if(tl.a<0) { sign = -1; tl.a=sign*tl.a; } tl.b=sign*(p1.x-p2.x); tl.c=sign*(p1.y*p2.x-p1.x*p2.y); return tl; } // 根据直线解析方程返回直线的斜率k,水平线返回 0,竖直线返回 1e200 double slope(LINE l) { if(abs(l.a) < 1e-20) return 0; if(abs(l.b) < 1e-20) return INF; return -(l.a/l.b); } // 返回直线的倾斜角alpha ( 0 - pi) double alpha(LINE l) { if(abs(l.a)< EP) return 0; if(abs(l.b)< EP) return PI/2; double k=slope(l); if(k>0) return atan(k); else return PI+atan(k); } // 求点p关于直线l的对称点 POINT symmetry(LINE l,POINT p) { POINT tp; tp.x=((l.b*l.b-l.a*l.a)*p.x-2*l.a*l.b*p.y-2*l.a*l.c)/(l.a*l.a+l.b*l.b); tp.y=((l.a*l.a-l.b*l.b)*p.y-2*l.a*l.b*p.x-2*l.b*l.c)/(l.a*l.a+l.b*l.b); return tp; } // 如果两条直线 l1(a1*x+b1*y+c1 = 0), l2(a2*x+b2*y+c2 = 0)相交，返回true，且返回交点p bool lineintersect(LINE l1,LINE l2,POINT &p) // 是 L1，L2 { double d=l1.a*l2.b-l2.a*l1.b; if(abs(d)
            
             =0 ) bc[(index+1)%vcount]=1; else bc[(index+1)%vcount]=0; index++; count--; } }// 返回值：多边形polygon是凸多边形时，返回true bool isconvex(int vcount,POINT polygon[]) { bool bc[MAXV]; checkconvex(vcount,polygon,bc); for(int i=0;i
             
              0; } // 另一种判断多边形顶点排列方向的方法 bool isccwize(int vcount,POINT polygon[]) { int i,index; POINT a,b,v; v=polygon[0]; index=0; for(i=1;i
              
               0; } /******************************************************************************************** 射线法判断点q与多边形polygon的位置关系，要求polygon为简单多边形，顶点逆时针排列 如果点在多边形内： 返回0 如果点在多边形边上： 返回1 如果点在多边形外： 返回2 *********************************************************************************************/ int insidepolygon(int vcount,POINT Polygon[],POINT q) { int c=0,i,n; LINESEG l1,l2; bool bintersect_a,bonline1,bonline2,bonline3; double r1,r2; l1.s=q; l1.e=q; l1.e.x=double(INF); n=vcount; for (i=0;i
               
                0) || (bonline3=online(l1,Polygon[(i+2)%n]))&& /* 下一条边是水平线，前一个端点和后一个端点在射线两侧 */ ((r2=multiply(Polygon[i],Polygon[(i+2)%n],l1.s)*multiply(Polygon[(i+2)%n], Polygon[(i+3)%n],l1.s))>0) ) ) ) c++; } if(c%2 == 1) return 0; else return 2; } //点q是凸多边形polygon内时，返回true；注意：多边形polygon一定要是凸多边形 bool InsideConvexPolygon(int vcount,POINT polygon[],POINT q) // 可用于三角形！ { POINT p; LINESEG l; int i; p.x=0;p.y=0; for(i=0;i
                
                 0 || // 极角更小 (multiply(PointSet[j],PointSet[k],PointSet[0])==0) && /* 极角相等，距离更短 */ dist(PointSet[0],PointSet[j])
                 
                  =0) top--; ch[++top]=PointSet[i]; } len=top+1; } // 卷包裹法求点集凸壳，参数说明同graham算法 void ConvexClosure(POINT PointSet[],POINT ch[],int n,int &len) { int top=0,i,index,first; double curmax,curcos,curdis; POINT tmp; LINESEG l1,l2; bool use[MAXV]; tmp=PointSet[0]; index=0; // 选取y最小点，如果多于一个，则选取最左点 for(i=1;i
                  
                   curmax || fabs(curcos-curmax)<1e-6 && dist(l2.s,l2.e)>curdis) { curmax=curcos; index=i; curdis=dist(l2.s,l2.e); } } use[first]=false; //清空第first个顶点标志，使最后能形成封闭的hull use[index]=true; ch[top++]=PointSet[index]; l1.s=ch[top-2]; l1.e=ch[top-1]; l2.s=ch[top-1]; } len=top-1; } /********************************************************************************************* 判断线段是否在简单多边形内(注意：如果多边形是凸多边形，下面的算法可以化简) 必要条件一：线段的两个端点都在多边形内； 必要条件二：线段和多边形的所有边都不内交； 用途： 1. 判断折线是否在简单多边形内 2. 判断简单多边形是否在另一个简单多边形内 **********************************************************************************************/ bool LinesegInsidePolygon(int vcount,POINT polygon[],LINESEG l) { // 判断线端l的端点是否不都在多边形内 if(!insidepolygon(vcount,polygon,l.s)||!insidepolygon(vcount,polygon,l.e)) return false; int top=0,i,j; POINT PointSet[MAXV],tmp; LINESEG s; for(i=0;i
                   
                    PointSet[j].x || fabs(PointSet[i].x-PointSet[j].x)
                    
                     PointSet[j].y ) { tmp=PointSet[i]; PointSet[i]=PointSet[j]; PointSet[j]=tmp; } } } for(i=0;i
                     
                       x1 ) { AddPosPart ((x2+x1)/2, y1/2, (x2-x1) * y1,xtr,ytr,wtr); /* rectangle 全局变量变化处 */ AddPosPart ((x1+x2+x2)/3, (y1+y1+y2)/3, (x2-x1)*(y2-y1)/2,xtr,ytr,wtr); // triangle 全局变量变化处 } else { AddNegPart ((x2+x1)/2, y1/2, (x2-x1) * y1,xtl,ytl,wtl); // rectangle 全局变量变化处 AddNegPart ((x1+x2+x2)/3, (y1+y1+y2)/3, (x2-x1)*(y2-y1)/2,xtl,ytl,wtl); // triangle 全局变量变化处 } } POINT cg_simple(int vcount,POINT polygon[]) { double xtr,ytr,wtr,xtl,ytl,wtl; //注意： 如果把xtr,ytr,wtr,xtl,ytl,wtl改成全局变量后这里要删去 POINT p1,p2,tp; xtr = ytr = wtr = 0.0; xtl = ytl = wtl = 0.0; for(int i=0;i
                      
                       =4的简单多边形至少有一条对角线 结论 ： x坐标最大，最小的点肯定是凸顶点 y坐标最大，最小的点肯定是凸顶点 ************************************************/ POINT a_point_insidepoly(int vcount,POINT polygon[]) { POINT v,a,b,r; int i,index; v=polygon[0]; index=0; for(i=1;i
                       
                        =0; // p is to the left of pre edge bln=multiply(en.e,p,en.s)>=0; // p is to the left of next edge if(!blp&&bln) { if(multiply(polygon[i],rp,p)>0) // polygon[i] is above rp rp=polygon[i]; } if(blp&&!bln) { if(multiply(lp,polygon[i],p)>0) // polygon[i] is below lp lp=polygon[i]; } } return ; } // 如果多边形polygon的核存在，返回true，返回核上的一点p.顶点按逆时针方向输入 bool core_exist(int vcount,POINT polygon[],POINT &p) { int i,j,k; LINESEG l; LINE lineset[MAXV]; for(i=0;i
                        
                         0) //多边形顶点按逆时针方向排列，核肯定在每条边的左侧或边上 break; } if(k == vcount) //找到了一个核上的点 break; } } if(j
                         
                           r1+r2 ) return 1; if( d < fabs(r1-r2) ) return 5; if( fabs(r1-r2) < d && d < r1+r2 ) return 3; return 0; // indicate an error! } //判断圆是否在矩形内：// 判定圆是否在矩形内，是就返回true（设矩形水平，且其四个顶点由左上开始按顺时针排列） // 调用ptoldist函数，在第4页 bool CircleRecRelation(POINT pc, double r, POINT pr1, POINT pr2, POINT pr3, POINT pr4) { if( pr1.x < pc.x && pc.x < pr2.x && pr3.y < pc.y && pc.y < pr2.y ) { LINESEG line1(pr1, pr2); LINESEG line2(pr2, pr3); LINESEG line3(pr3, pr4); LINESEG line4(pr4, pr1); if( r
                          
                            0; }//镜面反射线：// 已知入射线、镜面，求反射线。 // a1,b1,c1为镜面直线方程(a1 x + b1 y + c1 = 0 ,下同)系数; //a2,b2,c2为入射光直线方程系数; //a,b,c为反射光直线方程系数. // 光是有方向的，使用时注意：入射光向量:<-b2,a2>；反射光向量:
                           
                            . // 不要忘记结果中可能会有"negative zeros" void reflect(double a1,double b1,double c1,double a2,double b2,double c2,double &a,double &b,double &c) { double n,m; double tpb,tpa; tpb=b1*b2+a1*a2; tpa=a2*b1-a1*b2; m=(tpb*b1+tpa*a1)/(b1*b1+a1*a1); n=(tpa*b1-tpb*a1)/(b1*b1+a1*a1); if(fabs(a1*b2-a2*b1)<1e-20) { a=a2;b=b2;c=c2; return; } double xx,yy; //(xx,yy)是入射线与镜面的交点。 xx=(b1*c2-b2*c1)/(a1*b2-a2*b1); yy=(a2*c1-a1*c2)/(a1*b2-a2*b1); a=n; b=-m; c=m*yy-xx*n; } //矩形包含： // 矩形2（C，D）是否在1（A，B）内bool r2inr1(double A,double B,double C,double D) { double X,Y,L,K,DMax; if (A < B) { double tmp = A; A = B; B = tmp; } if (C < D) { double tmp = C; C = D; D = tmp; } if (A > C && B > D) // trivial case return true; else if (D >= B) return false; else { X = sqrt(A * A + B * B); // outer rectangle's diagonal Y = sqrt(C * C + D * D); // inner rectangle's diagonal if (Y < B) // check for marginal conditions return true; // the inner rectangle can freely rotate inside else if (Y > X) return false; else { L = (B - sqrt(Y * Y - A * A)) / 2; K = (A - sqrt(Y * Y - B * B)) / 2; DMax = sqrt(L * L + K * K); if (D >= DMax) return false; else return true; } } } //两圆交点： // 两圆已经相交（相切） void c2point(POINT p1,double r1,POINT p2,double r2,POINT &rp1,POINT &rp2) { double a,b,r; a=p2.x-p1.x; b=p2.y-p1.y; r=(a*a+b*b+r1*r1-r2*r2)/2; if(a==0&&b!=0) { rp1.y=rp2.y=r/b; rp1.x=sqrt(r1*r1-rp1.y*rp1.y); rp2.x=-rp1.x; } else if(a!=0&&b==0) { rp1.x=rp2.x=r/a; rp1.y=sqrt(r1*r1-rp1.x*rp2.x); rp2.y=-rp1.y; } else if(a!=0&&b!=0) { double delta; delta=b*b*r*r-(a*a+b*b)*(r*r-r1*r1*a*a); rp1.y=(b*r+sqrt(delta))/(a*a+b*b); rp2.y=(b*r-sqrt(delta))/(a*a+b*b); rp1.x=(r-b*rp1.y)/a; rp2.x=(r-b*rp2.y)/a; } rp1.x+=p1.x; rp1.y+=p1.y; rp2.x+=p1.x; rp2.y+=p1.y; } //两圆公共面积：// 必须保证相交 double c2area(POINT p1,double r1,POINT p2,double r2) { POINT rp1,rp2; c2point(p1,r1,p2,r2,rp1,rp2); if(r1>r2) //保证r2>r1 { swap(p1,p2); swap(r1,r2); } double a,b,rr; a=p1.x-p2.x; b=p1.y-p2.y; rr=sqrt(a*a+b*b); double dx1,dy1,dx2,dy2; double sita1,sita2; dx1=rp1.x-p1.x; dy1=rp1.y-p1.y; dx2=rp2.x-p1.x; dy2=rp2.y-p1.y; sita1=acos((dx1*dx2+dy1*dy2)/r1/r1); dx1=rp1.x-p2.x; dy1=rp1.y-p2.y; dx2=rp2.x-p2.x; dy2=rp2.y-p2.y; sita2=acos((dx1*dx2+dy1*dy2)/r2/r2); double s=0; if(rr
                            
                             0) return -1; else return 1; } /*公式： 球坐标公式： 直角坐标为 P(x, y, z) 时，对应的球坐标是(rsinφcosθ, rsinφsinθ, rcosφ),其中φ是向量OP与Z轴的夹角，范围[0，π]；是OP在XOY面上的投影到X轴的旋角，范围[0，2π] 直线的一般方程转化成向量方程： ax+by+c=0 x-x0 y-y0 ------ = ------- // (x0,y0)为直线上一点，m,n为向量 m n 转换关系： a=n；b=-m；c=m·y0-n·x0； m=-b; n=a; 三点平面方程： 三点为P1，P2，P3 设向量 M1=P2-P1; M2=P3-P1; 平面法向量： M=M1 x M2 （） 平面方程： M.i(x-P1.x)+M.j(y-P1.y)+M.k(z-P1.z)=0