• 计算几何算法1. 三角形和多边形的面积（2D及3D）
• 计算几何算法2. 关于线和点到线的距离（二维和三维）
• 计算几何算法3. 包含于多边形内的点的快速绕数
• 计算几何算法4. 关于平面以及点到平面的距离
• 计算几何算法5. 直线、线段和平面相交（2D和3D）

三角形
三角形面积的古老计算方法
三角形面积的现代计算方法
四边形
多边形
二维多边形
三维平面多边形
标准公式
四边形分解
投影到二维
实现
isLeft()
orientation2D_Triangle()
area2D_Triangle()
orientation2D_Polygon()
area2D_polygon()
area3D_Polygon()
参考资料

计算平面中多边形的面积是一个基本的几何计算问题，可以发现有许多的介绍性文字。但依据不同的给定条件可以有许多不同的计算方法。

三角形

古代三角形面积的计算方法

在毕达哥拉斯之前，平行四边形的面积（包括长方形和正方形）等于底乘以高。而且，两个同样的三角形可以合并在一起，形成一个平行四边形。因此，三角形的面积等于底乘以高的一半。因此，我们有：

平行四边形：
三角形：

然而，除了在特殊情况下，找到一个三角在任意方向上的的高度，通常需要计算从顶点到底的垂直距离。

例如，如果我们知道三角形两条边的长度a和b，和它们之间的夹角，那么就足以决定一个三角形的面积。用三角学的知识，相对于底b的高可以通过 求得，那么这样三角形的面积就是：

早在约公元前300年前的欧几里德发现三角形的面积是它的三条边的函数。这个发现在约公元50年达到顶峰，即发现了海伦公式：

其中a，b，c是边的长度， 和 s是半周长。下面是这个公式经过代数变换的版本：

通过避免计算3个平方根从而显著地提高了效率。

如果知道三个角中的两个角q 和φ和这个两个角的夹边b，那么这个三角形的面积是：

现代三角形面积的计算方法

从17世纪的笛卡尔和费马开始，线性代数的进展产生了新的面积公式。在3d空间，平面中平行四边形或三角形的面积可以通过计算两条边向量的叉乘的数量来计算，既然

其中θ是v和w向量的夹角。对于3d三角形有顶点，让和，可以得到：

2d向量(x,y)可以看做嵌在3d向量中且z组件为0的特例，可以通过计算2d向量的叉乘来计算面积。给定一个三角形有顶点其中i=0,1,2，能够计算：

z分量的绝对值是三角形面积绝对值的两倍，这里绝对值没有什么实际用途，反而让面积是一个有符号的数值是有意义的。

这是一个面积公式三角函数涉及非常有效的计算，只有2次乘法和5次加法，只有1个除法（这有时是可以避免的）。

请注意，如果三角型 V₀V₁V₂是逆时针方向面积将是负数，如果三角型顺时针方向是正数。所以它可以用来测试一个三角形的方向。测试点V₂是在线段V₀V₁的左侧（正数）还是右侧（负数），所以它是非常有用的。

四边形

古希腊人挑选特殊的四边形，其中包括正方形、长方形、平行四边形和梯形。那么给定一个任意四边形，他们给出了构造相同面积的平行四边形或正方形的方法。平行四边形的面积等于底乘以高，但这不是计算平行四边形的通用公式。

古印度数学家婆罗门笈多把海伦三角形面积公式扩展到四边形，但它仅仅工作在圆内接四边形的情形。对于圆内接多边形，让四个边的边长为a,b,c,d，半周长为s=(a+b+c+d)/2，那么Q 的面积为：

这是一个迷人的对称的公式。如果一边长为0，比如d = 0，那么可以退化为三角形的海伦公式。

在现代的线性代数，正如已经提到的那样，平行四边形的面积是两个相邻边向量的叉乘。因此，对于3d平面中的平行四边形，我们有：

在2d空间，有顶点其中i=0~3，有：

这又是一个有符号的面积，它也可以用来表示方向。

接下来，考虑任意四边形的情形。对于给定任意的四边形，有其四条边的中点构成一个新的平行四边形，称为伐利农平行四边形，它精确等于原始四边形的一半。对于任意四边形，每条边的中点显示如下：

从初等几何我们知道，三角形V₀V₁V₂ 的中线M₀M₁ 平行于底V₀V₂。在三角形V₀V₂V₃ 的中线M₂M₃ 也平行于V₀V₂。这样M₀M₁和M₂M₃互相平行。相似地，M₀M₃和M₁M₂ 也互相平行，所以M₀M₁M₂M₃是一个平行四边形。这样我们能计算它的面积：

它等于四边形对角线叉乘数量的一半。当限制在2d向量有，有下面任意四边形的有符号面积：

这个公式对于任意四边形而言是有效率的，只需要2个乘法和5个加法。对于简单四边形，但顶点是逆时针方向时它的面积值为正值，顺时针方向是面积为负值。但它对非简单四边形也是有效的，等于四边形边界的两个区域面积的差值。例如，在下图中顶点I是非简单四边形的自交叉点，我们有：

多边形

2d多边形

一个二维多边形可以分解成三角形区域，对于计算面积而言，有比较容易的方法来分解简单多边形（没有交叉点）。多边形Ω有顶点，其中i=0,n且V_n = V₀。 让P为任意一点，对于Ω中的每条边构成三角形。那么，Ω的面积等于所有三角形（其中i=0,n-1）有符号面积的和，我们有：

注意到对于逆时针的多边形，即点P在边的内左侧，那么的面积是正值；相反，当点P在边的外右侧，那么的面积是负值。

例如，在上面的图中，三角形和有正值的面积，同时也会发现在和中有超出面积的部分。另一方面在三角形和有负值的面积，它会抵消掉超出的部分。最终多余的部分都会被抵消掉，剩下的部分就是多边形Ω的实际面积。

通过将点P赋值为P=(0,0)使公式更加明确化，每个三角形的面积公式推导为。这样：

通过通过简单的代数转换就可以使公式变得更有效率。对于有n个顶点的多边形，第一个公式需要2n次乘法和(2n-1)次加法运算；第二个公式需要n次乘法和(3n-1)次加法；第三个公式仅仅需要n次乘法和(2n-1)次加法。所以第三个公式更有效率。

这个这个公式给出一个多边形的有符号面积，和三角形的有符号面积一样，当多边形的顶点是逆时针排列的它的面积为正值，当顶点是顺时针排列的面积是负值。所以，这个公式可以用来决定多边形的方向，但有更具效率的方法来判断多边形的方向。最容易的方法是查找多边形中最右下的顶点，测试进入和离开这个顶点的边的方向。如果离开边的最后一个顶点位于进入边的左侧，则意味着是逆时针的；否则是顺时针。

3d平面的多边形

我们已经给出3d空间三角形的面积是两个边向量的叉乘的模。也就是，

标准公式

对于3d多边形的面积有一个经典的标准公式，它源于斯托克斯定理，继承了三角形的叉乘公式。这里我们使用3d三角形分解的方法来推导这个公式，这个方法在集合意义上是更直观的。

在一个3d空间的平面多边形有顶点其中且。多边形中的所有顶点都位于3d空间的平面P 中，它有法向量n。让P成为3d空间中的任意一点（通常不在平面P 中）。对于中的每条边构成三角形。现在有一个金字塔形的锥体点P是这个锥体的顶点面是锥体的底。我们打算投射这个锥体的三角形到平面P 中，然后计算这些投射三角形有符号的面积。那么这些投射面积的总和就等于多边形的面积。

为了达到这个目标，首先连接每个三角形，它的面向量为垂直于，它的模等于三角形的面积。接下来，从点P向平面P 画垂线至P₀，是投射三角形。从P₀到B_i画垂线P₀B_i与边相交与B_i。既然PP₀垂直于e_i，三个点PP₀B_i定义的平面也垂直于e_i，也就是PB_i垂直与e_i。这样|PB_i|是的高度，|PB_i|是T_i的高度。这两个高度的夹角θ等于n和a之间的夹角。因为在PP₀B_i面旋转了90度，这样：

如果Ti是逆时针方向的，那么它的有符号面积就是正值。将所有的三角形Ti相加就得到了多边形的面积。则有：

最后，通过选择P值(0,0,0)，则有PV_i=V_i，这样可以生成更简洁的公式：

它使用了6n+3次乘法和4n+2次加法。

四边形分解

[凡盖德，1995]显示使用分解成四边形而不是三角形从而大大加快了计算过程，前面已经介绍过，3d空间的四边形能够通过计算对角线的叉乘的方法来计算它的面积。也就是，它减少了昂贵的叉乘计算。

那么多边形(n>4个顶点)能分解成V₀和其它三个序列顶点V_2i-1,V_2i和V_2i+1（其中i=1,h，这里h等于≤(n-1)/2的最大整数），当n为奇数时，分解到最后将是一个三角形。即：

其中n为奇数时k = 0，n为偶数时 k = n-1。这个公式减少了昂贵的叉乘计算。总共有3n+3次乘法和5n+1次加法运算，大约是经典公式的两倍快。

[Van Gelder, 1995]也指出，这种方法可以适用于二维多边形面积的计算。但他没有公布细节。利用这种方法生成的公式需要n次乘法和3n-1次加法，它没有之前给出的2d公式快，这里只是简单地提到并没有更进一步地说明。

投影到二维

通过投影多边形到2d平面可以更加加快3d空间平面多边形面积的计算。多边形的面积能够使用2d中最快的公式来计算，通过使用面积缩放因子来恢复3d的多边形面积。这个方法通过投射多边形到轴对齐的2d子平面来进行计算。为了得到正确的面积符号，2d自平面的基础向量必须排序。如果有(x,y,z)3d坐标，通过忽略x,y和z分量将它投射到yz,zx和xz子平面。也就是投射为：Proj_x(x,y,z)=(y,z)，Proj_y(x,y,z)=(z,x)和Proj_z(x,y,z)=(x,y)。这个被保护的多边形的符号和投射子平面的方向匹配。

接下来，选择投影最好避免退化和优化健壮性。多边形的法向量n，选择有最大绝对值的分量作为忽略的分量。让Projc()作为忽略的分量，其中c=x,y或z。那么投影多边形和原始多边形的面积比率为：

这个公式的结果为：

这个面积向量a是3d多边形的法向量，它的长度是3d多边形的面积。所以相对应地，如果3d多边形面积A()已知，n是面的单位法向量，那么A()n是面积向量，它的分量是投影的面积。这样，通过计算一个投影多边形的面积，就可以快速地得到另一个投影多边形的面积。

作为结果，这个算法需要n+5次乘法，2n+1次加法,1次平方根（当n不是单位法向量），加上选择忽略坐标的小开销，这个算法相对于传统3d多边形面积公式有惊人的效率提升，达到了6倍的性能。

实现

这里包括一些算法的C++实现代码。仅仅给出有整数坐标的2d情况，对于点、三角形和多边形只给出了最简单的结果。我们表示多边形只是坐标点的数组，但它常常使用顶点的链表结构是更为方便的（在绘制操作过程中允许插入和删除操作）

// This code may be freely used and modified for any purpose
// providing that this copyright notice is included with it.
// SoftSurfer makes no warranty for this code, and cannot be held
// liable for any real or imagined damage resulting from its use.
// Users of this code must verify correctness for their application.
// a Point (or vector) is defined by its coordinates
typedef struct {int x, y, z;} Point;
// exclude z for 2D
// a Triangle is given by three points: Point V0, V1, V2
// a Polygon is given by:
//        int n = number of vertex points
//        Point* V[] = an array of points with V[n]=V[0], V[n+1]=V[1]
// Note: for efficiency low-level functions are declared to be inline.
// isLeft(): tests if a point is Left|On|Right of an infinite lin
//    Input:  three points P0, P1, and P2
//    Return: >0 for P2 left of the line through P0 and P1
//            =0 for P2 on the line
//            0 for counterclockwise
 //            =0 for none (degenerate)
//            0 for counterclockwise
 //            =0 for none (degenerate)
//            //    Note: this algorithm is faster than computing the signed area.
int orientation2D_Polygon( int n, Point* V )
{
    // first find rightmost lowest vertex of the polygon
    int rmin = 0;
    int xmin = V[0].x;
    int ymin = V[0].y;

    for (int i=1; i<n; i++) { if (V[i].y > ymin)
            continue;
        if (V[i].y == ymin) {    // just as low
            if (V[i].x < xmin)   // and to left
                continue;
        }
        rmin = i;          // a new rightmost lowest vertex
        xmin = V[i].x;
        ymin = V[i].y;
    }

    // test orientation at this rmin vertex
    // ccw <=> the edge leaving is left of the entering edge
    if (rmin == 0)
        return isLeft( V[n-1], V[0], V[1] );
    else
        return isLeft( V[rmin-1], V[rmin], V[rmin+1] );
}
//===============================================================
// area2D_Polygon():computes the area of a 2D polygon
//    Input:  int n = the number of vertices in the polygon
//            Point* V = an array of n+2 vertices
 //                       with V[n]=V[0] and V[n+1]=V[1]
//    Return: the (float) area of the polygon
float area2D_Polygon( int n, Point* V )
{
     float area = 0;
     int   i, j, k;     // indices

     for (i=1, j=2, k=0; i<=n; i++, j++, k++) {
         area += V[i].x * (V[j].y - V[k].y);
     }
     return area / 2.0;
}

//===============================================================
// area3D_Polygon():computes the area of a 3D planar polygon
//    Input:  int n = the number of vertices in the polygon
//            Point* V = an array of n+2 vertices in a plane
//                       with V[n]=V[0] and V[n+1]=V[1]
//            Point N = unit normal vector of the polygon's plane
//    Return: the (float) area of the polygon</pre>

float area3D_Polygon( int n, Point* V, Point N )
{
     float area = 0;
     float an, ax, ay, az;  // abs value of normal and its coords
     int   coord;           // coord to ignore: 1=x, 2=y, 3=z
     int   i, j, k;         // loop indices

     // select largest abs coordinate to ignore for projection
     ax = (N.x>0 ? N.x : -N.x);     // abs x-coord
     ay = (N.y>0 ? N.y : -N.y);     // abs y-coord
     az = (N.z>0 ? N.z : -N.z);     // abs z-coord

     coord = 3;                     // ignore z-coord
     if (ax > ay) {
         if (ax > az) coord = 1;    // ignore x-coord
     }
     else if (ay > az) coord = 2;   // ignore y-coord

     // compute area of the 2D projection
     for (i=1, j=2, k=0; i<=n; i++, j++, k++)
         switch (coord) {
         case 1:
             area += (V[i].y * (V[j].z - V[k].z));
             continue;
         case 2:
             area += (V[i].x * (V[j].z - V[k].z));
             continue;
         case 3:
             area += (V[i].x * (V[j].y - V[k].y));
             continue;
         }

     // scale to get area before projection
     an = sqrt( ax*ax + ay*ay + az*az);  // length of normal vector
     switch (coord) {
     case 1:
         area *= (an / (2*ax));
         break;
     case 2:
         area *= (an / (2*ay));
         break;
     case 3:
         area *= (an / (2*az));
     }
     return area;
}
//==================================================

参考文献

Kevin Brown, MathPages : Heron’s Formula
Donald Coxeter & Samuel Greitzer,Geometry Revisited (1967)
Ronald Goldman, “Area of Planar Polygons and Volume of Polyhedra” in Graphics Gems II (1994)
Joseph O’Rourke, Computational Geometry in C (2nd Edition), Section 1.3 “Area of a Polygon” (1998)
J.P. Snyder and A.H. Barr, “Ray Tracing Complex Models Containing Surface Tessellations”, ACM Comp Graphics 21, (1987)
Dan Sunday, “Fast Polygon Area and Newell Normal Computation”, journal of graphics tools (jgt) Vol 7(2), (2002)
Allen Van Gelder, “Efficient Computation of Polygon Area and Polyhedron Volume” in Graphics Gems V (1995)
Eric Weisstein, MathWorld Site: Triangle or in The CRC Concise Encyclopedia of Mathematics (1998)

转载请注明：陈童的博客 » 计算几何算法1. 三角形和多边形的面积（2D及3D）