Construct a polyhedron from the coordinates of its vertices and calculate the area of each face

Question

I have the 3D coordinates for a set of points. I want to construct the convex polyhedron with those points as vertices.

I know I can use functions like DelaunayMesh or TetGenConvexHull (ConvexHullMesh doesn't work for me because it sometimes generates wrong polyhedra) but when I plot the polyhedron there are always “extra edges” on each face, deriving from the triangulation. I wonder if I can get rid of them.

I also would like to get access to each face (without those extra edges) individually to calculate their area. How could I achieve this?

q={{-10., -7.5, -6.25}, {-10., -7.5, 6.25}, {-10., 0., -10.}, {-10., 0.,
   10.}, {-10., 7.5, -6.25}, {-10., 7.5, 
  6.25}, {-1.66667, -11.6667, -8.33333}, {-1.66667, -11.6667, 
  8.33333}, {2.5, -13.75, 0.}, {-2.5, 
  0., -13.75}, {1.66667, -8.33333, -11.6667}, {-2.5, 0., 
  13.75}, {1.66667, -8.33333, 11.6667}, {-1.66667, 
  11.6667, -8.33333}, {1.66667, 8.33333, -11.6667}, {-1.66667, 
  11.6667, 8.33333}, {1.66667, 8.33333, 11.6667}, {2.5, 13.75, 
  0.}, {10., -10., 0.}, {10., -6.25, -7.5}, {10., -6.25, 7.5}, {10., 
  6.25, -7.5}, {10., 6.25, 7.5}, {10., 10., 0.}};

DelaunayMesh[q]

Needs["TetGenLink`"]
{vertex, surface} = TetGenConvexHull[q[[1]]];
Graphics3D[{EdgeForm[{Thick, Black}], GraphicsComplex[vertex, Polygon[surface]]}, Boxed -> False]

Answer 1

We can use the function Region`Mesh`MeshCellNormals to get face normals and group triangles by face normals:

ClearAll[combineCoplanarFaces]
combineCoplanarFaces[t_: 10^-3][bmr_ ] :=  Module[{faces = MeshPrimitives[bmr, 2],
    normals = Round[Region`Mesh`MeshCellNormals[bmr, 2], t]}, 
 Values @ GroupBy[Transpose[{normals, faces}], First -> Last, 
    Polygon@#[[Last @ FindShortestTour @ #]] & @ MeshCoordinates[RegionUnion@##] &]]

Example:

bdm = BoundaryDiscretizeRegion@DelaunayMesh[q];


Graphics3D[{EdgeForm[Thick], Hue @ RandomReal[], #} & /@ combineCoplanarFaces[][bdm], 
 Boxed -> False, ImageSize -> Large]

Update: Region`Mesh`MeshCellNormals is undocumented.

What we know about it is limited to what Information @ Region`Mesh`MeshCellNormals returns; that is, it takes an option Method:

Information @ Region`Mesh`MeshCellNormals

Through trial/error, we find that it takes two arguments; the first argument is aMeshRegion or a BoundaryMeshRegion object and the second argument is (1) an integer indicating mesh cell dimension (0 for Points, 1 for Lines and 2 for faces) or (2) a pair of integers {dim, index} indicating mesh cell index or (3) a list of integer pairs {{dim1, index1}...}.

So for mesh cells with indices {0, 5}, {1, 4}and {2, 15} we get

Region`Mesh`MeshCellNormals[bdm, {{0, 5}, {1, 4}, {2, 15}}]

{{-3.71093, 2.60535, -1.63299},
 {0.486229, 0.464079, -0.740413}, 
 {0.40825, 0.816495, 0.408249}}

indices = {{0, 5}, {1, 4}, {2, 15}}; 

Show[bdm, 
 Graphics3D[{AbsoluteThickness[5], AbsolutePointSize[15], #} & /@ 
   MapThread[{Darker@#3, FaceForm[Lighter@#3], #, 
      Line[{If[Head[#] === Point, #[[1]], Mean@#[[1]]], 
        If[Head[#] === Point, #[[1]], Mean@#[[1]]] + 3 Normalize @ #2}]} &, 
    {MeshPrimitives[bdm, indices], 
     Region`Mesh`MeshCellNormals[bdm, indices], 
    {Red, Green, Blue}}]],
  ImageSize -> Large]

Row[Table[Show[bdm, 
   Graphics3D[{AbsoluteThickness[5], AbsolutePointSize[15], 
       Darker[rc = RandomColor[]], FaceForm[Lighter@rc], #} & /@ 
     MapThread[{#, Line[{If[i == 0, #[[1]], Mean@#[[1]]], 
          If[i == 0, #[[1]], Mean@#[[1]]] + 3 Normalize @ #2}]} &, 
       {MeshPrimitives[bdm, i], Region`Mesh`MeshCellNormals[bdm, i]}]], 
   ImageSize -> Medium], 
  {i, 0, 2}]]

Answer 2

[Update: combinepolys now works when face meshes have interior vertices.]

Start by getting the bounding polygons:

reg = DelaunayMesh[q]

bdypolys = Cases[Normal@Show[
    BoundaryMesh[reg]
    ], _Polygon, Infinity]

We group coplanar polygons into faces (Tolerance may need adjustment in other cases):

coplanarQ[pts_?MatrixQ] := 
  MatrixRank[Transpose@pts - pts[[1]], Tolerance -> 10^(-4)] == 2;

faces = RegionUnion @@@ 
   Gather[bdypolys, 
    coplanarQ[Flatten[{##} /. Polygon -> Identity, 1]] &];

We can get the areas:

Area /@ faces

(*
  {243.75, 174.693, 117.372, 117.372, 117.372, 117.372, 243.75, 
  117.372, 174.693, 174.693, 117.372, 174.693, 117.372, 117.372}
*)

The function combinepolys probably has horrible time complexity. On this small example, it's no problem. The trouble is that a common edge might be at the end points of the list of vertices of a polygon or consecutive entries. To combine two adjacent polygons, we need to match each of the four possible combinations.

combinepolys = # //. {
     {x___,
       {p___, a_Integer, b_Integer, q___}, y___,
       {s___, b_, a_, r___}, z___} :> {x, {p, a, r, s, b, q}, y, z},
     {x___,
       {p___, a_Integer, b_Integer, q___}, y___,
       {a_, r___, b_}, z___} :> {x, {p, a, r, b, q}, y, z},
     {x___,
       {b_Integer, p___, a_Integer}, y___,
       {s___, b_, a_, r___}, z___} :> {x, {b, p, a, r, s}, y, z},
     {x___,
       {b_Integer, p___, a_Integer}, y___,
       {a_, r___, b_}, z___} :> {x, {b, p, a, r}, y, z},
     (* update: cut out singular edges *)
     {x___, a_Integer, b_Integer, a_, y___} :> {x, a, y},
     {b_Integer, a_Integer, x___, a_} :> {a, x},
     {a_Integer, x___, a_, b_Integer} :> {a, x}
    } &;

Show /@ faces // combinepolys (* shows each individual face *)

Graphics3D[{EdgeForm[{Thick, Black}],
     {RandomColor[], Cases[Normal@#, _Polygon, Infinity]}
     }] & /@ combinepolys[Show /@ faces];
Show[%]

Answer 3

We use FindShortestTour to order the coordinates for each group of co-planar polygons:

ClearAll[coPQ, triangleCombine]

coPQ[x_, y_, t_: 10^-4] := Max[RegionDistance[InfinitePlane[x[[1]]], #] & /@ y[[1]]] <= t

triangleCombine[t_: 10^-4] := Module[{mc = MeshCoordinates[RegionUnion@##] & @@@ 
   Gather[MeshPrimitives[#, 2], coPQ[##, t] &]}, 
 Polygon @ #[[FindShortestTour[#][[2]]]] & /@ mc] &;

Example:

bdm = BoundaryDiscretizeRegion@DelaunayMesh@q;

gathered = Gather[MeshPrimitives[bdm, 2], coPQ];

Tally[Length /@ gathered]

 {{4, 2}, {3, 12}}

Graphics3D[{EdgeForm[AbsoluteThickness[3]], RandomColor[], #}& /@ triangleCombine[][bdm],
  ImageSize -> Large, Boxed -> False]

Note: We could have also defined coPQ using MichaelE2's coplanarQ as

coPQ = coplanarQ[Join[#[[1]], #2[[1]]]]&

Answer 4

It would be nice to be able to get the region boundary as a polygon or a (single) line in terms of the region. BoundaryMesh does not work on surface regions in 3D. In 2D, BoundaryMesh drops unused vertices, so you have to map the points in the boundary back to the points in original mesh to get the relation between them. This is needed if you project the face onto 2D coordinate system (parallel to the face plane, for instance.

The advantage in this approach is that it works when there are points in the interior of the faces in the original mesh. Here the setup is similar to my other answer except for the subdivided mesh:

reg = DiscretizeRegion[ConvexHullMesh[q]]; (* subdivides OP's example *)

bdypolys = Cases[Normal@Show[BoundaryMesh[reg]], _Polygon, Infinity];

Length@bdypolys
(*  754  <-- DelaunayMesh gives 44 *)

coplanarQ[pts_?MatrixQ] := (* same as my other answer *)    
  MatrixRank[Transpose@pts - pts[[1]], Tolerance -> 10^(-4)] == 2;
faces // Length
(*  14  <-- same as before *)

See below for the utilities boundingPolygon and findNormalProjection, which I think I might find useful in other work:

Graphics3D[{
  EdgeForm[{Thick, Black}],
  {RandomColor[], boundingPolygon[#]} & /@ faces
  }]

Utilities:

boundingPolygon[mr_MeshRegion] := 
  Module[{coords, proj, mesh2d, bmesh, idcs, invproj},
   coords = MeshCoordinates[mr];
   proj = findNormalProjection[coords];
   invproj = Nearest[coords.proj -> coords];
   mesh2d = MeshRegion[coords.proj, MeshCells[mr, 2]];
   bmesh = BoundaryMesh[mesh2d];
   GraphicsComplex[
    Flatten[invproj@MeshCoordinates[bmesh], 1],
    MeshCells[bmesh, 2]
    ]
   ];

ClearAll[findNormalProjection];
findNormalProjection // Options = {Tolerance -> 0.001}; (*
 * In an exact solution the matrix rank will be 2
 * In an approximate solution there may be two large
 *   and several small nonzero singular values 
 * If more than two singular values are large
 *   then the points are not coplanar
 *)
findNormalProjection[pts_, OptionsPattern[]] /; Length[pts] >= 3 :=
  Module[{u, w, v, centroid, mat, proj, orientation},
  centroid = Mean[pts];
  mat = Transpose@pts - centroid;
  {u, w, v} = SingularValueDecomposition[mat];
  proj = u[[All, 1 ;; 2]]; (* u[[All,3]] = normal to plane *)
  orientation = Cross[mat[[All, 1]].proj].(mat[[All, 2]].proj);
  If[orientation < 0,
   proj = proj.{{0, 1}, {1, 0}} (* reflect *)
   ];
  proj /;
   Count[Diagonal[w], s_ /; s > OptionValue@Tolerance] == 2
  ]

Answer 5

A long time ago, 25/04/1998, I wrote my own ConvexHull3D package. It returns the hull as:

{{1, 3, 10, 11, 7}, {1, 7, 9, 8, 2}, {2, 8, 13, 12, 4}, {3, 5, 14, 15,
   10}, {4, 12, 17, 16, 6}, {5, 6, 16, 18, 14}, {7, 11, 20, 19, 
  9}, {8, 9, 19, 21, 13}, {10, 15, 22, 20, 11}, {12, 13, 21, 23, 
  17}, {14, 18, 24, 22, 15}, {16, 17, 23, 24, 18}, {1, 2, 4, 6, 5, 
  3}, {19, 20, 22, 24, 23, 21}}

and, since your coordinates are integers divided by 12, you can even get exact results for volume, area(s) and Steven Finch's Mean Width.
Volume = 16125/2, Areas are

{575/(2 Sqrt[6]), (625 Sqrt[5])/8, 575/(2 Sqrt[6]), 575/(
 2 Sqrt[6]), 575/(2 Sqrt[6]), (625 Sqrt[5])/8, 575/(2 Sqrt[6]), 575/(
 2 Sqrt[6]), (625 Sqrt[5])/8, (625 Sqrt[5])/8, 575/(2 Sqrt[6]), 575/(
 2 Sqrt[6]), 975/4, 975/4}

and MeanWidth is an ugly expression numerically equal to 26.7631.

For those who don't abhor 'ugly programming', the link is still http://users.telenet.be/Wouter.Meeussen/ConvexHull3D.m