计算几何算法1. 三角形和多边形的面积(2D及3D)

计算几何 everyinch 61598℃ 0评论

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

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

三角形

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

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

平行四边形: clip_image002[48] clip_image003[8]
 三角形: clip_image004[8] clip_image005[8]

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

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

clip_image008[8] clip_image009[8]

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

clip_image010[8] clip_image001[47] clip_image012[4]

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

clip_image013[8]

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

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

clip_image014[8] clip_image015[8]
现代三角形面积的计算方法

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

clip_image016[8]其中θ是v和w向量的夹角。对于3d三角形有顶点clip_image017[8],让clip_image018[8]clip_image019[8],可以得到:

clip_image020[8] clip_image001[48] clip_image021[8]

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

clip_image023[12]

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

clip_image024[8] clip_image001[49] clip_image025[8]

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

请注意,如果三角型 V0V1V2是逆时针方向面积将是负数,如果三角型顺时针方向是正数。所以它可以用来测试一个三角形的方向。测试点V2是在线段V0V1的左侧(正数)还是右侧(负数),所以它是非常有用的。

四边形

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

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

clip_image027[10]

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

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

clip_image029[8] clip_image001[50] clip_image030[8]

在2d空间,有顶点clip_image022[11]其中i=0~3,有:

clip_image031[8]

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

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

clip_image034[8]

从初等几何我们知道,三角形V0V1V2 的中线M0M1 平行于底V0V2。在三角形V0V2V3 的中线M2M3 也平行于V0V2。这样M0M1M2M3互相平行。相似地,M0M3M1M2 也互相平行,所以M0M1M2M3是一个平行四边形。这样我们能计算它的面积:

clip_image035[8]

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

clip_image037[10]

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

clip_image038[8] clip_image001[51] clip_image039[8]

多边形

2d多边形

一个二维多边形可以分解成三角形区域,对于计算面积而言,有比较容易的方法来分解简单多边形(没有交叉点)。多边形Ω有顶点clip_image036[11],其中i=0,n且Vn = V0。 让P为任意一点,对于Ω中的每条边clip_image040[12]构成三角形clip_image041[16]。那么,Ω的面积等于所有三角形clip_image042[20](其中i=0,n-1)有符号面积的和,我们有:

clip_image043[18] clip_image001[52] clip_image044[14]

注意到对于逆时针的多边形,即点P在边clip_image045[12]的内左侧,那么clip_image042[21]的面积是正值;相反,当点P在边clip_image045[13]的外右侧,那么clip_image042[22]的面积是负值。

例如,在上面的图中,三角形clip_image046[14]clip_image047[12]有正值的面积,同时也会发现在clip_image046[15]clip_image047[13]中有超出面积的部分。另一方面在三角形clip_image048[10]clip_image049[8]有负值的面积,它会抵消掉超出的部分。最终多余的部分都会被抵消掉,剩下的部分就是多边形Ω的实际面积。

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

clip_image051[8] clip_image001[53] clip_image052[8]

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

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

3d平面的多边形

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

标准公式

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

clip_image055[26]在一个3d空间的平面多边形clip_image056[26]有顶点clip_image057[26]其中clip_image058[8]clip_image059[10]。多边形中的所有顶点都位于3d空间的平面P 中,它有法向量n。让P成为3d空间中的任意一点(通常不在平面P 中)。对于clip_image056[27]中的每条边clip_image040[13]构成三角形clip_image041[17]。现在有一个金字塔形的锥体点P是这个锥体的顶点面clip_image056[28]是锥体的底。我们打算投射这个锥体的三角形到平面P 中,然后计算这些投射三角形有符号的面积。那么这些投射面积的总和就等于多边形的面积。

clip_image060[10]为了达到这个目标,首先连接每个三角形clip_image041[18],它的面向量为clip_image061[8]垂直于clip_image042[23],它的模等于三角形的面积。接下来,从点P向平面P 画垂线至P0clip_image062[8]是投射三角形。从P0到Bi画垂线P0Bi与边clip_image040[14]相交与Bi。既然PP0垂直于ei,三个点PP0Bi定义的平面也垂直于ei,也就是PBi垂直与ei。这样|PBi|是clip_image042[24]的高度,|PBi|是Ti的高度。这两个高度的夹角θ等于n和a之间的夹角。因为在PP0Bi面旋转了90度,这样:

clip_image063[8]

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

clip_image064[8]

最后,通过选择P值(0,0,0),则有PVi=Vi,这样可以生成更简洁的公式:

clip_image065[8] clip_image001[54] clip_image066[8]

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

四边形分解

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

那么多边形clip_image056[30](n>4个顶点)能分解成V0和其它三个序列顶点V2i-1,V2i和V2i+1(其中i=1,h,这里h等于≤(n-1)/2的最大整数),当n为奇数时,分解到最后将是一个三角形。即:

clip_image069[8]

其中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子平面。也就是投射为:Projx(x,y,z)=(y,z),Projy(x,y,z)=(z,x)和Projz(x,y,z)=(x,y)。这个被保护的多边形的符号和投射子平面的方向匹配。

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

clip_image070[8]

这个公式的结果为:

clip_image071[8]

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

作为结果,这个算法需要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)

喜欢 (11)
发表我的评论
取消评论

表情

Hi,您需要填写昵称和邮箱!

  • 昵称 (必填)
  • 邮箱 (必填)
  • 网址
(20)个小伙伴在吐槽
  1. 略微高深了,有时间一定看看!
    别闹大湿2014-07-14 16:47 回复
    • 慢慢学习,多多交流,谢谢支持!
      everyinch2014-07-17 13:15 回复
    • Thanks for spending time on the computer (wgirint) so others don't have to.
      Janae2016-06-21 01:49 回复
    • Many people are working hard to get traffic from search engines like Google, Yahoo and MSN. Traffic from search engines are of top quality, which means that the traffic is highly targeted and is the most likely to convert and earn commissions for affiliates. In this affiliate marketing tips article, we are going to look at 3 major reasons why an affiliate should not depend on search engine traffic alone and should diversify into leveraging social media when building traffic for their affiliate sites.BerniceH recently posted..
    • I totally agree! I was beginning to feel that possibly too much might be stressed about Pastor Luter being African-American and not enough about his stellar godly character and his outstanding qulifications to lead Southern Baptists. Let us as a aconvention, be wary of patting ourselves on the back. After saying this about patting ourselves on the back, I must say that I am proud to be a Southern Baptist.Bill WilliamsPeoria, IL
    • ohlala, elles tuent tes photos…elles sont vraiment magnifiques…et bien joué les retouches abdos sur cédric, lolla 2ème est vraiment géniale, la grossesse (même presque à terme) te va super bien
      http://www./2016-07-21 01:15 回复
  2. 留名以后看
    田策2014-07-15 15:24 回复
    • 嗯嗯,多多交流,谢谢支持
      everyinch2014-07-17 13:14 回复
      • This piece was a likceajfet that saved me from drowning.
        Kalyn2016-06-21 02:02 回复
      • easily passed on to others when it…is in writing.the ability for your message to be permanent, compelling, and shareable is very powerful. don’t dilute your message by slacking on the basic rules of spelling, grammar, and punctuation. rather, master these fundamentals to make your mess…
        http://www.kibo.us/2016-07-10 00:11 回复
    • Filyanl! This is just what I was looking for.
      Monkey2016-06-21 01:55 回复
    • I like the fact that your blog is so honest, and it’s one of the reasons I come back. You’re a talented writer and storyteller. You inspire me to get off my ass when I see you posting four times a day. That says a lot.
    • It's really great that people are sharing this information.
      http://www./2016-08-12 07:01 回复
  3. This info is the cat's paasmja!
    Jenn2016-06-21 02:10 回复
  4. Stevens-Henager College’s Associate of Occupational Studies in Medical Specialties program at its Logan, Provo/Orem, and Ogden/West Haven campuses is accredited by the CAHEP. Was this answer helpful?
  5. 博主可能未加以考证即简单转载引用。古希腊没有乘法,所以关于计算面积的种种纯属西方人杜撰。另外,中世纪以前没有代数符号,所以所谓的海伦公式也是欧洲人后来杜撰的
    江河2020-09-26 14:25 回复
'; } if( dopt('d_footcode_b') ) echo dopt('d_footcode'); ?>