几何算法.cpp

一些常用的几何算法的C++实现,包括两点确定一条直线,两点间距离,点直线距离, 海伦公式求三角形面积,判断直线相交等

#include <iostream>
#include <math.h>

using namespace std;

const double INF=1e100;
const double ZERO=1e-6;
const double PI=2*asin(1.0);

struct XYpoint{  //(x,y)
   double x;
   double y;
};

struct XYline{   // Ax+By+C=0;
   double A;
   double B;
   double C;
};

struct XYsegment{
   XYpoint a,b;
};

struct XYround{  // 圆
   double x;
   double y;
   double r;
};

inline double min(double a,double b)
{
   return a<b?a:b;
}
inline double max(double a,double b)
{
   return a>b?a:b;
}


/********************************************/
/*            两点确定一条直线              */
/********************************************/
XYline makeLine(double x1,double y1,double x2,double y2)
{
   XYline line;
   line.A=(y2-y1);
   line.B=(x1-x2);
   line.C=y1*(x2-x1)+x1*(y1-y2);
   return line;
}

/********************************************/
/*               两点间距离                 */
/********************************************/
double dis_PP(double x1,double y1,double x2,double y2)
{
   return sqrt( (x1-x2)*(x1-x2) + (y1-y2)*(y1-y2) );
}

/********************************************/
/*               点直线距离                 */
/********************************************/
double dis_PL(double x,double y,XYline line)
{
   return fabs(line.A*x + line.B*y + line.C) / sqrt(line.A*line.A + line.B*line.B);
}

/********************************************/
/*          海伦公式求三角形面积            */
/********************************************/
double triangleArea(double x1,double y1,double x2,double y2,double x3,double y3)
{
   double a=sqrt( (x1-x2)*(x1-x2) + (y1-y2)*(y1-y2) );
   double b=sqrt( (x1-x3)*(x1-x3) + (y1-y3)*(y1-y3) );
   double c=sqrt( (x3-x2)*(x3-x2) + (y3-y2)*(y3-y2) );
   double s=(a+b+c)/2;
   return sqrt( s * (s-a) * (s-b) * (s-c) );	
}

/********************************************/
/*              判断直线相交                */
/*    1 - 有交点  0 - 无交点  -1 - 重合     */
/********************************************/
int inter_LL( XYline line1 , XYline line2 , XYpoint &point )
{
   if(line1.A*line2.B == line1.B*line2.A)
   {
      if(line1.A*line2.C == line1.C*line2.A && line1.B*line2.C == line1.C*line2.B)
         return -1;// 重合
      else
         return 0;// 平行
   }
   else
   {
      point.x= -1 * (line1.C*line2.B-line2.C*line1.B) / (line1.A*line2.B-line2.A*line1.B);
      point.y= -1 * (line1.C*line2.A-line2.C*line1.A) / (line1.B*line2.A-line2.B*line1.A);
      // 交点
      return 1;
   }
}

/********************************************/
/*            两直线夹角 （弧度）           */
/********************************************/
double angle_LL(XYline line1 , XYline line2)
{
   double t;
   if( (line1.A*line2.A + line1.B*line2.B)==0 )
      t=INF;
   else
      t=fabs( (line1.A*line2.B - line2.A*line1.B) / (line1.A*line2.A + line1.B*line2.B) );
   return atan(t);
}

/********************************************/
/*          过某点和直线垂直的直线          */
/********************************************/
XYline ver_PL( double x ,double y , XYline line )
{
   XYline temp;
   temp.A=line.B;
   temp.B=-1 * line.A;
   temp.C=temp.A*y - temp.B*x;
   return temp;
}

/********************************************/
/*           多边形面积 ( v[0]=v[n] )       */
/*          n                               */
/*  a=0.5*sigma{ X[i-1]*Y[i]-X[i]*Y[i-1] }  */
/*         i=1                              */
/********************************************/
double polygonArea( int vcount , XYpoint ver[] )
{
   int i;
   double s=0;
   XYpoint t;
   if( vcount<3 )
      return 0;
   t=ver[vcount];
   ver[vcount]=ver[0];
   for(i=1;i<=vcount;i++)
      s+=ver[i-1].x*ver[i].y - ver[i].x*ver[i-1].y ;
   ver[vcount]=t;
   return s/2 ;
}

/********************************************/
/* 返回(P1-P0)*(P2-P0)的叉积。              */
/* 若为正，则<P0,P1>在<P0,P2>的顺时针方向   */
/* 若为零，则<P0,P1><P0,P2>共线；           */
/* 若为负，则<P0,P1>在<P0,P2>的在逆时针方向 */
/* 这个函数可以确定两条线段在交点处的转向   */
/* 比如确定p0p1和p1p2在p1处是左转还是右转， */
/* 只要求 (p2-p0)*(p1-p0)，                 */
/* 若<0则左转，>0则右转，=0则共线           */
/********************************************/
double multiply( XYpoint p1 , XYpoint p2 , XYpoint p0 )
{
	return((p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y));	
}

/********************************************/
/*             判断直线段相交               */
/********************************************/
bool inter_SS( XYsegment u , XYsegment 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));
}


/********************************************/
/*             判断两点是否共点             */
/********************************************/
bool equal_PP(XYpoint a,XYpoint b)
{
   if( fabs(a.x-b.x)<ZERO && fabs(a.y-b.y)<ZERO )
      return true;
   else
      return false;
}

/********************************************/
/*           判断点p是否在线段l上           */
/********************************************/
bool online( XYsegment l , XYpoint p )
{
   if( equal_PP(p,l.a) || equal_PP(p,l.b) )
      return true;
	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 )) );
}

/********************************************/
/*          判断点是否在多边形内            */
/********************************************/
bool insidePolygon( int vcount , XYpoint ver[] , XYpoint point )
{
   int i,c=0;
   XYpoint t;
   XYsegment line,edge;
   t=ver[vcount];
   ver[vcount]=ver[0];
   line.a=line.b=point;
   line.b.x=INF;
   for(i=0;i<vcount;i++)
   {
      edge.a=ver[i]; edge.b=ver[i+1];
      if( online(edge,point) )
         return true;
      if( ver[i].y==ver[i+1].y )
         continue;
      if( min(ver[i].y,ver[i].y)!=point.y && inter_SS(line,edge) )
         c++;
   }
   ver[vcount]=t;
   return c%2 ;
}

/********************************************/
/*             得到线段的中点               */
/********************************************/
XYpoint middle( XYpoint p1 , XYpoint p2 )
{
   XYpoint point;
   point.x= (p1.x+p2.x)/2;
   point.y= (p1.y+p2.y)/2;
   return point;
}

/********************************************/
/*          判断线段是否在多边形内          */
/*    1 - 在内部  0 - 在外部  -1 - 相交     */
/********************************************/
int segmentPolygon( int vcount , XYpoint ver[] , XYsegment seg )
{
   int i;
   XYpoint point,t;
   XYsegment edge;
   t=ver[vcount];
   ver[vcount]=ver[0];
   for(i=0;i<vcount;i++)
   {
      edge.a=ver[i]; edge.b=ver[i+1];
      if( inter_SS(seg,edge) )//&& !online(edge,seg.a) && !online(edge,seg.b) )
         return -1;
   }
   return insidePolygon( vcount , ver , middle(seg.a,seg.b) ) ;
}

/********************************************/
/*      计算线段（p0->p1)的幅角  [0,2*PI)   */
/********************************************/
double angle_PP( XYpoint p0 , XYpoint p1 )
{
   double x,y;
   x=p1.x-p0.x;
   y=p1.y-p0.y;
   if( equal_PP( p0 , p1 ) )
      return -1;
   if( x==0 ) {
      if( y>0 )
         return PI/2;
      else
         return PI+PI/2;
   }
   if( x>0 && y>=0 ) return atan(y/x);
   if( x>0 && y<0 ) return 2*PI+atan(y/x);
   if( x<0 && y>=0 ) return PI+atan(y/x);
   if( x<0 && y<0 ) return PI+atan(y/x);
}

/********************************************/
/*              卷包裹法寻找凸包            */
/********************************************/
int convexHull( int pointNum , XYpoint pointSet[] , XYpoint convex[] )
{
   int i;
   int verNum;
   int lastPoint,nowPoint;
   double lastAngle,nowAngle,tempAngle,t1,t2;
   // 找到 最右 最下的点
   for(i=0,lastPoint=0;i<pointNum;i++)
   {
      if( pointSet[i].y<pointSet[lastPoint].y )
         lastPoint=i;
      else if( pointSet[i].y==pointSet[lastPoint].y &&
               pointSet[i].x<pointSet[lastPoint].x )
         lastPoint=i;
   }
   // 寻找凸包
   lastAngle=PI;
   convex[0]=pointSet[lastPoint];
   verNum=1;
   do {
      nowAngle=lastAngle-2*PI;
      for(i=0;i<pointNum;i++)
      {
         tempAngle=angle_PP( convex[verNum-1] , pointSet[i] );
         if( tempAngle!=-1 )
         {
            t1=lastAngle-tempAngle; t2=lastAngle-nowAngle;
            if(t1<0) t1+=2*PI;
            if(t2<0) t2+=2*PI;
            if( t1<t2 ) {
               nowPoint=i; nowAngle=tempAngle;
            }
            else if( t1==t2 ) {
               if( dis_PP(convex[verNum-1].x,convex[verNum-1].y,pointSet[i].x,pointSet[i].y)
                  <dis_PP(convex[verNum-1].x,convex[verNum-1].y,pointSet[nowPoint].x,pointSet[nowPoint].y) )
                   nowPoint=i; nowAngle=tempAngle;
            }
         }
      }
      lastAngle=nowAngle;
      convex[verNum]=pointSet[nowPoint];
      verNum++;
   }
   while( !equal_PP(convex[0],convex[verNum-1]) ) ;
   return verNum-1;
}

int main()
{
   XYpoint pointSet[1000];
   XYpoint convex[1000];
   int pointNum,verNum,i;
   int cas;
   double len;

   cout.setf(ios::fixed);
   cin>>cas;
   cout<<cas<<endl;
   while(cas--)
   {
      cin>>pointNum;
      if( pointNum==0 )
         break;
      for(i=0;i<pointNum;i++)
         cin>>pointSet[i].x>>pointSet[i].y;
      if(cas) cin>>i;
      verNum=convexHull(pointNum-1,pointSet,convex);
      cout.precision(0);
      cout<<verNum+1<<endl;
      cout<<convex[0].x<<" "<<convex[0].y<<endl;
      for(i=verNum-1;i>=0;i--)
         cout<<convex[i].x<<" "<<convex[i].y<<endl;
      if(cas)
         cout<<"-1"<<endl;
   }
   return 0;
}

/*
int main()
{
   int cas,i,j,r;
   int segNum,plyNum;
   int vcount[100];
   XYpoint ver[100][30];
   XYsegment seg[100];

// input
   cin>>segNum;
   for(i=0;i<segNum;i++)
      cin>>seg[i].a.x>>seg[i].a.y>>seg[i].b.x>>seg[i].b.y;
   cin>>plyNum;
   for(i=0;i<plyNum;i++)
   {
      cin>>vcount[i];
      for(j=0;j<vcount[i];j++)
         cin>>ver[i][j].x>>ver[i][j].y;
   }
// output
   for(i=0;i<segNum;i++)
   {
      for(j=0;j<plyNum;j++)
      {
         r=segmentPolygon(vcount[j],ver[j],seg[i]);
         if(r==-1) {
            cout<<"Submarine "<<i+1<<" is partially on land."<<endl;
            break;
         }
         else if(r==1) {
            cout<<"Submarine "<<i+1<<" is completely on land."<<endl;
            break;
         }
      }
      if(j==plyNum)
         cout<<"Submarine "<<i+1<<" is still in water."<<endl;
   }
   system("pause");
   return 0;
}
*/