Construct a polyhedron from the coordinates of its vertices and calculate the area of each face
[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[%]
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 Point
s, 1
for Line
s 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}]]
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]]]]&