几何计算方法-优快云博客

初始定义：（基本都要用到，做题直接全打上）

struct Point
{
  double x,y;
  Point(double x=0,double y=0):x(x),y(y){} //构造函数，以后直接用Point(x,y)就表示点
};
struct Polygon
{
  double x,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)
{
  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;
}
*/
//有精度要求时采用上面的写法
bool operator == (const Point& a, const Point& b){
    return a.x == b.x && a.y == b.y;
}
struct Line 
{
  Point P;
  Vector v;
  double ang;
  Line(){}
  Line(Point P,Vector v):P(P),v(v){ ang=atan2(v.y,v.x);}
  bool operator < (const Line& L)const{
    return ang<L.ang;
  }
};

1.点积：两个向量的长度的乘积再乘上他们的夹角的余弦。

double Dot(Vector a,Vector b) { return a.x*b.x+a.y*b.y;}

2.求向量的长度：

double Length(Vector a)
{
  return sqrt(Dot(a,a));
}

3.叉积：等于两个向量所组成三角形的有向面积两倍。

向量A在B的下方时Cross>0，反之小于0，共线等于0.

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

4.向量的旋转：

Vector Rotate(Vector a,double rad)//rad是弧度
{
  return Vector(a.x*cos(rad)-a.y*sin(rad),a.x*sin(rad)+a.y*cos(rad));
}

5.向量的法向量：

Vector Normal(Vector a)
{
  double L=Length(A);
  return Vector(-a.y/L,a.x/L);
}

6.两条直线的交点：

前提是直线的单位向量得知道，然后在直线上取一个点，组成参数方程，下面的函数调用前必须保证有交点。

Point GLI(Point P,Vector v,Point Q,Vector w)
{
  Vector u=P-Q;
  double t=Cross(w,u)/Cross(v,w);
  return P+v*t;
}

7.点到直线的距离：

double DTL(Point p,Point a,Point b)
{
  Vector v=b-a,u=p-a;
  return fabs(Cross(v,u))/Length(v);//不取绝对值的话是有向距离
}

8.点到线段的距离：

注意区分点的投影在不在线段上，这个可看代码。

double DTS(Point p,Point a,Point b)
{
  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);
}

9.求投影点：

这里提及下点积有分配率即

Dot(v,p-(a+t0*v))

可以得到

Dot(v,p-a)-t0*Dot(v,v)

Point GLP(Point p,Point a,Point b)
{
  Vector v=b-a;
  return a+v*(Dot(v,p-a)/Dot(v,v));
}

10.线段相交的判定：

bool SPI(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 OS(Point p,Point a,Point b)
{
  return dcmp(Cross(a-p,b-p))==0&&dcmp(Dot(a-p,b-p))<0;
}

·包含端点的代码（注意要将一条线段的两个端点都调用判断下）：

bool OS(Point p,Point a,Point b)
{
  if(min(a.x,b.x)<=p.x<=max(a.x,b.x)&&(min(a.y,b.y)<=p.y<=max(a.y,b.y)))
  {
    return true;
  }
  else return false;
}

细心的人会发现我上面的基本都是考虑精度的代码，都用了dcmp函数，在不要求的情况下可以不用。

11.多边形的面积计算：

适用于凸多边形也适用于凹多边形，就是将多边形切成n-2个三角形。

double CPA(Point *p,int n)
{
  double area=0;
  for(int i=0;i<n;i++)
  {
    area+=Cross(p[i]-p[0],p[i+1]-p[0]);
  }
  return area/2;
}

12.点在多边形内判定

假想有一条向右的射线，统计多边形穿过这条射线正反多少次，把这个数记为绕数wn，逆时针穿过时，wn加1，顺时针穿过时，wn减1。

int IPIP(Point p,Polygon poly)
{
  int wn=0;
  int n=v.size();
  for(int i=0;i<n;i++)
  {
    if(OS(p,poly[i],poly[(i+1)%n])) return -1;//在边界上
    int k=dcmp(Cross(poly[(i+1)%n]-poly[i],p-poly[i]));
    int d1=dcmp(poly[i].y-p.y);
    int d2=dcmp(poly[(i+1)%n].y-p.y);
    if(k>0&&d1<=0&&d2>0) wn++;
    if(k<0&&d2<=0&&d1>0) wn--;
  }
  if(wn!=0) return 1;//在内部
  return 0;//在外部
}

下面附上代码部分总模板：

#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
#define INF 0x3f3f3f3f
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)
{
  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;
}
/*bool operator == (const Point& a, const Point& b){
    return a.x == b.x && a.y == b.y;
}*/
double Cross(Vector A,Vector B) { return A.x*B.y - A.y*B.x;}
Point GLI(Point P,Vector v,Point Q,Vector w)
{
  Vector u=P-Q;
  double t=Cross(w,u)/Cross(v,w);
  return P+v*t;
}
double Dot(Vector a,Vector b) { return a.x*b.x+a.y*b.y;}
double Length(Vector a)
{
  return sqrt(Dot(a,a));
}
Vector Rotate(Vector a,double rad)
{
  return Vector(a.x*cos(rad)-a.y*sin(rad),a.x*sin(rad)+a.y*cos(rad));
}
Vector Normal(Vector a)
{
  double L=Length(a);
  return Vector(-a.y/L,a.x/L);
}
double DTL(Point p,Point a,Point b)
{
  Vector v=b-a,u=p-a;
  return fabs(Cross(v,u))/Length(v);
}
double DTS(Point p,Point a,Point b)
{
  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 GLP(Point p,Point a,Point b)
{
  Vector v=b-a;
  return a+v*(Dot(v,p-a)/Dot(v,v));
}
bool SPI(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 OS(Point p,Point a,Point b)
{
  return dcmp(Cross(a-p,b-p))==0&&dcmp(Dot(a-p,b-p))<0;
}*/
bool OS(Point p,Point a,Point b)
{
  if(min(a.x,b.x)<=p.x<=max(a.x,b.x)&&(min(a.y,b.y)<=p.y<=max(a.y,b.y)))
  {
    return true;
  }
  else return false;
}
double CPA(Point *p,int n)
{
  double area=0;
  for(int i=0;i<n;i++)
  {
    area+=Cross(p[i]-p[0],p[i+1]-p[0]);
  }
  return area/2;
}

计算几何模板