3D Polygon Areas

We continue the discussion initiated in the recent post on

two-dimensional polygon area calculation
.
We can use these areas to determine which of the edge loops is the outer loop versus inner loops defining holes.
The polygon with the largest area is the outer loop, all others are inner holes.
This is a continuation of the analysis for determining the profile boundary loop polygons for

floor slabs

and

walls
.
In this instalment, we expand our analysis to include three-dimensional polygons, i.e. planar polygons oriented arbitrarily in 3D space.

To determine the area of a two-dimensional polygon in the previous post, we used a simple formula which basically works by calculating the area of the triangles spanned by an arbitrary point in the plane and any two consecutive polygon vertices, and summing up all of those triangle areas.
The polygons returned for the floor and wall boundary loops are three dimensional.
We have two choices for calculating the area of a three-dimensional planar polygon in space.
One approach is to transform the polygon appropriately onto the XY plane and make use of the formula we already have.
Another approach makes use of the same algorithm adapted for use on the three-dimensional vertices directly.
It calculates each polygon’s plane data in space, i.e. its distance from the origin and normal vector.
The polygon area is by-product of these calculations, because it is equal to half of the length of the non-normalized normal vector.
I would like to demonstrate both approaches, so we can compare them and ensure that both return the same results.
We will begin with the implementation of the three-dimensional algorithm, including optimised implementations for the special cases of triangles and four-sided polygons:


bool GetPolygonPlane(
  List<XYZ> polygon,
  out XYZ normal,
  out double dist,
  out double area )
{
  normal = XYZ.Zero;
  dist = area = 0.0;
  int n = ( null == polygon ) ? 0 : polygon.Count;
  bool rc = ( 2 < n );
  if( 3 == n )
  {
    XYZ a = polygon[0];
    XYZ b = polygon[1];
    XYZ c = polygon[2];
    XYZ v = b - a;
    normal = v.Cross( c - a );
    dist = normal.Dot( a );
  }
  else if( 4 == n )
  {
    XYZ a = polygon[0];
    XYZ b = polygon[1];
    XYZ c = polygon[2];
    XYZ d = polygon[3];
 
    normal.X = ( c.Y - a.Y ) * ( d.Z - b.Z )
      + ( c.Z - a.Z ) * ( b.Y - d.Y );
    normal.Y = ( c.Z - a.Z ) * ( d.X - b.X )
      + ( c.X - a.X ) * ( b.Z - d.Z );
    normal.Z = ( c.X - a.X ) * ( d.Y - b.Y )
      + ( c.Y - a.Y ) * ( b.X - d.X );
 
    dist = 0.25 *
      ( normal.X * ( a.X + b.X + c.X + d.X )
      + normal.Y * ( a.Y + b.Y + c.Y + d.Y )
      + normal.Z * ( a.Z + b.Z + c.Z + d.Z ) );
  }
  else if( 4 < n )
  {
    XYZ a;
    XYZ b = polygon[n - 2];
    XYZ c = polygon[n - 1];
    XYZ s = XYZ.Zero;
 
    for( int i = 0; i < n; ++i ) {
      a = b;
      b = c;
      c = polygon[i];
 
      normal.X += b.Y * ( c.Z - a.Z );
      normal.Y += b.Z * ( c.X - a.X );
      normal.Z += b.X * ( c.Y - a.Y );
 
      s += c;
    }
    dist = s.Dot( normal ) / n;
  }
  if( rc )
  {
    double length = normal.Length;
    rc = !Util.IsZero( length );
    Debug.Assert( rc );
 
    if( rc )
    {
      normal /= length;
      dist /= length;
      area = 0.5 * length;
    }
  }
  return rc;
}

We define a variable n for the number of polygon vertices.
We have implemented specific code for triangles, i.e. the case n = 3, and 4-sided polygons.
The polygon area is half of the length of the non-normalized normal vector of the plane.

Here is the main section of the external command making use of this function.
It lists the areas of the various loops of the selected walls, or all walls in the model, if none were explicitly selected, and highlights the largest area.
Just like in the floor area calculation, the code assumes that only one wall is picked.
If more than one is processed, then the determination of the outer loop will only work for the wall with the largest surface area, since all loops of all walls are included in one single list.
We repackaged the functionality for determining the wall profile polygons into a separate utility method GetWallProfilePolygons() in the CmdWallProfile class, so that we can reuse that functionality from the new command as well:


List<List<XYZ>> polygons
  = CmdWallProfile.GetWallProfilePolygons(
    app, walls );
 
int i = 0, n = polygons.Count;
double[] areas = new double[n];
double d, a, maxArea = 0.0;
XYZ normal;
foreach( List<XYZ> polygon in polygons )
{
  GetPolygonPlane( polygon,
    out normal, out d, out a );
  if( Math.Abs( maxArea ) < Math.Abs( a ) )
  {
    maxArea = a;
  }
  areas[i++] = a;
}
 
Debug.WriteLine( string.Format(
  "{0} boundary loop{1} found.",
  n, Util.PluralSuffix( n ) ) );
 
for( i = 0; i < n; ++i )
{
  Debug.WriteLine( string.Format(
    "  Loop {0} area is {1} square feet{2}",
    i,
    Util.RealString( areas[i] ),
    ( areas[i].Equals( maxArea )
      ? ", outer loop of largest wall"
      : "" ) ) );
}

Here is a small example model with two walls selected and highlighted in red. Their boundary loops are represented by model lines in green, offset outward from the outer wall face by one foot:

Wall boundary loops

This is the list of the boundary loop areas of the two selected walls.
The first wall has an outer loop and two inner ones for the two windows.
The second wall has only one single simple rectangular outer loop:


4 boundary loops found.
Loop 0 area is 288.98 square feet, outer loop of largest wall
Loop 1 area is 2.67 square feet
Loop 2 area is 2.67 square feet
Loop 3 area is 172.22 square feet

There are still some aspects of this topic left that we would like to address in future posts.
As mentioned above, we would like to transform the 3D polygons into the 2D XY plane
so we can use the two-dimensional GetSignedPolygonArea() on those as well, instead of performing the area calculation in 3D.
Then we can compare the 2D and 3D results to ensure that they are equal.
It would also be interesting to compare the relative speed of transforming the 3D situation to 2D and using 2D area calculation versus direct 3D area calculation.

Here is an updated
version 1.0.0.15
of the complete Visual Studio solution,
including the two new commands CmdSlabBoundaryArea and CmdWallProfileArea as well as all other commands discussed so far.


Comments

5 responses to “3D Polygon Areas”

  1. Jarvis Aguero

    A big thank you for the blog article.Thanks Again.

  2. Hi Jeremy. I find your posts very interesting. I use them almost every day!
    Talking about the outer loop of a floor perimeter. Contrary to roofs, revit lets you create more than one loop with floors. There is a possibility that you find several floor perimeters and none of them is a hole. Therefore, you have to take all of them as outer loops.
    I think it is simpler to test if one of the vertex is inside of any of the other perimeters. Revit won´t let you create intersecting loops, so if a single point of the loop is inside the polygon, then all of them will be inside.
    To determine if a point is inside of a polygon you have to sum all angles between the point and the polygon vertexes. If the sum is 2PI (or zero), then it is inside. Graphically is easy to verify: you need complete the whole circle to join the point with all polygon’s vertexes.
    I attach here a simple routine.
    //Obtaining perimeters
    List<List> Poligonos
    = CmdSlabBoundary.GetFloorBoundaryPolygons(fl, opt);
    List<List> ExteriorPoligon = new List<List>();
    foreach (List poligon in Poligonos)
    {
    ExteriorPoligon.Add(poligon);
    }
    //I need to take out the interior loops
    foreach (List poligono in Poligonos)
    {
    //I need to check a point in each poligon
    XYZ ChkP = poligono[0];
    //With each other poligon.
    foreach (List pol2 in ExteriorPoligon)
    {
    if (pol2 == poligono) continue;
    double totalangle = 0;
    for (int i = 0; i < pol2.Count – 1; i++)
    {
    double angle =
    GetAngle(ChkP, pol2[i], pol2[i + 1]);
    totalangle += angle;
    }
    //And the last point is the first one again
    double angle2 =
    GetAngle(ChkP, pol2[pol2.Count – 1], pol2[0]);
    totalangle += angle2;
    winforms.MessageBox.Show(totalangle.ToString());
    //if the sum of all angles is 2
    Pi, then the point is interior to the poligon
    if (totalangle == 2 * Math.PI
    || totalangle ==0)
    {
    //If it is an interior poligon exclude it
    ExteriorPoligon.Remove(poligono);
    }
    }
    }
    private static double GetAngle(XYZ Or, XYZ P1, XYZ P2)
    {
    XYZ v0 = new XYZ(P1.X – Or.X, P1.Y – Or.Y, P1.Z – Or.Z);
    XYZ v1 = new XYZ(P2.X – Or.X, P2.Y – Or.Y, P2.Z – Or.Z);
    double scalarp = (v0.X * v1.X) + (v0.Y * v1.Y) + (v0.Z + v1.Z);
    double Mv0 = Math.Sqrt(Math.Pow(v0.X, 2) + Math.Pow(v0.Y, 2) + Math.Pow(v0.Z, 2));
    double Mv1 = Math.Sqrt(Math.Pow(v1.X, 2) + Math.Pow(v1.Y, 2) + Math.Pow(v1.Z, 2));
    double Cos = (scalarp) / (Mv0 * Mv1);
    double angle = Math.Acos(Cos);
    return angle;
    }

  3. Dear José,
    Thank you for the idea, code and appreciation.
    I published a similar polygon point containment algorithm here as well:
    http://thebuildingcoder.typepad.com/blog/2010/12/point-in-polygon-containment-algorithm.html
    It is actually not the same as what you suggest, because the version there is significantly optimised to use quadrants instead of the complete exact angle between points, significantly reducing the computational requirements and speeding up the algorithm.
    Cheers, Jeremy.

  4. Hi Jeremy,
    Thanks for the code that you have shown. Do you know if this supports concave polygons, or polygons that may have a hole? I am seeking a general solution for this calculation and would like to use yours. Also, are there any degenerate polygon plane rotations that would cause the this code to fail?
    Thanks,
    Jay

  5. Dear Jay,
    Thank you for your appreciation.
    This algorithm does support concave polygons.
    I order to support holes, I would suggest performing the calculation separately for the outer and inner loops and subtracting the latter from the former.
    Alternatively, use a full-fledged polygon library, e.g. the Generic Polygon Clipper GPC:
    http://thebuildingcoder.typepad.com/blog/2013/09/boolean-operations-for-2d-polygons.html
    That is definitely most suited for the general solution.
    Nope, I am not aware of any degenerate cases.
    I would suggest you add unit tests for all the possible types of cases that might appear in you usage scenarios.
    Cheers, Jeremy.

Leave a Reply

Discover more from Autodesk Developer Blog

Subscribe now to keep reading and get access to the full archive.

Continue reading