ACM经典算法之计算几何

语法：result=polygonarea(Point *polygon,int N);

参数：

*polygon：

多变形顶点数组

N：

多边形顶点数目

返回值：

多边形面积

注意：

支持任意多边形，凹、凸皆可

多边形顶点输入时按顺时针顺序排列

源程序：

typedef struct {
    double x,y;
} Point; 

double polygonarea(Point *polygon,int N)
{
    int i,j;
    double area = 0;

    for (i=0;i<N;i++) {
        j = (i + 1) % N;
        area += polygon[i].x * polygon[j].y;
        area -= polygon[i].y * polygon[j].x;
        }

    area /= 2;
    return(area < 0 ? -area : area);
}

语法：result=area3(float x1,float y1,float x2,float y2,float x3,float y3);

参数：

x1～3：

三角形3个顶点x坐标

y1～3：

三角形3个顶点y坐标

返回值：

三角形面积

注意：

需要 math.h

源程序：

float area3(float x1,float y1,float x2,float y2,float x3,float y3)
{
    float a,b,c,p,s;
    a=sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
    b=sqrt((x1-x3)*(x1-x3)+(y1-y3)*(y1-y3));
    c=sqrt((x3-x2)*(x3-x2)+(y3-y2)*(y3-y2));
    p=(a+b+c)/2;
    s=sqrt(p*(p-a)*(p-b)*(p-c));
    return s;
}

语法：result=angle(double x1, double y1, double x2, double y2);

参数：

x/y1～2：

两矢量的坐标

返回值：

两的角度矢量

注意：

返回角度为弧度制，并且以逆时针方向为正方向

需要 math.h

源程序：

#define PI 3.1415926

double angle(double x1, double y1, double x2, double y2)
{
    double dtheta,theta1,theta2; 
    theta1 = atan2(y1,x1);
    theta2 = atan2(y2,x2);
    dtheta = theta2 - theta1;
    while (dtheta > PI)
        dtheta -= PI*2;
    while (dtheta < -PI)
        dtheta += PI*2; 
    return(dtheta);
}

语法：result=distance_2d(float x1,float x2,float y1,float y2);

参数：

x/y/z1～2：

各点的x、y、z坐标

返回值：

两点之间的距离

注意：

需要 math.h

源程序：

float distance_2d(float x1,float x2,float y1,float y2) 
{
    return(sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)));
}


float distance_3d(float x1,float x2,float y1,float y2,float z1,float z2)
{
    return(sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2)+(z1-z2)*(z1-z2)));
}

语法：result=insidepolygon(Point *polygon,int N,Point p);

参数：

*polygon：

多边形顶点数组

N：

多边形顶点个数

p：

被判断点

返回值：

0：点在多边形内部；1：点在多边形外部

注意：

若p点在多边形顶点或者边上，返回值不确定，需另行判断

需要 math.h

源程序：

#define MIN(x,y) (x < y ? x : y)
#define MAX(x,y) (x > y ? x : y)

typedef struct {
    double x,y;
} Point;

int insidepolygon(Point *polygon,int N,Point p)
{
    int counter = 0;
    int i;
    double xinters;
    Point p1,p2;

    p1 = polygon[0];
    for (i=1;i<=N;i++) {
        p2 = polygon[i % N];
        if (p.y > MIN(p1.y,p2.y)) {
            if (p.y <= MAX(p1.y,p2.y)) {
                if (p.x <= MAX(p1.x,p2.x)) {
                    if (p1.y != p2.y) {
                        xinters = (p.y-p1.y)*(p2.x-p1.x)/(p2.y-p1.y)+p1.x;
                        if (p1.x == p2.x || p.x <= xinters)
                            counter++;
                        }
                    }
                }
            }
            p1 = p2;
        }

    if (counter % 2 == 0)
        return(OUTSIDE);
    else
        return(INSIDE);
}

语法：result=Pointonline(Point p1,Point p2,Point p);

参数：

p1、p2：

线段的两个端点

p：

被判断点

返回值：

0：点在不在线段上；1：点在线段上

注意：

若p线段端点上返回1

需要 math.h

源程序：

#define MIN(x,y) (x < y ? x : y)
#define MAX(x,y) (x > y ? x : y) 

typedef struct {
double x,y;
} Point;

int FC(double x1,double x2)
{
    if (x1-x2<0.000002&&x1-x2>-0.000002) return 1; else return 0;
}


int Pointonline(Point p1,Point p2,Point p)
{
    double x1,y1,x2,y2;
    x1=p.x-p1.x;
    x2=p2.x-p1.x;
    y1=p.y-p1.y;
    y2=p2.y-p1.y;
    if (FC(x1*y2-x2*y1,0)==0) return 0;
    if ((MIN(p1.x,p2.x)<=p.x&&p.x<=MAX(p1.x,p2.x))&&
            (MIN(p1.y,p2.y)<=p.y&&p.y<=MAX(p1.y,p2.y)))
        return 1; else return 0;
}

语法：result=sectintersect(Point p1,Point p2,Point p3,Point p4);

参数：

p1～4：

两条线段的四个端点

返回值：

0：两线段不相交；1：两线段相交；2两线段首尾相接

注意：

p1!=p2;p3!=p4;

源程序：

#define MIN(x,y) (x < y ? x : y)
#define MAX(x,y) (x > y ? x : y) 

typedef struct {
    double x,y;
} Point;

int lineintersect(Point p1,Point p2,Point p3,Point p4)
{
    Point tp1,tp2,tp3;
    if ((p1.x==p3.x&&p1.y==p3.y)||(p1.x==p4.x&&p1.y==p4.y)||(p2.x==p3.x&&p2.y==p3.y)||(p2.x==p4.x&&p2.y==p4.y))
        return 2;
//快速排斥试验
    if ((MIN(p1.x,p2.x)<p3.x&&p3.x<MAX(p1.x,p2.x)&&MIN(p1.y,p2.y)<p3.y<MAX(p1.y,p2.y))||
            (MIN(p1.x,p2.x)<p4.x&&p3.x<MAX(p1.x,p2.x)&&MIN(p1.y,p2.y)<p3.y<MAX(p1.y,p2.y)))
        ;else return 0;

//跨立试验
    tp1.x=p1.x-p3.x;
    tp1.y=p1.y-p3.y;
    tp2.x=p4.x-p3.x;
    tp2.y=p4.y-p3.y;
    tp3.x=p2.x-p3.x;
    tp3.y=p2.y-p3.y;
    if ((tp1.x*tp2.y-tp1.y*tp2.x)*(tp2.x*tp3.y-tp2.y*tp3.x)>=0) return 1; else return 0;
}

语法：result=lineintersect(Point p1,Point p2,Point p3,Point p4);

参数：

p1、p2：

线段的两个端点

p3、p4：

直线上的两个点

返回值：

0：线段直线不相交；1：线段和直线相交

注意：

如线段在直线上，返回 1

源程序：

typedef struct {
    double x,y;
} Point;

int lineintersect(Point p1,Point p2,Point p3,Point p4)
{
    Point tp1,tp2,tp3;
    tp1.x=p1.x-p3.x;
    tp1.y=p1.y-p3.y;
    tp2.x=p4.x-p3.x;
    tp2.y=p4.y-p3.y;
    tp3.x=p2.x-p3.x;
    tp3.y=p2.y-p3.y;
    if ((tp1.x*tp2.y-tp1.y*tp2.x)*(tp2.x*tp3.y-tp2.y*tp3.x)>=0) return 1; else return 0;
}

语法：result=mindistance(Point p1,Point p2,Point q);

参数：

p1、p2：

线段的两个端点

q：

判断点

返回值：

点q到线段p1p2的距离

注意：

需要 math.h

源程序：

#define MIN(x,y) (x < y ? x : y)
#define MAX(x,y) (x > y ? x : y)

typedef struct {
    double x,y;
} Point;

double mindistance(Point p1,Point p2,Point q)
{
    int flag=1;
    double k;
    Point s;
    if (p1.x==p2.x) {s.x=p1.x;s.y=q.y;flag=0;}
    if (p1.y==p2.y) {s.x=q.x;s.y=p1.y;flag=0;}
    if (flag)
        {
        k=(p2.y-p1.y)/(p2.x-p1.x);
        s.x=(k*k*p1.x+k*(q.y-p1.y)+q.x)/(k*k+1);
        s.y=k*(s.x-p1.x)+p1.y;
        }
    if (MIN(p1.x,p2.x)<=s.x&&s.x<=MAX(p1.x,p2.x))
        return sqrt((q.x-s.x)*(q.x-s.x)+(q.y-s.y)*(q.y-s.y));
    else 
        return MIN(sqrt((q.x-p1.x)*(q.x-p1.x)+(q.y-p1.y)*(q.y-p1.y)),sqrt((q.x-p2.x)*(q.x-p2.x)+(q.y-p2.y)*(q.y-p2.y))); }
       

语法：result=mindistance(Point p1,Point p2,Point q);

参数：

p1～p4：

直线上不相同的两点

*p：

通过指针返回结果

返回值：

1：两直线相交；2：两直线平行

注意：

如需要判断两线段交点，检验k和对应k1（注释中）的值是否在0～1之间，用在0～1之间的那个求交点

源程序：

typedef struct {
   double x,y;
} Point;

int linecorss(Point p1,Point p2,Point p3,Point p4,Point *p)
{
   double k;

   //同一直线

  if ((p4.x-p3.x)*(p1.y-p3.y)-(p4.y-p3.y)*(p1.x-p3.x)==0&&

        (p2.x-p1.x)*(p1.y-p3.y)-(p2.y-p1.y)*(p1.x-p3.x)==0) return 2;

   //平行，不同一直线

  if ((p4.y-p3.y)*(p2.x-p1.x)-(p4.x-p3.x)*(p2.y-p1.y)==0) return 0;

    k=((p4.x-p3.x)*(p1.y-p3.y)-(p4.y-p3.y)*(p1.x-p3.x))/((p4.y-p3.y)*(p2.x-p1.x)-(p4.x-p3.x)*(p2.y-p1.y));

//k1=((p2.x-p1.x)*(p1.y-p3.y)-(p2.y-p1.y)*(p1.x-p3.x))/((p4.y-p3.y)*(p2.x-p1.x)-(p4.x-p3.x)*(p2.y-p1.y));

   (*p).x=p1.x+k*(p2.x-p1.x);

   (*p).y=p1.y+k*(p2.y-p1.y);

   return 1;//有交点}

语法：result=convex(Point *p,int n);

参数：

*p：

封闭曲线顶点数组

n：

封闭曲线顶点个数

返回值：

1：凸集；-1：凹集；0：曲线不符合要求无法计算

注意：

默认曲线为简单曲线：无交叉、无圈

源程序：

typedef struct {
    double x,y;
} Point;

int convex(Point *p,int n)
{
    int i,j,k;
    int flag = 0;
    double z;

    if (n < 3)
        return(0);

    for (i=0;i<n;i++) {
        j = (i + 1) % n;
        k = (i + 2) % n;
        z = (p[j].x - p[i].x) * (p[k].y - p[j].y);
        z -= (p[j].y - p[i].y) * (p[k].x - p[j].x);
        if (z < 0)
            flag |= 1;
        else if (z > 0)
            flag |= 2;
        if (flag == 3)
            return －1; //CONCAVE
        }
    if (flag != 0)
        return 1; //CONVEX
    else
    return 0;
}

语法：Graham_scan(Point PointSet[],Point ch[],int n,int &len);

参数：

PointSet[]：

输入的点集

ch[]：

输出的凸包上的点集，按照逆时针方向排列

n：

PointSet中的点的数目

len：

输出的凸包上的点的个数

返回值：

null

源程序：

struct Point{
    float x,y;
}; 

float multiply(Point p1,Point p2,Point p0)
{
    return((p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y)); 
}

float distance(Point p1,Point p2)
{
    return(sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y))); 
}

void Graham_scan(Point PointSet[],Point ch[],int n,int &len)
{
    int i,j,k=0,top=2;
    Point tmp;

   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; 
    for (i=1;i<n-1;i++)
        {
        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;
}