首先定义编写计算几何代码时必须用到的概念,并观察其实现方法。 向量的实现 二维平面上的两个不同点的相对位置该如何表示呢?对于这个问题,应当先求出两点之间的距离,然后求出从一个点到另外一个点的方向。求出方向和距离后,就能准确表示出两点之间的相对位置。这种既有大小又有方向的距离就是向量(vector),向量是数学、物理学、工学等领域中广泛应用的概念,也是计算几何的最基础的工具。 向量最直观的表示方法是带有箭头的线段。向量由方向和长度两个因素组成,所以箭头的起点并不重要。即使改变了起点,向量也不会发生改变。因此,可以把向量的起点定义为坐标空间的原点,那么就可以把向量表示成箭头的终点位置(x, y)。经过这些处理后,就能把向量表示为带有两个成员变量的类函数。 代码15-1是实现这种方法的vector2类函数。代码中不仅包含了两个点的坐标,还实现了多种运算符的重载以及成员函数。vector2类的运算符重载函数实现了对实际向量的运算。向量的加法、减法,以及与实数的乘法可从几何学角度定义为如下形式。 两个向量a和b的相加向量a+b等于,将b的起点移到a的终点时,在b的终点结束的向量。把两个向量的坐标分别相加就能得到相加后的向量。 从一个向量a中减去另一个向量b后的向量ab等于,从b的终点到a终点的向量。把两个向量的坐标分别相减就能得到相减后的向量。这样定义时可知,(ab)+b=a。 一个向量a乘以实数r之后的向量a•r等于,将a的长度增加到r倍后的向量。向量的坐标分别乘以r就能得到相乘后的向量。 代码15-1 表示二维向量的vector2类函数 const double PI = 2.0 * acos(0.0); // 表示二维向量。 struct vector2 { double x, y; // 把构造函数指定为explicit,就可以防止发生 // 在vector2的位置中误加入实数的错误。 explicit vector2(double x_ = 0, double y_ = 0) : 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() const { return hypot(x, y); } // 返回相同方向的单位向量(unit vector)。 // 调用0向量时返回值无定义。 vector2 normalize() const { return vector2(x / norm(), y / norm()); } // 从x轴的正方向逆时针到达当前向量时的角度 double polar() const { return fmod(atan2(y, x) + 2*PI, 2*PI); } // 内积/叉积的实现 double dot(const vector2& rhs) const { return x * rhs.x + y * rhs.y; } double cross(const vector2& rhs) const { return x * rhs.y - rhs.x * y; } // 将当前向量映射到rhs的结果。 vector2 project(const vector2& rhs) const { vector2 r = rhs.normalize(); return r * r.dot(*this); } }; 不过,比较向量的大小较为特殊。该函数将终点的x坐标值越小的向量,以及x坐标值相同时y坐标值越小的向量判断成相对较小的向量。这种方法虽然没有数学意义,但很常用。本章的许多算法都将使用这种方法比较向量的大小。 代码中还需要关注计算极角(polar angle)的polar()函数。向量的极角指的是,该向量与x轴正方向之间逆时针算出的角度。为了计算此角度,函数中使用了atan2()函数。atan2()函数在给出(y, x)坐标时,根据各坐标的符号自动计算出此角度。该函数返回值的取值范围是(π, π],因此,polar()函数会将此返回值修正为[0, 2π)的一个值。fmod()函数会求出两个实数相除之后的余数。所以把返回的角度先加上2π,再求出与2π相除之后的余数,就能得到[0, 2π)的一个角度。 点与直线、线段的表示方法 利用向量就能很方便地表示出点、直线和线段。本书将二维平面上的点表示成以该点为终点的向量,将线段表示为以其两个端点为终点的两个向量,将直线表示为包含于该直线的任意一条线段。图15-1给出了这几种表示方式。 其实这种表示方法较为繁琐。点可以用两个坐标值的实数变量表示,线段可以用两个端点坐标的4个实数变量表示,而直线可以用直线方程ax+by+c=0的系数这3个实数变量表示。既然可以用这么简便的方式表示出点、线段和直线,那么为什么非要用向量表示呢?这是因为,利用向量表示这3个几何图形,能够让编写的代码更加简洁,而且向量可以成为解题的强大工具。例如,编写函数,使其在给出3个点p、a、b时,返回a到p距离和b到p距离之差。利用vector2()就可以编写成如下形式。 图15-1 利用向量表示的点、线段和直线 double howMuchCloser(vector2 p, vector2 a, vector2 b) { return (b - p).norm() - (a - p).norm(); } 反之,以两个实数表示各点时,需要编写如下形式的代码。 double howMuchCloser(double px, double py, double ax, double ay, double bx, double by) { return sqrt((bx - px) * (bx - px) + (by - py) * (by - py)) - sqrt((ax - px) * (ax - px) + (ay - py) * (ay - py)); } 当然,本书中的表示方法也不总是最好的方法。根据个人喜好或问题需要,也可以采用其他更好的表示方法。例如,可以把线段表示成一个端点和另一个坐标的相对坐标。将一个点表示成极坐标的话,有时对解决问题会更有利。重要的是,同一程序中选择的表示方式要始终如一,不应该混用不同的表示方法。 向量的内积和叉积 除加法和减法外,向量中还有几种常用的计算方式。其中的内积和叉积都具有几何意义,所以经常应用于计算几何算法。首先介绍向量的内积(inner product)。两个二维向量a=(ax, ay)和b=(bx, by)的内积a•b是利用如下两种方法算出的实数。 可以看到,第二个定义中突然出现θ,这个θ表示两个向量之间的角度。向量内积的意义就在于通过其求出两个向量之间的角度。虽然向量内积的意义在于等式的第二行,而实际计算内积时,用到的却是等式的第一行。代码15-1的dot()就是利用等式的第一行求出了向量的内积。 图15-2 向量的内积和叉积 向量的内积常用于以下用途。 求两个向量的夹角:内积的最基本的用途就是求两个向量的夹角。对等式的第二行整理后得出如下等式。 判断两个向量是否垂直:如果cosθ=0,那么θ只能是π/2或3π/2。因此,如果两个长度非0的向量a和b的内积等于0,那么二者只能是相互垂直的关系。虽然也可以利用前面介绍过的方法,但是此方法不必调用运算速度较慢的acos函数,而且具有数值稳定性。 向量投影:在向量b上垂直照射光线时,向量a在向量b上的投影就是a的向量投影(vector projection)。通过图15-2(a)可知,向量投影的长度是|a|cosθ。利用公式|a|cosθ=a•b/|b|就能很容易地计算出向量投影的结果。代码15-1中的project()函数实现了向量的投影,代码中先把rhs转换为单位向量,使代码更加简洁。 接下来介绍向量的叉积(cross product)。向量的叉积运算其实是对两个三维向量进行定义的计算,这种计算在给出两个三维向量a和b时,返回垂直于这两个向量的另一个向量。不过,所有二维向量都可以视为z坐标为0的三维向量,所以对二维向量也可以定义向量的叉积。 当然,我们并不是为了找出与两个二维向量垂直的向量而进行向量叉积的计算 ,叉积中最有用的是返回的向量长度和方向。可以利用如下两种方法计算二维向量a=(ax, ay)和b=(bx, by)的叉积a×b的长度。 为便于讲解,后文会把两个二维向量a和b的叉积长度直接称为a和b的叉积。与内积定义类似,此公式的第一行是计算叉积的方法,第二行是计算叉积的目的。代码15-1中的cross()函数会利用公式的第一行求出两个向量的叉积。 向量的叉积通常用于以下用途。 计算面积:如图15-2(b)所示,a和b的叉积绝对值等于将a和b用作两边的平行四边形的面积。此平行四边形的面积相当于把原点和两个向量a和b作为顶点的三角形面积的两倍,因此,求出叉积的绝对值后再除以2就能求出三角形的面积。此公式可扩展到三角形以外的其他多边形,所以可利用于求多边形面积的许多问题。15.6节将介绍利用叉积求简单多边形面积的方法。 判别两个向量的方向:计算平行四边形面积时需要对叉积取绝对值的原因在于,叉积的值有可能是负值。叉积的结果值会根据计算方向的不同而发生变化。请看下列计算式。 发生这种现象的原因在于,第二行定义中的θ并不只是两个向量之间的角度,而是表示从向量a到向量b按照逆时针方向需要旋转的角度。因为sinθ=sinθ,所以只要知道叉积的符号,就能判断出θ是正值还是负值。图15-3(a)中,θ是正值,所以叉积也为正;而(b)中,θ是负值,所以叉积也为负。 图15-3 利用叉积判断两个向量的方向 叉积的主要用途就是利用这种性质判断两个向量的方向。两个向量a和b的叉积a×b如果是正数,就能判断出向量b在向量a的逆时针方向的180度以内;如果叉积是负数,则可以判断出向量b在向量a的顺时针方向的180度以内。图15-3(c)表示给出a×b的符号时,判断b在a的什么方向。代码15-2是ccw()函数 ,给出两个向量时,该函数返回一个向量在另一个向量的顺时针方向还是逆时针方向。 代码15-2 判断两个向量方向的ccw()函数 // 原点向量b在向量a的逆时针方向上,则返回正数;若在顺时针方向上,则返回负数; // 如果相互平行,则返回0。 double ccw(vector2 a, vector2 b) { return a.cross(b); } // 把点p视为基准点时,向量b在向量a的逆时针方向上,则返回正数;若在顺时针方向上,则返回负数; // 如果相互平行,则返回0。 double ccw(vector2 p, vector2 a, vector2 b) { return ccw(a-p, b-p); } 利用ccw()就可以解决只用坐标值无法解决的问题。假设有3个表示点a、b、c的向量a、b、c,此时,按顺序连接a、b、c的两条线段会在b点向左拐还是向右拐呢?只用坐标值会很难解决此问题,不过,利用ccw(a, b, c)就很容易了。如果此函数的返回值是一个正数,则c在b的逆时针方向,所以线段会向左拐,否则会右拐或直行。 其实,可以利用计算向量极角的polar()编写ccw()函数。不过,利用叉积的方法不仅使代码更加简洁,而且运算速度也更快,数值上也更稳定。