Face vertices of Polyhedron objects are not actually coplanar?

Question

I'm trying to draw "random" polyhedra with exactly $n$ faces. I figured the easiest way to do so was to get an "evenly-spaced" set of $n$ points on the sphere, get their ConvexHullMesh , then get the DualPolyhedron of the result. There's one problem:

The coordinates of the resulting DualPolyhedron make it so more often than not, its faces are not coplanar.

convexPoly = ConvexHullMesh[myPoints];
dualPoly = DualPolyhedron[convexPoly];
dualPolyMesh = BoundaryMeshRegion[
         dualPoly[[1]], Map[Polygon[#] &, dualPoly[[2]]] ]

This will give :coplnr: The vertices in the polygon Polygon[...] are not coplanar. error. Polyhedron objects are more lenient towards non-coplanar face vertices than BoundaryMeshRegion objects.

Is there any way to fix this so the resulting face-vertices of the DualPolyhedron will be coplanar?

thorimur · Answer 1 · 2021-07-31T23:44:08.930

Here's a way to get a random polyhedron with n faces by a slightly different method: take n random points on the sphere, then take the half-space bounded by the plane tangent to the sphere at each point, and intersect them all.

First, though, we need to be aware that we might encounter a pitfall when randomly choosing our points on the sphere: namely, we might end up with an unbounded region if all the points are too close together. The following function takes three points (assuming they will not be collinear, since the chance of that is infinitesimal), then chooses the fourth from a region that will "close off" the resulting bounded volume. I get that region with Coconic3 here. The remaining points after the first 4 are chosen randomly.

Coconic3[pts : (Repeated[{_, _, _}, 3])] := 
 RegionUnion @@ (Apply[HalfSpace[-Projection[#1, Cross[#2, #3]]] &] /@
     Partition[pts, 3, 1, {1, 1}])
RandomTetrahedronPoints[] := (pts |-> 
    Append[pts, 
     RandomPoint[
      DiscretizeRegion@RegionDifference[Sphere[], Coconic3[pts]]]])@
  RandomPoint[Sphere[], 3]
RandomSpherePoints[n_Integer] := 
 Join[RandomTetrahedronPoints[], RandomPoint[Sphere[], n - 4]] /; 
  n >= 4

I tried at first to use the built-in HalfSpace for the next part, but couldn't figure out a way to get a polygonal mesh out of it. (Any comments detailing how are welcome!) So, I opted for the homebrewed route. This looks at all intersection points of each subset of 3 planes, determines if the resulting point lies inside the volume bounded by the remaining planes, and if so, collects that point via Sow/Reap.

Since we know the resulting polyhedron will be convex by our construction, we then apply ConvexHullMesh.

PartitionPlanes[pts : {{_, _, _} ...}] := {#, Complement[pts, #]} & /@
   Subsets[pts, {3}]
InsideHull[solpt : {, _, _}, restpts : {{, _, _} ...}] := 
 AllTrue[restpts, solpt . # <= 1 &]
GetPoint[{subset : {Repeated[{_, _, _}, 3]}, comp : {{_, _, _} ...}}] := 
 If[InsideHull[#, comp], Sow[#]] &@
  Solve[And @@ ((pt |-> pt . {x, y, z} == 1) /@ subset), {x, y, 
     z}][[1, All, 2]]
GetVertices[planepts : {{_, _, _} ...}] := 
 First@Last@Reap[GetPoint /@ PartitionPlanes[planepts];]
RandomSolid[n_Integer] := ConvexHullMesh @ GetVertices @ RandomSpherePoints[n]
(* Test: *)
RandomSolid[6]
Length@MeshCells[RandomSolid[6], 2]
(* Out: 6 *)

One downside is that the constructed polyhedra can be very long. This could be rectified by changing our point-choosing procedure.

score 2 · Answer 2 · answered Aug 01 '21 at 01:18

Kinda looks like a bug in DualPolyhedron. You should report it, and see what WRI says.

Borrowing & adapting from my answer, https://mathematica.stackexchange.com/a/132920:

ClearAll[dual, sortvertices, reciprocate];
sortvertices[coords_, normal_, face_] := 
  With[{proj = DeleteCases[
       Orthogonalize[Join[{normal}, N@IdentityMatrix[3]]],
       {0., 0., 0.}][[2 ;; 3]],
    centroid = Mean[coords[[face]]]}, 
   SortBy[face, ArcTan @@ (proj . (coords[[#]] - centroid)) &]];
reciprocate[face_?MatrixQ, r_ : 1] /; Length[face] >= 3 := 
  r^2 {1, -1, 1} Most[#]/Last[#] &@
   Reverse@ Last@ Minors@
      Join[{{0, 0, 0, 0}},(dummy row)
       PadRight[face[[;; 3]], {Automatic, 4}, 1]];
dual[Graphics3D@GraphicsComplex[coords_, Polygon[faces_]], r_ : 1] := 
  With[{nvertices = Max@faces, nfaces = Length@faces}, 
   With[{mat = SparseArray@ Flatten@
        Table[{v, f} -> 1, {f, nfaces}, {v, faces[[f]]}], 
     dualcoords = reciprocate[coords[[#]], r] & /@ faces}, 
    With[{dualfaces = mat["AdjacencyLists"]}, 
     Graphics3D@ GraphicsComplex[dualcoords, 
       Polygon[Table[
         sortvertices[dualcoords, coords[[v]], dualfaces[[v]]],
         {v, Length@dualfaces}]]]]]];
dual[Graphics3D@GraphicsComplex[coords_, g_, ___], r_ : 1] :=
  With[{res = dual[
      Graphics3D@ GraphicsComplex[coords, 
        Replace[Cases[g, _Polygon, Infinity], 
         p : {__Polygon} :> Polygon[Join @@ p[[All, 1]]]]], r]},
   res /; Head@res === Graphics3D];

Example. Dual in red. A good radius r for the polar reciprocation depends on the polyhedron; nearly coplanar faces makes the corresponding vertices far from the center (e.g. random seeds 41, 45, 50, 66).

SeedRandom[0];
myPoints = Transpose[Transpose@# - Mean[#] &@RandomReal[{-10, 10}, {10, 3}]];
convexPoly = ConvexHullMesh[myPoints, WorkingPrecision -> Infinity];
convexPolyG = 
  Graphics3D[
   GraphicsComplex[MeshCoordinates@convexPoly, 
    First@convexPoly@"IndexedBoundaryPolygons"]];
Show[
 MapAt[{Opacity[0.5], #} &, convexPolyG, {1, 2}],
 (* dual *)
 MapAt[
  {Opacity[0.4], Red, #} &, 
  dual[convexPolyG,
   Max[Norm /@ myPoints]/Sqrt@3], (* often a nice radius *)
  {1, 2}]
 ]

Example. OP's sphere points.

myPoints = SpherePoints[11];
convexPoly = ConvexHullMesh[myPoints, WorkingPrecision -> Infinity];
convexPolyG = 
  Graphics3D[
   GraphicsComplex[MeshCoordinates@convexPoly, 
    First@convexPoly@"IndexedBoundaryPolygons"]];
Show[
 MapAt[{Opacity[0.5], #} &, convexPolyG, {1, 2}],
 (* dual *)
 MapAt[
  {Opacity[0.4], Red, #} &,
  dual[convexPolyG, 0.85],
  {1, 2}]
 ]

Face vertices of Polyhedron objects are not actually coplanar?

2 Answers2