计算几何的相关知识

我有预感以后还会碰到关于平面几何或者立体几何题的编程实现，所以要未雨绸缪。。。以下转载一篇关于这方面知识的好帖子，其实没有细看，因为很多东西还没用到，可是刚好碰到了就不想放它走。。。

给出文章的出处 计算几何所用知识回收站

以后可以慢慢补充

1.向量的叉积

向量a=(x1,y1),b=(x2,y2)；

向量的叉积|a×b|=x1*y2-x2*y1;

(1)当|a×b|>0时，b在a的逆时针方向，当a×b=0时，b与a共线，当a×b<0时，b在a的顺时针方向。

(2)|a×b|=|a|
|b|sin<a,b>，三角形ABC，边向量a=AB，b=AC，面积S=1/2*|a×b|=1/2*(x1*y2-x2*y1)

（1）向量积判断线段相交

[cpp] view
plaincopyprint?





struct node{

double x1,y1,x2,y2;

}e[maxn],f[maxn];

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

{

return x1*y2-x2*y1;

}

int find(node a,node b) //判断是否相交

{

double c[4];

c[0]=cross(a.x2-a.x1,a.y2-a.y1,b.x1-a.x1,b.y1-a.y1);

c[1]=cross(a.x2-a.x1,a.y2-a.y1,b.x2-a.x1,b.y2-a.y1);

c[2]=cross(b.x2-b.x1,b.y2-b.y1,a.x1-b.x1,a.y1-b.y1);

c[3]=cross(b.x2-b.x1,b.y2-b.y1,a.x2-b.x1,a.y2-b.y1);

if(c[0]*c[1]<=0&&c[2]*c[3]<=0)return 1;

return 0;

}

感觉原著里的代码写的不如我在实际解决一道题时候想的周全啊，而且只有生硬的代码让和我一样数学不好的朋友看起来会很头疼吧，这里再搬点好料过来！！！

第1 步：快速排斥试验，如果分别以P1P2 ，P3P4 为对角线做矩形，而这两个矩形不相交，则这两条线段肯定不相交，如下左图；即使两个矩形相交，这两条线段也不一定相交，如下右图，这时再用第2 步判断；



表示成语句，即两个矩形相交当且仅当下列式子为真：

（max(x1,x2)≥min(x3,x4)）∧ (max(x3,x4)≥min(x1,x2)) ∧（max(y1,y2)≥min(y3,y4)）∧(max(y3,y4)≥min(y1,y2))

两个矩形相交必须在两个方向上都相交，式子的前半部分判断在x 方向上是否相交，后半部分判断在y 方向上是否相交。

第2 步：确定每条线段是否“跨立”另一条线段所在的直线。

跨立：如果点P1 处于直线P3P4的一边，而P2处于该直线的另一边，则我们说线段p1p2跨立直线P3P4，如果P1 或P2 在直线P3P4 上，也算跨立。

两条线段相交当且仅当它们能够通过第1 步的快速排斥试验，并且每一条线段都跨立另一条线段所在的直线。



具体第2 步的实现，只要用叉积去做就可以了，即只要判断矢量p1p3和p1p4是否在p1p2的两边相对的位置上，如果这样，则线段p1p2跨立直线P3P4。也即检查叉积（P3-P1）×（P2-P1）与（P4-P1）×（P2-P1）的符号是否相同，相同则不跨立，线段也就不相交，否则相交。

当然也有一些特殊情况需要处理，如任何一个叉积为0，则P3 或P4 在直线P1P2 上，又因为通过了快速排斥试验，所以下图左边的情况是不可能出现的，只会出现右边的两种情况。当然，还会出现一条或两条线段的长度为0，如果两条线段的长度都是0，则只要通过快速排斥试验就能确定；如果仅有一条线段的长度为0，如p3p4的长度为0，则线段相交当且仅当叉积（P3-P1）×（P2-P1）。



还有一些情况要考虑：





然后给出我学习之后的代码

#define eps 1e-6
struct Coord{
double x,y;
};

struct Segment{
Coord a,b;
}str[101];

bool Judge(Coord &a,Coord &b,Coord &c,Coord &d){
if(!(min(a.x,b.x)<=max(c.x,d.x)&&max(a.x,b.x)>=min(c.x,d.x)
&&min(a.y,b.y)<=max(c.y,d.y)&&max(a.y,b.y)>=min(c.y,d.y))){
return false;
}
double u,v,w,z;
u = (c.x-a.x)*(b.y-a.y)-(b.x-a.x)*(c.y-a.y);
v = (d.x-a.x)*(b.y-a.y)-(b.x-a.x)*(d.y-a.y);
w = (a.x-c.x)*(d.y-c.y)-(d.x-c.x)*(a.y-c.y);
z = (b.x-c.x)*(d.y-c.y)-(d.x-c.x)*(b.y-c.y);
return (u*v<=eps&&w*z<=eps);
}

2.重心

(1)三角形重心：x=(x0+x1+x3)/3，y=(y0+y1+y2)/3;

(2)质量集中在顶点上。n个顶点坐标为(xi,yi)，质量为mi，则重心

　　X = ∑( xi×mi ) / ∑mi

　　Y = ∑( yi×mi ) / ∑mi

　　特殊地，若每个点的质量相同，则

　　X = ∑xi / n

　　Y = ∑yi / n

(3)已知一多边形没有边相交，质量分布均匀。顺序给出多边形的顶点坐标，求其重心。

将n+2多变形分成n个三角形，总面积为S，分重心为(xi,yi),分面积为si，则重心为X=(∑si*xi)/S,Y=(∑si*yi)/S;(凹多边形和凸多边形，均适用)

这里的si=0.5*(a×b),有正负值，表示边是按顺时针还是逆时针，用于求凹多边形时加减多余部分。

3.平面几何预算模板

[cpp] view
plaincopyprint?





#include <cstdio>

#include <cstdlib>

#include <cmath>

#include <cstring>

#include <iostream>

#include <algorithm>

#include <map>

#include <set>

#include <queue>

using namespace std;

//基础点和向量运算

struct Point{

double x,y;

Point(double x=0,double y=0):x(x),y(y){}

};

typedef Point Vector;

Vector operator + (Vector A,Vector B){return Vector(A.x+B.x,A.y+B.y);}

Vector operator - (Vector A,Vector B){return Vector(A.x-B.x,A.y-B.y);}

Vector operator * (Vector A,double p){return Vector(A.x*p,A.y*p);}

Vector operator / (Vector A,double p){return Vector(A.x/p,A.y/p);}

bool operator <(const Point& a, const Point& b)

{

return a.x<b.x||(a.x==b.x&&a.y<b.y);

}

const double eps=1e-10;

int dcmp(double x)//判断正负，或者等于0

{

if(fabs(x)<eps)return 0;else return x<0?-1:1;

}

bool operator==(const Point& a,const Point &b)

{

return dcmp(a.x-b.x)==0&&dcmp(a.y-b.y)==0;

}

double Dot(Vector A, Vector B){return A.x*B.x+A.y*B.y;}//点积

double Length(Vector A){return sqrt(Dot(A,A));}//OA长

double Angle(Vector A,Vector B){return acos(Dot(A,B)/Length(A)/Length(B));}//OA和OB的夹角

double Cross(Vector A,Vector B){return A.x*B.y-A.y*B.x;}//叉积

double Area2(Point A,Point B,Point C){return Cross(B-A,C-A);}//三角形面积

Vector Rotate(Vector A,double rad)//rad为弧度，旋转rad度

{

return Vector(A.x*cos(rad)-A.y*sin(rad),A.x*sin(rad)+A.y*cos(rad));

}

Vector Normal(Vector A)//A的单位法向量，A不能为零向量

{

double L=Length(A);

return Vector(-A.y/L,A.x/L);

}

//点和直线

//P+tv表示一条直线，P为点，tv为方向向量

Point GetLineIntersection(Point P,Vector v,Point Q,Vector w)//求直线交点，确保存在交点，即Cross(v,w)非0

{

Vector u=P-Q;

double t=Cross(w,u)/Cross(v,w);

return P+v*t;

}

double DistanceToLine(Point P,Point A,Point B)//P点到直线AB的距离

{

Vector v1=B-A,v2=P-A;

return fabs(Cross(v1,v2)/Length(v1));

}

double DistanceToSegment(Point P,Point A,Point B)//点P到线段AB的距离

{

if(A==B)return Length(P-A);

Vector v1=B-A,v2=P-A,v3=P-B;

if(dcmp(Dot(v1,v2))<0)return Length(v2);

else if(dcmp(Dot(v1,v3))>0)return Length(v3);

else return fabs(Cross(v1,v2)/Length(v1));

}

Point GetLineProjection(Point P,Point A,Point B)//点在直线上的投影

{

Vector v=B-A;

return A+v*(Dot(v,P-A)/Dot(v,v));

}

bool SegmentProperIntersection(Point a1,Point a2,Point b1,Point b2)//判断线段相交,不在端点相交

{

double c1=Cross(a2-a1,b1-a1),c2=Cross(a2-a1,b2-a1),c3=Cross(b2-b1,a1-b1),c4=Cross(b2-b1,a2-b1);

return dcmp(c1)*dcmp(c2)<0&&dcmp(c3)*dcmp(c4)<0;

}

bool OnSegment(Point p,Point a1,Point a2)//判断点是否在线段上（不包括端点）

{

return dcmp(Cross(a1-p,a2-p))==0&&dcmp(Dot(a1-p,a2-p))<0;

}

//多边型

double ConvexPolygonArea(Point* p,int n)//多边形面积,，点按顺序

{

double area=0;

for(int i=1;i<n-1;i++)

area+=Cross(p[i]-p[0],p[i+1]-p[0]);

return area/2;

}

int main()

{

}

圆：

[cpp] view
plaincopyprint?





const double PI=atan(1.0)*4;

//圆

struct Circle{//定义圆

Point c;

double r;

Circle(Point c,double r):c(c),r(r){}

Point point(double a){//根据圆心角计算圆上的坐标

return Point(c.x+cos(a)*r,c.y+sin(a)*r);

}

};

struct Line{//定义线

Point p,v;

Line(Point p,Point v):p(p),v(v){}

Point point(double a){

return p+(v-p)*a;

}

};

int getLineCircleIntersection(Line L,Circle C,double &t1,double &t2,vector<Point> &sol)//圆与线相交

{

double a=L.v.x,b=L.p.x-C.c.x,c=L.v.y,d=L.p.y-C.c.y;

double e=a*a+c*c,f=2*(a*b+c*d),g=b*b+d*d-C.r*C.r;

double delta=f*f-4*e*g;

if(dcmp(delta)<0)return 0;

if(dcmp(delta)==0)

{

t1=t2=-f/(2*e);

sol.push_back(L.point(t1));

return 1;

}

t1=(-f-sqrt(delta))/(2*e);sol.push_back(L.point(t1));

t2=(-f+sqrt(delta))/(2*e);sol.push_back(L.point(t2));

return 2;

}

double angle(Vector v){return atan2(v.y,v.x);}//计算向量极角

int getCircleCircleIntersection(Circle C1,Circle C2,vector<Point> &sol)//圆与圆相交

{

double d=Length(C1.c-C2.c);

if(dcmp(d)==0)

{

if(dcmp(C1.r-C2.r)==0)return -1;//两圆重合

return 0;

}

if(dcmp(C1.r+C2.r-d)<0)return 0;

if(dcmp(fabs(C1.r-C2.r)-d)>0)return 0;

double a=angle(C2.c-C1.c); //向量C1C2的极角

double da=acos((C1.r*C1.r+d*d-C2.r*C2.r)/(2*C1.r*d));

//C1C2到C1P1的角

Point p1=C1.point(a-da),p2=C1.point(a+da);

sol.push_back(p1);

if(p1==p2)return 1;

sol.push_back(p2);

return 2;

}

int getTangents(Point p,Circle C,Vector *v)//过定点作圆的切线

{

Vector u=C.c-p;

double dist=Length(u);

if(dist<C.r)return 0;

else if(dcmp(dist-C.r)==0)

{

v[0]=Rotate(u,PI/2);

return 1;

}

else

{

double ang=asin(C.r/dist);

v[0]=Rotate(u,-ang);

v[1]=Rotate(u,+ang);

return 2;

}

}

//求两圆的切线

//返回切线的条数，-1表示无穷条切线

//a[i]和b[i]分别是第i条切线在圆A和圆B上的切点

int getTangents(Circle A,Circle B,Point *a,Point *b)

{

int cnt=0;

if(A.r<B.r){swap(A,B);swap(a,b);}

double d2=(A.c.x-B.c.x)*(A.c.x-B.c.x)+(A.c.y-B.c.y)*(A.c.y-B.c.y);

double rdiff=A.r-B.r;

double rsum=A.r+B.r;

if(dcmp(d2-rdiff*rdiff)<0)return 0;//内含

double base=atan2(B.c.y-A.c.y,B.c.x-A.c.x);

if(d2==0&&A.r==B.r)return -1;//重合

if(dcmp(d2-rdiff*rdiff)==0) //内切

{

a[cnt]==A.point(base);b[cnt]=B.point(base);cnt++;

return 1;

}

//外公切线

double ang=acos((A.r-B.r)/sqrt(d2));

a[cnt]=A.point(base+ang);b[cnt]=B.point(base+ang);cnt++;

a[cnt]=A.point(base-ang);b[cnt]=B.point(base-ang);cnt++;

if(dcmp(d2-rsum*rsum)==0)

{

a[cnt]=A.point(base);b[cnt]=B.point(PI+base);cnt++;

}

else if(dcmp(d2-rsum*rsum)>0)

{

double ang=acos(A.r+B.r)/sqrt(d2);

a[cnt]=A.point(base+ang);b[cnt]=B.point(PI+base+ang);cnt++;

a[cnt]=A.point(base-ang);b[cnt]=B.point(PI+base-ang);cnt++;

}

return cnt;

}

//三角形外接圆

Circle CircumscribedCircle(Point p1,Point p2,Point p3)

{

double Bx=p2.x-p1.x,By=p2.y-p1.y;

double Cx=p3.x-p1.x,Cy=p3.y-p1.y;

double D=2*(Bx*Cy-By*Cx);

double cx=(Cy*(Bx*Bx+By*By)-By*(Cx+Cx+Cy*Cy))/D+p1.x;

double cy=(Bx*(Cx*Cx+Cy*Cy)-Cx*(Bx*Bx+By*By))/D+p1.y;

Point p=Point(cx,cy);

return Circle(p,Length(p1-p));

}

//三角形内切圆

Circle InscribedCircle(Point p1,Point p2,Point p3)

{

double a=Length(p2-p3);

double b=Length(p3-p1);

double c=Length(p1-p2);

Point p=(p1*a+p2*b+p3*c)/(a+b+c);

return Circle(p,DistanceToLine(p,p1,p2));

}

int ConvexHull(Point *p,Point *ch,int n)//求凸包

{

sort(p,p+n);

int i,m=0,k;

for(i=0;i<n;i++)

{

while(m>1&&Cross(ch[m-1]-ch[m-2],p[i]-ch[m-2])<=0)m--;

ch[m++]=p[i];

}

k=m;

for(i=n-2;i>=0;i--)

{

while(m>k&&Cross(ch[m-1]-ch[m-2],p[i]-ch[m-2])<=0)m--;

ch[m++]=p[i];

}

if(n>1)m--;

return m;

}

三维几何:

[cpp] view
plaincopyprint?





#include <cstdio>

#include <cstdlib>

#include <cmath>

#include <cstring>

#include <iostream>

#include <algorithm>

#include <map>

#include <set>

#include <queue>

using namespace std;

const double pi=atan(1.0)*4;

const double r=6371009;

const double eps=1e-10;

struct Point3{//定义三维点

double x,y,z;

Point3(double x=0,double y=0,double z=0):x(x),y(y),z(z){}

};

typedef Point3 Vector3;

Vector3 operator +(Vector3 A,Vector3 B){return Vector3(A.x+B.x,A.y+B.y,A.z+B.z);}

Vector3 operator -(Vector3 A,Vector3 B){return Vector3(A.x-B.x,A.y-B.y,A.z-B.z);}

Vector3 operator *(Vector3 A,double p){return Vector3(A.x*p,A.y*p,A.z*p);}

Vector3 operator /(Vector3 A,double p){return Vector3(A.x/p,A.y/p,A.z/p);}

double Dot(Vector3 A,Vector3 B){return A.x*B.x+A.y*B.y+A.z*B.z;}//点积

double Length(Vector3 A){return sqrt(Dot(A,A));}

double Angle(Vector3 A,Vector3 B){return acos(Dot(A,B)/Length(A)/Length(B));}

int dcmp(double x)

{

if(fabs(x)<eps)return 0;

else return x<0?-1:1;

}

bool operator ==(Vector3 A,Vector3 B){return dcmp(A.x-B.x)==0&&dcmp(A.y-B.y)==0&&dcmp(A.z-B.z);}

//点p到平面p0-n的距离。n必须为单位向量

double DistanceToPlane(const Point3 &p,const Point3 &p0,const Vector3 &n)

{

return fabs(Dot(p-p0,n));

}

//点在平面上的投影。n必须是单位向量

Point3 GetPlaneProjection(const Point3 &p,const Point3 &p0,const Vector3 &n)

{

return p-n*Dot(p-p0,n);

}

//直线p1-p2和平面的交点。假定点唯一存在

Point3 LinePlaneIntersection(Point3 p1,Point3 p2,Point3 p0,Vector3 n)

{

Vector3 v=p2-p1;

double t=(Dot(n,p0-p1)/Dot(n,p2-p1));//判断分母是否为0

return p1+v*t;//如果是线段，判断t是否在0和1之间。

}

//三维叉积

Vector3 Cross(Vector3 A,Vector3 B)

{

return Vector3(A.y*B.z-A.z*B.y,A.z*B.x-A.x*B.z,A.x*B.y-A.y*B.x);

}

double Area2(Point3 A,Point3 B,Point3 C){return Length(Cross(B-A,C-A));}

//判断点P0是否在三角形△P0P1P2中

bool PointInTri(Point3 P,Point3 P0,Point3 P1,Point3 P2)

{

double area1=Area2(P,P0,P1);

double area2=Area2(P,P1,P2);

double area3=Area2(P,P2,P0);

return dcmp(area1+area2+area3-Area2(P0,P1,P2))==0;

}

//△P0P1P2是否和线段AB相交

bool TriSegIntersection(Point3 P0,Point3 P1,Point3 P2,Point3 A,Point3 B,Point3 &P)

{

Vector3 n=Cross(P1-P0,P2-P0);

if(dcmp(Dot(n,B-A))==0)return false;//线段AB和平面P0P1P2平行或者共面

else

{

double t=Dot(n,P0-A)/Dot(n,B-A);

if(dcmp(t)<0||dcmp(t-1)>0)return false;//交点不在线段AB上

P=A+(B-A)*t; //计算交点

return PointInTri(P,P0,P1,P2);//判断交点是否在三角形内

}

}

//点P到直线AB的距离

double DistanceToLine(Point3 P,Point3 A,Point3 B)

{

Vector3 v1=B-A,v2=P-A;

return Length(Cross(v1,v2))/Length(v1);

}

//点P到线段AB的距离

double DistanceToSegment(Point3 P,Point3 A,Point3 B)

{

if(A==B)return Length(P-A);

Vector3 v1=B-A,v2=P-A,v3=P-B;

if(dcmp(Dot(v1,v2))<0)return Length(v2);

else if(dcmp(Dot(v1,v3))>0)return Length(v3);

else return Length(Cross(v1,v2))/Length(v1);

}

//返回AB,AC,AD的混合积。它也等于四面体ABCD的有向体积的6倍。

double Volume6(Point3 A,Point3 B,Point3 C,Point3 D)

{

return Dot(D-A,Cross(B-A,C-A));

}

int main()

{

}