計算幾何算法總集
#include<stdio.h>
#include<math.h>
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<b)return a; else return b;}
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<n;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<n;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 <stdlib.h>
#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)的叉積。
若結果為正，則<P0,P1>在<P0,P2>的順時針方向；
若為0則<P0,P1><P0,P2>共線；
若為負則<P0,P1>在<P0,P2>的在逆時針方向;
可以根據這個函數確定兩條線段在交點處的轉向,
比如確定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)<EP)&&(abs(p1.y-p2.y)<EP));
}
//一種線段相交判斷函數，當且僅當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<vcount;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<vcount;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<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]從近到遠進行排序
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;
}
////////////////////////////////////////////////////////////////////////////
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;
}
凸包算法(nlogn)
#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;
}