9
$\begingroup$

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]

3d graphical output of DelaunayMesh

3D graphical output from the results of TetGenConvexHull

$\endgroup$
7
$\begingroup$

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]

enter image description here

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

enter image description here

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]

enter image description here

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

enter image description here

$\endgroup$
  • $\begingroup$ That is a perfect approach. Thanks for helping! $\endgroup$ – Dennis Nov 30 '20 at 2:36
  • $\begingroup$ But it seems that there is no documentation for function "RegionMeshMeshCellNormals". Could you please explain how it works? $\endgroup$ – Dennis Nov 30 '20 at 2:45
  • $\begingroup$ @Dennis, it is undocumented. I updated with some info that i found thru trial/error. $\endgroup$ – kglr Nov 30 '20 at 9:05
  • $\begingroup$ That's clear and very helpful! Thanks so much! $\endgroup$ – Dennis Nov 30 '20 at 14:12
  • $\begingroup$ 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? $\endgroup$ – Dennis Nov 30 '20 at 20:47
10
$\begingroup$

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

enter image description here

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Dennis Nov 29 '20 at 0:55
  • $\begingroup$ @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. $\endgroup$ – Michael E2 Nov 29 '20 at 1:46
  • $\begingroup$ 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 :-) $\endgroup$ – MarcoB Nov 29 '20 at 3:35
  • $\begingroup$ I obtained somewhat different results, making use of Maple. See here. The same is obtained with three digits. $\endgroup$ – user64494 Nov 29 '20 at 8:58
  • $\begingroup$ 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. $\endgroup$ – user64494 Nov 29 '20 at 11:35
7
$\begingroup$

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]

enter image description here

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

coPQ = coplanarQ[Join[#[[1]], #2[[1]]]]&
$\endgroup$
4
$\begingroup$

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

enter image description here

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
  ]
$\endgroup$
0
$\begingroup$

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

$\endgroup$
  • $\begingroup$ The PolyhedralSets package of Maple works in higher dimensions too (see here for more info). $\endgroup$ – user64494 Dec 1 '20 at 5:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.