计算几何基础

一、精度控制

1.有时候我们直接按大于，小于，等于比大小会错，因为精度问题，计算几何经常牵扯到浮点数的运算，所以就会产生精度误差，因此我们需要设置一个eps（偏差值），一般取1e-7到1e-10之间，并用下面的函数控制精度。

int sgn(double d){
    if( fabs(d)<eps ) return 0;
    else if( d>0 ) return 1;
    else return -1;
}

2.尽量少用三角函数、除法、开方、求幂、取对数运算

3.尽量把公式整理成使用的以上操作最少的形式

(e.g.) (1.0/2.0*3.0/4.0*5.0/6.0)最好转化为(1.0*3.0*5.0/left ( 2.0*4.0*6.0 ight ))

4.在不溢出的情况下将除式比较转化为乘式比较

(e.g.) (frac{a}{b}> cLeftrightarrow a> b*c)

5.不要直接用等号判断浮点数是否相等！

(e.g.) 同样的一个数(1.5)如果从不同的途径计算得来，他们的存储方式可能会是(a=1.4999999)与(b=1.5000001)，此时使用(a==b)判断将返回(false)

解决方法1：误差判别法：

#define eps 1e-6
int dcmp(double x,double y){
    if( fabs(x-y)<eps ) return 0;
    if( x>y ) return 1;
    return -1;
}

解决方法2：化浮为整法：

在不溢出整数范围的情况下，可以通过乘上(10^{x})转化为整数运算，最后再将结果转化为浮点数输出

6.负零

输出时注意不要输出(-0)，(e.g)

int main()
{
    double a=-0.000001;
    printf("%.4f
",a);
    return 0;
}

7.使用反三角函数时，要注意定义域的范围，比如，经过计算 (x=1.000001)

double x = 1.000001;
double acx = acos(x);
//可能会返回runtime error

//此时我们可以加一句判断
double x = 1.000001;
if(fabs(x - 1.0) < eps || fabs(x + 1.0) < eps)
    x = round(x);
double acx = acos(x);

8.取整函数

int main(){
    double x=1.56666;
    int fx=floor(x);    //floor()向下取整函数
    int cx=ceil(x);     //ceil()向上取整函数
    int rx=round(x);    //四舍五入函数
    cout<<x<<' '<<fx<<' '<<cx<<' '<<rx<<endl;
    return 0;
}

 二：点与向量

1.定义

点的定义：二维平面下的点的表示只需要两个实数，即横坐标与纵坐标即可

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

向量定义：既有大小又有方向的量叫做向量，在计算机中我们常用坐标表示

这样看来，向量这个结构体貌似与点没有任何区别，因此我们可以

typedef Point Vector;

2.运算

1.加法运算

点与点之间的加法运算没有意义

点与向量相加得到另一个点

向量与向量相加得到另外一个向量

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

2.减法运算

两个点之间作差将得到一个向量，(A-B)将得到由(B)指向(A)的向量(overrightarrow{BA})

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

 3.乘法运算

向量与实数相乘得到等比例缩放的向量

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

4.除法运算

向量与实数相除得到等比例缩放的向量

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

5.小于运算（Left Then Low运算）

有时我们需要将点集按照(x)坐标升序排序，若(x)相同，则按照(y)坐标升序排序

bool operator < (const Point& a, const Point& b){
    if( a.x==b.x ) return a.y<b.y;
    return a.x<b.x;
}

此比较器将在Andrew算法中用到

而Graham Scan算法用到的比较器基于极角排序

 6.内积运算：又称数量积，点积

(alpha cdot eta = left | alpha ight |left | eta ight |cos heta)

对加法满足分配律

几何意义：(向量alpha 在向量eta 的投影{alpha }'（带有方向性）与eta 的长度乘积 )

(若alpha 与eta 的夹角为锐角，则其内积为正)

(若alpha 与eta 的夹角为钝角，则其内积为正)

(若alpha 与eta 的夹角为直角，则其内积为0)

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

7.外积运算：又称向量积，叉积

(alpha imes eta =left | alpha ight |left | eta ight |sin heta )

( heta) 表示向量(alpha) 旋转到向量(eta) 所经过的夹角

对加法满足分配律

几何意义：(向量alpha 与eta 所张成的平行四边形的有向面积)

判断外积符号：右手定则：

(alpha imes eta ：若eta 在alpha 的逆时针方向，则为正值;顺时针则为负值;两向量共线则为0)

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

8.常用函数：

取模（长度）

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

计算两向量夹角：返回值为弧度制下的夹角

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

计算两向量构成的平行四边形有向面积

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

计算向量逆时针旋转后的向量

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

计算向量逆时针旋转九十度后的单位法向量

Vector Normal(Vector A){//向量A左转90°的单位法向量
    double L=Length(A);
    return Vector(-A.y/L,A.x/L);
}

ToLeftTest

(判断折线overrightarrow{bc}是不是向overrightarrow{ab}的逆时针方向（左边）转向)

凸包构造时将会频繁用到此公式

bool ToLeftTest(Point a,Point b,Point c){
    return Cross(b-a,c-b)>0;
}

 总模板：

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 - (Point A, Point 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){
    if(a.x == b.x)
        return a.y < b.y;
    return a.x < b.x;
}
const double eps = 1e-6;
int sgn(double x){
    if(fabs(x) < eps)
        return 0;
    if(x < 0)
        return -1;
    return 1;
}
bool operator == (const Point& a, const Point& b){
    if(sgn(a.x-b.x) == 0 && sgn(a.y-b.y) == 0)
        return true;
    return false;
}
double Dot(Vector A, Vector B){
    return A.x*B.x + A.y*B.y;
}
double Length(Vector A){
    return sqrt(Dot(A, A));
}
double Angle(Vector A, Vector B){
    return acos(Dot(A, B)/Length(A)/Length(B));
}
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为弧度 且为逆时针旋转的角
    return Vector(A.x*cos(rad)-A.y*sin(rad), A.x*sin(rad)+A.y*cos(rad));
}
Vector Normal(Vector A){//向量A左转90°的单位法向量
    double L = Length(A);
    return Vector(-A.y/L, A.x/L);
}
bool ToLeftTest(Point a, Point b, Point c){
    return Cross(b - a, c - b) > 0;
}

 参考博客：here