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

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]



We can use the function RegionMeshMeshCellNormals to get face normals and group triangles by face normals:

ClearAll[combineCoplanarFaces]
combineCoplanarFaces[t_: 10^-3][bmr_ ] :=  Module[{faces = MeshPrimitives[bmr, 2],
normals = Round[RegionMeshMeshCellNormals[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: RegionMeshMeshCellNormals is undocumented.

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

Information @ RegionMeshMeshCellNormals


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

RegionMeshMeshCellNormals[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], #} & /@
If[Head[#] === Point, #[[1]], Mean@#[[1]]] + 3 Normalize @ #2}]} &,
{MeshPrimitives[bdm, indices],
RegionMeshMeshCellNormals[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], RegionMeshMeshCellNormals[bdm, i]}]],
ImageSize -> Medium],
{i, 0, 2}]]


• That is a perfect approach. Thanks for helping! Commented Nov 30, 2020 at 2:36
• But it seems that there is no documentation for function "RegionMeshMeshCellNormals". Could you please explain how it works? Commented Nov 30, 2020 at 2:45
• @Dennis, it is undocumented. I updated with some info that i found thru trial/error.
– kglr
Commented Nov 30, 2020 at 9:05
• That's clear and very helpful! Thanks so much! Commented Nov 30, 2020 at 14:12
• I just implemented this method and also try to use function Area to find out the area of each face after grouping the faces with the same normal together. However, the Area (or RegionMeasure) function only works on some of the faces (polygons). A lot of the polygons won't be taken by Area to generate the result. I can't see why. For example, this issue pops up by using the vertex coordinates posted in the original question. Could you also clarify why this happens? Commented Nov 30, 2020 at 20:47

[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[%]


• That's exactly what I expected. Thank so much! I didn't know about "combinepolys" before and I will need to try understanding how this function works. Commented Nov 29, 2020 at 0:55
• @Dennis You're welcome. combinepolys works on indices found in MeshRegion/GraphicsComplex. In combinepolys[{{1, 2, 3}, {3, 2, 5}}] there are two triangles that share an edge between vertices 2 and 3. The result is a quadrilateral {{1, 2, 5, 3}}. I thought I wrote such a function years ago or discovered a built-in one, but I couldn't find it. Commented Nov 29, 2020 at 1:46
• Michael, Loved your solution! I’d been thinking about a way of identifying and merging co-planar faces since I read the question, but I came up empty. I’m relieved :-) Commented Nov 29, 2020 at 3:35
• I obtained somewhat different results, making use of Maple. See here. The same is obtained with three digits. Commented Nov 29, 2020 at 8:58
• A typo in the code (It should be [[-10., -7.5, -6.25]] instead of [-10., -7.5, -6.2].) causes the difference. Correcting the typo, the results by Maple coinside with the results of the answer under consideration. See Dropbox. Commented Nov 29, 2020 at 11:35

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]]]]&


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
]


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

• The PolyhedralSets package of Maple works in higher dimensions too (see here for more info). Commented Dec 1, 2020 at 5:51