计算几何算法源码【转】

计算几何常用算法，包括凸包算法，平面点集最接远对算法，平面点集最接近对算法，计算任意多边形的面积，三角形相关算法。
                              转自http://blog.csdn.net/ivcan/archive/2008/06/11/2537075.aspx

 

/*
  三角形常用算法模块
*/
const eps=1e-8;
struct TPoint{
   double x,y;
};
struct TLine{
//表示一个直线 a,b,c是参数 a*x+b*y+c=0; 
   double a,b,c;
};
struct TCircle{
//表示圆
   double r;
   TPoint Centre;
};
//三角形描述
typedef TPoint [3] TTriangle //////////////////////////////
bool same(double x,double y)
{ return fabs(x-y)<eps ? 1:0;}
double distance(TPoint p1,TPoint p2)
{ return sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y) )}
TLine lineFromSegment(TPoint p1,TPoint p2)
{//求两点形成的直线
   TPoint tem;
   tem.a=p2.y-p1.y;
   tem.b=p1.x-p2.x;
   tem.c=2*p1.x*p1.y-p1.x*p2.y+p2.x*p1.y;
   return tem;
}
double triangleArea(TTriangle t)
{//三角形面积 
 return abs(t[0].x*t[1].y+t[1].x*t[2].y+t[2].x*t[0].y-t[1].x*t[0].y-t[2].x*t[1].y-t[0].x*t[2].y)/2;
}
TCircle circumcirlceOfTriangle(TTringle t)
{//三角形的外接圆
   TCircle tem;
   double a,b,c,c1,c2;
   double xa,ya.xb,yb,xc,yc;
   a=distance(t[0],t[1]);
   b=distance(t[1],t[2]);
   b=distance(t[0],t[2]);
   tem.r=a*b*c/triangleArea(t)/4;
   xa=t[0].x;ya=t[0].y;
   xb=t[1].x;yb=t[1].y;
   xc=t[2].x;yc=t[2].y;
   c1=(xa*xa+ya*ya-xb*xb-yb*yb)/2;
   c2=(xa*xa+ya*ya-xc*xc-yc*yc)/2;
   tem.centre.x=(c1*(ya-yc)-c2*(ya-yb))/((xa-xb)*(ya-yc)-(xa-xc)*(ya-yb));
   tem.centre.y=(c1*(xa-xc)-c2*(xa-xb))/((ya-yb)*(xa-xc)-(ya-yc)*(xa-xb));
   return tem;
}
TCircle incircleOfTriangle(TTriangle t)
{//三角形的内接圆
   TCircle tem;
   double a,b,c,angleA,angleB,angleC,p,p2,p3,f1,f2;
   double xa,ya,xb,yb,xc,yc;
   a=distance(t[0],t[1]);
   b=distance(t[1],t[2]);
   c=distance(t[2],t[0]);
   tem.r=2*triangleArea(t)/(a+b+c);
   angleA=arccos((b*b+c*c-a*a)/(2*b*c));
   angleB=arccos((a*a+c*c-b*b)/(2*a*c));
   angleC=arccos((b*b+a*a-c*c)/(2*b*c));
   p=sin(angleA/2.0);
   p2=sin(angleB/2.0);
   p3=sin(angleC/2.0);
   xa=t[0].x;ya=t[0].y;
   xb=t[1].x;yb=t[1].y;
   xc=t[2].x;yc=t[2].y;
   f1=((tem.r/p2)*(tem.r/p2)-(tem.r/p)*(tem.r/p)+xa*xa-xb*xb+ya*ya-yb*yb)/2;
   f1=((tem.r/p3)*(tem.r/p3)-(tem.r/p)*(tem.r/p)+xa*xa-xc*xc+ya*ya-yc*yc)/2;
   tem.centre.x=(f1*(ya-yc)-f2*(ya-yb))/((xa-xb)*(ya-yc)-(xa-xc)*(ya-yb));
   tem.centre.y=(f1*(xa-xc)-f2*(xa-xb))/((ya-yb)*(xa-xc)-(ya-yc)*(xa-xb));
   return tem;
}
bool isPointInTriangle(TPoint p,TTriangle t)
{//判断点是否在三角形内 
   TTriangle tem;
   double area;
   int i,j;
   area=0;
   for(i=0;i<=2;i++)
   {
      for(j=0;j<=2;j++)
      {
         if(i==j)tem[j]=p;
         else tem[j]=T[j];
      }
      area+=triangleArea(tem);
   }
   return same(area,triangleArea(t));
}
TPoint symmetricalPointofLine(TPoint p,TLine L)
{//求点p 关于直线 L 的对称点
   TPoint p2;
   double d;
   d=L.a*L.a+L.b*L.b;
   p2.x=(L.b*L.b*p.x-L.a*L.a*p.x-2*L.a*L.b*p.y-2*L.a*L.c)/d;
   p2.y=(L.a*L.a*p.y-L.b*L.b*p.y-2*L.a*L.b*p.x-2*L.b*L.c)/d;
   return p2;
}

 

 

#include <stdio.h>
#include <math.h>
#include <stdlib.h> 
struct POINT
{
   double x,y;
}X[50005];
struct A_POINT
{
   bool operator <= (A_POINT a) const
    {     return (y <= a.y); }
   int p;
   double x,y;
};
int cmp1(const void *a,const void *b)
{
    POINT *c,*d;c=(POINT *)a;d=(POINT *)b;
    return c->x>d->x? 1 :-1 ; 
}
int cmp2(const void *a,const void *b)
{
    A_POINT *c,*d;c=(A_POINT *)a;d=(A_POINT *)b;
    return c->y>d->y? 1 :-1 ; 
}
double dist(POINT a,POINT b)
{
       double dx=a.x-b.x,dy=a.y-b.y;
       return sqrt(dx*dx+dy*dy);
}
double distY(A_POINT a,A_POINT b)
{
       double dx=a.x-b.x,dy=a.y-b.y;
       return sqrt(dx*dx+dy*dy);
}
template <class Type>
void Merge(Type Y[], Type Z[], int l, int m, int r)
{// [l : m] [m + 1 : r]
    Type *a = &Z[l];
    Type *b = &Z[m + 1];
    Type *c = &Y[l];
    int anum = m-l+1, ap = 0;
    int bnum = r-m, bp = 0;
    int cnum = r-l+1, cp = 0; 
 
    while (ap < anum && bp < bnum)
    {
        if (a[ap] <= b[bp])
        {   c[cp++] = a[ap++];}
        else
        {    c[cp++] = b[bp++];}
    }
    while (ap < anum)
    {   c[cp++] = a[ap++]; }
    while (bp < bnum)
   {    c[cp++] = b[bp++];}
} 
void closest(POINT X[], A_POINT Y[], A_POINT Z[], int l, int r, 
             POINT &a, POINT &b, double &d)
{
    if ((r-l)== 1) // two node
    {
        a = X[l]; b = X[r];d = dist(X[l], X[r]);
        return;
    } 
    if ((r-l)== 2)
    {
        double d1 = dist(X[l], X[l + 1]);
        double d2 = dist(X[l + 1], X[r]);
        double d3 = dist(X[l], X[r]);       
        if (d1 <= d2 && d1 <= d3)
        {   a = X[l]; b = X[l + 1]; d = d1; }
        else if (d2 <= d3)
        {   a = X[l + 1]; b = X[r]; d = d2; }
        else
        { a = X[l]; b = X[r]; d = d3;}
        return;
    }
    int mid = (l + r) / 2;
    int f = l;
    int g = mid + 1; 
    int i;
    for (i=l; i<=r; i++)
    {
        if (Y[i].p > mid)
       { Z[g++] = Y[i]; }
        else
        { Z[f++] = Y[i]; }
    } 
    closest(X, Z, Y, l, mid, a, b, d); 
    double dr; POINT ar, br; 
    closest(X, Z, Y, mid+1, r, ar, br, dr);
    if (dr < d)
    {   a = ar; b = br; d = dr;}
    Merge(Y, Z, l, mid, r); // 重构数组Y 
    int k = l;
    for (i=l; i<=r; i++)
    {
        if (fabs(Y[mid].x - Y[i].x) < d)
        {   Z[k++] = Y[i]; }
    } 
    int j;
    for (i=l; i<k; i++)
    {
        for (j=i+1; j<k && Z[j].y-Z[i].y<d; j++)
        {
            double dp = distY(Z[i], Z[j]);
            if (dp < d)
            {   d = dp; a = X[Z[i].p]; b = X[Z[j].p];}
        }
    }
} 
bool closest_Pair(POINT X[], int n, POINT &a, POINT &b,double &d)
{
    if (n < 2) return false; 
    qsort(X,n,sizeof(X[0]),cmp1);
    A_POINT *Y = new A_POINT[n];
    int i;
    for (i=0; i<n; i++)
    {
        Y[i].p = i; // 同一点在数组X中的坐标
        Y[i].x = X[i].x;
        Y[i].y = X[i].y;
    }
    qsort(Y,n,sizeof(Y[0]),cmp2);
    A_POINT *Z = new A_POINT[n];
    closest(X, Y, Z, 0, n - 1, a, b, d);
    delete []Y;
    delete []Z; 
    return true;
} 
int main()
{ 
    int n;
    POINT a, b;
    double d;
    while(1)
    {
       scanf("%d",&n);
       for(int i=0;i<n;i++)
         scanf("%lf%lf",&X[i].x,&X[i].y);
       closest_Pair(X,n,a,b,d);
       printf("%lf",d);
    }
} 

 

 

#include<stdio.h> 
#include<math.h>
#define M 50009
const double INF=1E200; 
const double EP=1E-10; 
#define PI acos(-1) 
/*基本几何数据结构*/ 
//点(x,y) 
struct POINT 
{   int x,y; }; 
// 返回两点之间欧氏距离
double distance(POINT p1, POINT p2) 
{   return sqrt( (p1.x - p2.x) * (p1.x - p2.x) + (p1.y - p2.y) * (p1.y - p2.y)); } 
inline int sqd(POINT a,POINT b)
{//距离平方
      return (a.x-b.x)*(a.x-b.x)+(a.y-b.y)*(a.y-b.y);
}
//叉积就是2向量形成的平行四边形的面积
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)); } 
int partition(POINT a[],int p,int r) 
{ 
   int i=p,j=r+1,k; 
   double ang,dis; 
   POINT R,S; 
   k=(p+r)/2;//防止快排退化
   R=a[p];  a[p]=a[k]; a[k]=R; R=a[p]; 
   dis=distance(R,a[0]); 
   while(1) 
   { 
      while(1) 
      { 
         ++i; 
         if(i>r) 
         {   i=r; break; } 
         ang=multiply(R,a[i],a[0]); 
         if(ang>0) 
            break; 
         else if(ang==0) 
         {   if(distance(a[i],a[0])>dis)     break; } 
      } 
      while(1) 
      {   --j; 
         if(j<p) 
         {   j=p;   break; } 
         ang=multiply(R,a[j],a[0]); 
         if(ang<0)   break; 
         else if(ang==0) 
         {   if(distance(a[j],a[0])<dis)    break; } 
      } 
      if(i>=j)break; 
      S=a[i]; a[i]=a[j]; a[j]=S; 
   } 
a[p]=a[j]; a[j]=R; return j; 
} 
void anglesort(POINT a[],int p,int r) 
{ 
   if(p<r) 
   { 
      int q=partition(a,p,r); 
      anglesort(a,p,q-1); 
      anglesort(a,q+1,r); 
   } 
} 
void Graham_scan(POINT PointSet[],POINT ch[],int n,int &len) 
{ 
      int i,k=0,top=2; 
      POINT 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] 
      /* 对顶点按照相对PointSet[0]的极角从小到大进行排序，极角相同
                  的按照距离PointSet[0]从近到远进行排序 */ 
      anglesort(PointSet,1,n-1); 
 
               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; 
} 
int main() 
{ 
   POINT a[M],b[M]; 
   int n,i,l; 
   double x,y;
   scanf("%d",&n);
   for(i=0;i<n;i++)scanf("%d%d",&a[i].x,&a[i].y);
   Graham_scan(a,b,n,l); 
   int max=0;
      for(int i=0;i<l;i++)for(int j=i+1;j<l;j++)
      {
           int tmp=sqd(b[i],b[j]);
           if(max<tmp)max=tmp;
      }
      printf("%d\n",max);
   //while(1);
return 0; 
}