计算几何算法源码【转】
计算几何常用算法(包括平面点集最接远对算法,平面点集最接近对算法,三角形相关算法)
计算几何常用算法,包括凸包算法,平面点集最接远对算法,平面点集最接近对算法,计算任意多边形的面积,三角形相关算法。
转自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;
}
/*
三角形常用算法模块
*/
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>
#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;
}
#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;
}