计算几何算法 java_计算几何算法总集 - Feng.Li's Java See - BlogJava

计算几何算法总集

#include

#include

struct Point{

double x,y;

};

int dblcmp(double d)

{

if(fabs(d)<0.000000001)return 0;

return (d>0)?1:-1;

}

double det(double x1,double y1,double x2,double y2)

{ return x1*y2-x2*y1;}

double cross(Point a,Point b,Point c)

{//叉积

return det(b.x-a.x, b.y-a.y, c.x-a.x, c.y-a.y);

}

int xycmp(double p,double mini,double maxi)

{ return dblcmp(p-mini)*dblcmp(p-maxi);}

double min(double a,double b)

{ if(a

double max(double a,double b)

{   if(a>b)return a; else return b;}

int betweencmp(Point a,Point b,Point c )

{//点a与b或c重合时为0

//点a在bc内部时为-1

//点a在bc外部时为1

if( fabs(b.x-c.x)>fabs(b.y-c.y))

return xycmp(a.x,min(b.x,c.x),max(b.x,c.x));

else

return xycmp(a.x,min(b.y,c.y),max(b.y,c.y));

}

int segscross(Point a,Point b,Point c,Point d,Point &p)

{//线段ab ,cd规范相交返回1 ,并求交点P ;不规范相交返回2 ;没有交点返回0

double s1,s2,s3,s4;

int d1,d2,d3,d4;

d1=dblcmp(s1=cross(a,b,c));

d2=dblcmp(s2=cross(a,b,d));

d3=dblcmp(s3=cross(c,d,a));

d4=dblcmp(s4=cross(c,d,b));

if( ((d1^d2)==-2) && ((d3^d4)==-2))

{

p.x=(c.x*s2-d.x*s1)/(s2-s1);

p.y=(c.y*s2-d.y*s1)/(s2-s1);

return 1;

}

if( ((d1==0)&& (betweencmp(c,a,b)<=0)) ||

((d2==0)&& (betweencmp(d,a,b)<=0)) ||

((d3==0)&& (betweencmp(a,c,d)<=0)) ||

((d4==0)&& (betweencmp(b,c,d)<=0)) )

return 2;

return 0;

}

double area(Point a,Point b)

{ return a.x*b.y-a.y*b.x;}

double areas(Point A[],int n)

{//求n个点的面积A中的点按逆时针方向存放

//多边形是任意的凸或凹多边形，

double re=0;

int i;

if(n<3)return 0;

for(i=0;i

re+=area(A[i],A[(i+1)%n]);

re=fabs(re/2);

}

void MidPoint(Point A[],int n,Point &p)

{//求多边形的重心A中的点按逆时针方向存放

int i;

double areass,px=0,py=0,tem;

areass=areas(A, n);

for(i=0;i

{ tem=area(A[i],A[(i+1)%n]);

px+=tem*(A[i].x+A[(i+1)%n].x);

py+=tem*(A[i].y+A[(i+1)%n].y);

}

p.x=fabs(px)/(6*areass);

p.y=fabs(py)/(6*areass);

}

int main()

{

Point a,b,c,d,p;

Point   A[10000];

int n;

a.x=0;a.y=0;

b.x=1;b.y=1;

c.x=0;c.y=2;

d.x=3;d.y=0;

int i=segscross(a,b,c,d,p);

printf("%d"n%lf %lf"n",i,p.x,p.y);

while(1);

}

//另一版本

#include

#define infinity 1e20

#define EP 1e-10

struct TPoint{

float x,y;

};

struct TLineSeg{

TPoint a,b;

};

//求平面上两点之间的距离

float distance(TPoint p1,TPoint p2)

{

return(sqrt((p1.x-p2.x)*(p1.x-p2.x)+(p1.y-p2.y)*(p1.y-p2.y)));

}

/********************************************************

返回(P1-P0)*(P2-P0)的叉积。

若结果为正，则在的顺时针方向；

若为0则共线；

若为负则在的在逆时针方向;

可以根据这个函数确定两条线段在交点处的转向,

比如确定p0p1和p1p2在p1处是左转还是右转，只要求

(p2-p0)*(p1-p0)，若<0则左转，>0则右转，=0则共线

*********************************************************/

float multiply(TPoint p1,TPoint p2,TPoint p0)

{

return((p1.x-p0.x)*(p2.y-p0.y)-(p2.x-p0.x)*(p1.y-p0.y));

}

//确定两条线段是否相交

int intersect(TLineSeg u,TLineSeg 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));

}

//判断点p是否在线段l上

int online(TLineSeg l,TPoint p)

{

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 )) );

}

//判断两个点是否相等

int Euqal_Point(TPoint p1,TPoint p2)

{

return((abs(p1.x-p2.x)

}

//一种线段相交判断函数，当且仅当u,v相交并且交点不是u,v的端点时函数为true;

int intersect_A(TLineSeg u,TLineSeg v)

{

return((intersect(u,v))&&

(!Euqal_Point(u.a,v.a))&&

(!Euqal_Point(u.a,v.b))&&

(!Euqal_Point(u.b,v.a))&&

(!Euqal_Point(u.b,v.b)));

}

/*===============================================

判断点q是否在多边形Polygon内，

其中多边形是任意的凸或凹多边形，

Polygon中存放多边形的逆时针顶点序列

================================================*/

int InsidePolygon(int vcount,TPoint Polygon[],TPoint q)

{

int c=0,i,n;

TLineSeg l1,l2;

l1.a=q;

l1.b=q;

l1.b.x=infinity;

n=vcount;

for (i=0;i

{

l2.a=Polygon[i];

l2.b=Polygon[(i+1)%n];

if ( (intersect_A(l1,l2))||

(

(online(l1,Polygon[(i+1)%n]))&&

(

(!online(l1,Polygon[(i+2)%n]))&&

(multiply(Polygon[i],Polygon[(i+1)%n],l1.a)*multiply(Polygon[(i+1)%n],Polygon[(i+2)%n],l1.a)>0)

||

(online(l1,Polygon[(i+2)%n]))&&

(multiply(Polygon[i],Polygon[(i+2)%n],l1.a)*multiply(Polygon[(i+2)%n],Polygon[(i+3)%n],l1.a)>0)

)

)

) c++;

}

return(c%2!=0);

}

/********************************************"

*计算多边形的面积*

*                                            *

*要求按照逆时针方向输入多边形顶点*

*可以是凸多边形或凹多边形*

"********************************************/

float area_of_polygon(int vcount,float x[],float y[])

{

int i;

float s;

if (vcount<3) return 0;

s=y[0]*(x[vcount-1]-x[1]);

for (i=1;i

s+=y[i]*(x[(i-1)]-x[(i+1)%vcount]);

return s/2;

}

void Graham_scan(TPoint PointSet[],TPoint ch[],int n,int &len)

{//寻找凸包graham扫描法

int i,j,k=0,top=2;

TPoint tmp;

//选取PointSet中y坐标最小的点PointSet[k]，如果这样的点右多个，则取最左边的一个

for(i=1;i

if ((PointSet[i].y

k=i;

tmp=PointSet[0];

PointSet[0]=PointSet[k];

PointSet[k]=tmp; //现在PointSet中y坐标最小的点在PointSet[0]

//对顶点按照相对PointSet[0]的极角从小到大进行排序，极角相同的按照距离PointSet[0]从近到远进行排序

for (i=1;i

{     k=i;

for (j=i+1;j

if ( (multiply(PointSet[j],PointSet[k],PointSet[0])>0)

||

( (multiply(PointSet[j],PointSet[k],PointSet[0])==0)

&&(distance(PointSet[0],PointSet[j])

) 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

{

while (multiply(PointSet[i],ch[top],ch[top-1])>=0) top--;

ch[++top]=PointSet[i];

}

len=top+1;

}

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)

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

#include

#include

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

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

{

for (j=i+1; 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

{

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

scanf("%lf%lf",&X[i].x,&X[i].y);

closest_Pair(X,n,a,b,d);

printf("%lf",d);

}

}

平面点集最接远对算法(完全正确)

#include

#include

#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

{   j=p;   break; }

ang=multiply(R,a[j],a[0]);

if(ang<0)   break;

else if(ang==0)

{   if(distance(a[j],a[0])

}

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

{

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

if ( PointSet[i].y

(PointSet[i].y==PointSet[k].y) && (PointSet[i].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

{

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

Graham_scan(a,b,n,l);

int max=0;

for(int i=0;i

{

int tmp=sqd(b[i],b[j]);

if(max

}

printf("%d"n",max);

//while(1);

return 0;

}

凸包算法(nlogn)

#include

#include

#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

{   j=p; break;   }

ang=multiply(R,a[j],a[0]);

if(ang<0) break;

else if(ang==0)

{ if(distance(a[j],a[0])

}

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

{

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

if ( PointSet[i].y

(PointSet[i].y==PointSet[k].y) && (PointSet[i].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

{

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

Graham_scan(a,b,n,l);

int max=0;

for(int i=0;i

for(int j=i+1;j

{

int tmp=sqd(b[i],b[j]);

if(max

}

printf("%d"n",max);

//while(1);

return 0;

}