最开头的话
与点,线,多边形,圆形等各种几何图形相关的算法称为计算几何算法。
计算几何已经称为3D图形,CAD,机器人等多中领域的计基础,在程序设计竞赛中也常常出现有关的问题。
计算几何包含很多的内容,涉及范围很广。但是,在竞赛中出现的题目大多集中于基础性的问题上。而解决这些问题的关键是,把本科线性代数和高中几何学的知识转换成代码。因此,这篇博客的主要内容是如何将基础数学理论转化为简洁而无异常的代码形式。
计算几何的工具
向量的实现
向量最直观的表现形式为带有方向的箭头。向量有长度和方向两个因素构成,所以箭头的起点并不重要。因此,可以将向量的起点定义为坐标空间的原点,那么就可以将向量表示成箭头的终位置(x,y)。经过这样的处理后,就能把向量变成带有两个成员变量的类。
以下的代码实现了这种方法的vector2类函数。并且实现了多种运算符的重载。
const double PI = 2.0 * acos(0.0);
struct vector2
{
double x, y;
//explicit 避免加入实数
explicit vector2(double x_, double y_):x(x_),y(y_){}
// 比较两个向量
bool operator == (const vector2& rhs) const
{
return (x == rhs.x && y == rhs.y);
}
bool operator < (const vector2& rhs) const
{
return x != rhs.x ? x < rhs.x : y < rhs.y;
}
//向量之间的加减
vector2 operator + (const vector2& rhs) const
{
return vector2(x+rhs.x, y+rhs.y);
}
vector2 operator - (const vector2& rhs) const
{
return vector2(x-rhs.x, y-rhs.y);
}
//向量乘实数
vector2 operator * (double rhs) const
{
return vector2(x*rhs, y*rhs);
}
//返回向量的长度
double norm(){return hypot(x, y);}
//返回单位向量
vector2 normalize()
{
return vector2(x/norm(), y/norm());
}
//返回从x轴正方向逆时针到达当前向量的角度
double polar() const
{
return fmod(atan2(y,x) + 2*PI, 2*PI);
}
//点乘
double dot(vector2& rhs) const
{
return x*rhs.x + y*rhs.y;
}
//叉乘
double cross(const vector2& rhs) const
{
return x * rhs.y - y * rhs.x;
}
};
其中可能存在大家不太了解的函数,麻烦自行百度。
点与直线、线段的表示方法
我们将把线段表示成以这两个端点为终点的两个向量。
将直线表示成包含于该直线的任意一条线段。
这表示的最直接的好处,是可以和已有的向量建立起联系,可以使用现成的函数。
向量的內积和叉积
这两者在计算几何中很常用,都有着各自的几何含义。
内积
定义
各种符号的,不好打,自行百度。
用处
- 向量的夹角
- 判断两个向量是否 相等
- 向量的投影
叉积
定义
略。。。
用处
- 计算面积
向量a和b的叉积绝对值等于将向量a和b用作两边的平行四边形的面积。 - 判断两个向量的方向
两个向量的叉积,如果为正数,就能判断a在b的逆时针方向的180度以内;
如果为负数,则在顺时针方向。
叉积的实现代码:
//原点向量在向量a的逆时针上,返回正数;顺时针 返回负数;平行 返回 0;
double ccw(vector2 a, vector2 b)
{
return a.cross(b);
}
//把点p视为基准点时,返回值情况同上;
double ccw(vector2 p, vector2 a, vector2 b)
{
return ccw(a-p, b-p);
}
相交、距离、面积
直线与直线相交
确认直线之间是否相交在几何题中很常见,但是,编写这样的代码却不太容易,以为要考虑各种特殊的情况。
表示相交直线的最好方式为,将直线表示成一个点和一个方向向量。
<因为不知道向量符号咋打的, 所以 下文中两个相同字符的视为向量符号加上单个字符>
例如:求解直线l1: aa + p * bb ; 直线 l2: cc + q * dd;
可得: aax + p * bbx = ccx + q * ddx; aay + p * bby = ccy + q * ddy;
所以:
化简得:
最后将p的值 带入直线l1 中即可。实现代码:
//计算两个直线(a,b)和直线(c,d)的交点
double lineIntersection(vector2 a, vector2 b, vector2 c, vector2 d, vector2& x)
{
double det = (b-a).cross(d-c);
if(fabs(det) < 1e-16) return false;
x = a + (b-a)*((c-a).cross(d-c)/det);
return true;
}
线段和线段相交,不需要求交点
自己画一下图,很好理解的。
//线段(a,b)和线段(c,d) 是否有交点。
bool segmentIntersection(vector2 a, vector2 b, vector2 c, vector2 d)
{
double ab = ccw(a,b,c)*ccw(a,b,d);
double cd = ccw(c,d,a)*ccw(c,b,d);
if(ab == 0 && cd == 0)
{
if(b < a) swap(a,b);
if(d < c) swap(d,c);
return !(b < c || d < a);
}
return ab <= 0 && cd <= 0;
}
凸包
先上张图片
这里写图片描述
先来了解一下,什么是凸包。
点集Q的凸包(convex hull)是指一个最小凸多边形,满足Q中的点或者在多边形边上或者在其内。
有人把这些点看成在木板上的钉子,而凸包就是一个橡皮筋套住所有钉子是的轮廓。
求解凸包的简单方法是卷包裹法。
卷包裹算法从一个必然在凸包上的点开始向着一个方向依次选择最外侧的点
很容易理解对不对。上代码:
//凸包(卷包裹法)
typedef vector<vector2> polygon;
polygon giftWrap(const vector<vector2>& points)
{
int n = points.size();
polygon hull;
//找出最左下角的点
vector2 pivot = *min_element(points.begin(), points.end());
hull.push_back(pivot);
while(true)
{
//找出最左边的点next
vector2 ph = hull.back(),next = points[0];
for(int i = 1; i < n; i++)
{
double cross = ccw(ph, next, points[0]);
double dist = (next - ph).norm() - (points[0] - ph).norm();
if(cross > 0 || (cross == 0 && dist < 0))
next = points[i];
}
if(next == pivot) break;
hull.push_back(next);
}
}