Generating convex polyhedron from face planes?

Question

Suppose I have lists of normals and points for planes. There's a convex polyhedron whose faces lie on these planes and are bounded by plane intersections. What would be the easiest way to produce an image of this polyhedron (preferably with the vertices known) in Mathematica?

I thought of one really ugly way to do it, but I'm looking for something more automated/efficient than:

1. Iterate through all triples of planes looking to see if they have a common intersection point, and throwing those intersection points into a list.
2. Using a Mathematica command to strip duplicates from the list.
3. ListPointPlot3D to display those points so I can manually decide which should be grouped together in a polygonal face.
4. Using a Mathematica command to build the polyhedron from properly oriented lists of vertices for faces.

Edit: As requested, here's an example: Suppose I have

normals = {{0, 0, 1}, {0, -2, 2}, {2, 0, 2}, {0, 2, 2}, {-2, 0, 2}}

and

pts = {{1, 1, 0}, {0, 0, 1}, {0, 0, 1}, {0, 0, 1}, {0, 0, 1}}

Then I can make a bunch of plane equations with, say

Table[Table[
  Dot[normals[[i]], {x, y, z} - pts[[i]]], {i, 1, 
  Dimensions[normals][[1]]}][[j]] == 0, {j, 1, 
  Dimensions[normals][[1]]}]

Three of these five planes intersect at a single point at exactly five points (which could be calculated in a number of ways): those in {{0, 0, 1}, {1, 1, 0}, {1, -1, 0}, {-1, 1, 0}, {-1, -1, 0}}. The final result I'd be looking for in this case would be something like

Graphics3D[
 Polygon[{{{0, 0, 1}, {1, -1, 0}, {1, 1, 0}}, {{0, 0, 1}, {1, 1, 
 0}, {-1, 1, 0}}, {{0, 0, 1}, {-1, 1, 0}, {-1, -1, 0}}, {{0, 0, 
 1}, {-1, -1, 0}, {1, -1, 0}}, {{1, 1, 0}, {-1, 1, 0}, {-1, -1, 
 0}, {1, -1, 0}}}]]

Answer 1

Because (a) RegionPlot3D does not render edges well and (b) detailed information about the vertices and faces could be worthwhile, I will offer a solution that finds this information and displays it clearly. (The first two lines of code produce the region plot, if you just want to stop there; the rest develop the improved solution.)

I am stuck at one thing: it is hard to find an efficient algorithm to determine the proper orientation of the normals. When you're just given a bunch of planes, they partition space into lots of polytopes. Some of those will be unbounded, so they can be neglected, but potentially there are many bounded polytopes. We could assume exactly one of them contacts every one of the planes nontrivially: this gives a criterion for finding the polytope that is being described. But this description, if carried out naively (e.g., through a brute-force examination) takes $2^N$ operations for $N$ planes, which is highly unsatisfactory except for small problems.

The following solution identifies a polytope by finding all mutual intersections of the planes (the "vertices"), then optionally reorienting each normal so that the number of vertices behind it is at least as great as the number of vertices in front of it. Although this does not always work, it may be of some service. Otherwise, if all normals are given in the proper (outward) orientations in the input, one can just delete the single line of code that does the re-orientation and get what was intended.

Step by step description

I will take you through the procedure step by step; the full listing is at the end. Begin with the data: parallel arrays of normals and points on the planes they define. The normals need to point constistently outward or inward of the polyhedron.

normals = {{0, 0, -1}, {0, -2, 2}, {2, 0, 2}, {0, 2, 2}, {-2, 0, 2}};
pts = {{1, 1, 0}, {0, 0, 1}, {0, 0, 1}, {0, 0, 1}, {0, 0, 1}};
dataGraphics = Graphics3D[{PointSize[0.015], Gray, Point[pts], Black, Thick, Arrowheads[Medium], 
  MapThread[Arrow[{#2, #1 + #2}] &, {normals, pts}]}]

It will be expeditious to exploit Mathematica's fast, compact matrix operations. Anticipating this, I represent each plane $(n_1,n_2,n_3)\cdot(x,y,z) = p$ as the four-vector $(-n_1,-n_2,-n_3,p)$:

planes = Union[MapThread[Append[-#1, #1.#2] &, {normals, pts}]]

At this point we can easily see the polyhedron by means of RegionPlot:

regionGraphic = RegionPlot3D[Min[planes . {x, y, z, 1}] >= 0, {x, -1, 1}, {y, -1, 1}, {z, 0, 1}, 
  PlotPoints -> 50, BoxRatios -> {1, 1, 1/2}, Mesh -> None, PlotStyle -> Opacity[0.85]]

We will eventually need to inspect all mutual intersections, which are obtained by taking planes three at a time:

nodes = Union[Append[#, 1] & /@ 
          Quiet[Cases[LinearSolve[Most /@ #, -Last /@ #] & /@ Subsets[planes, {3}], _List]]]

{{-1, -1, 0, 1}, {-1, 1, 0, 1}, {0, 0, 1, 1}, {1, -1, 0, 1}, {1, 1, 0, 1}}

Quiet suppresses messages when LinearSolve finds no intersections. The reason for appending $1$ to each node is that the oriented distance of a point $(x,y,z)$ from any plane given in the form $(-n_1,-n_2,-n_3,p)$ is proportional to the inner product $(x,y,z,1)\cdot(-n_1,-n_2,-n_3,p)$: this is what made the RegionPlot application so easy.

To reorient the normals (which is optional) we only have to count the signs of the oriented distances of every node:

planes = MapThread[Times, {planes, 2 UnitStep[Total[nodes . #]] - 1 & /@ planes}]

{{-2, 0, -2, 2}, {0, -2, -2, 2}, {0, 0, 1, 0}, {0, 2, -2, 2}, {2, 0, -2, 2}}

The vertices of the polytope, then, are those behind every plane, with some allowance for numerical imprecision:

vertices = Select[nodes, Chop[Min[planes.#]] >= 0 &];

(See the comments concerning this expression.)

The next steps assemble the vertices into faces in a form suitable for a GraphicsComplex. To do this, we first create the vertex-face incidence matrix, once again allowing for some imprecision, and square it to obtain the vertex-vertex adjacency matrix:

incidence = SparseArray[Outer[Boole[Chop[#1.#2] == 0] &, vertices, planes, 1]];
adjacency = Map[Boole[# >= 2] & , incidence . incidence\[Transpose], {2}];

The adjacency matrix determines the vertex graph:

At this point we can exploit Mathematica's graph algorithms to tie these vertices into polygons to represent the faces: it's a question of (easily) finding their order around each face.

faceNodes = Flatten[Position[# // Normal, 1]] & /@ (incidence\[Transpose]);
faceGraphs = (SimpleGraph[AdjacencyGraph[adjacency[[#, #]]]] & /@ faceNodes);
orderings = First /@ First[FindEulerianCycle[#]] &  /@ faceGraphs;
faces = MapThread[Part, {faceNodes, orderings}]

{{3, 5, 4}, {2, 5, 3}, {1, 4, 5, 2}, {1, 4, 3}, {1, 3, 2}}

The display is now easy:

polyGraphics = Graphics3D[{GraphicsComplex[Most /@ vertices, 
  {Opacity[0.5], Polygon[faces], 
   PointSize[0.015], Red, Opacity[1], Point[Range[Length[vertices]]]}]}];
Show[dataGraphics, polyGraphics, Boxed -> False]

It works pretty well on more complex convex polytopes, too. Here's one with 50 facets, found in 5 seconds:

---

Complete code listing

polyhedron[normals_, pts_] := 
 Module[{planes, nodes, vertices, incidence, adjacency, faceNodes, 
   faceGraphs, orderings, faces, result},
  planes = Union[MapThread[Append[-#1, #1.#2] &, {normals, pts}]];
  nodes = 
   Union[Append[#, 1] & /@ 
     Quiet[Cases[LinearSolve[Most /@ #, -Last /@ #] & /@ Subsets[planes, {3}], _List]]];
  (* planes = MapThread[Times, {planes, 2 UnitStep[Total[nodes . #]] - 1& /@ planes}];*)
  vertices = Select[nodes, Chop[Min[planes.#]] >= 0 &];
  incidence = SparseArray[Outer[Boole[Chop[#1.#2] == 0] &, vertices, planes, 1]];
  adjacency = Map[Boole[# >= 2] & , incidence . incidence\[Transpose], {2}];
  faceNodes = Flatten[Position[# // Normal, 1]] & /@ (incidence\[Transpose]);
  faceGraphs = (SimpleGraph[AdjacencyGraph[adjacency[[#, #]]]] & /@ faceNodes);
  orderings = First /@ First[FindEulerianCycle[#]] &  /@ faceGraphs;
  faces = MapThread[Part, {faceNodes, orderings}];
  result["vertices"] = Most /@ vertices; result["faces"] = faces;
  result
  ];
dataGraphics = Graphics3D[{PointSize[0.015], Blue, Point[pts], Blue, Arrowheads[Medium], 
  MapThread[Arrow[{#2, #1 + #2}] &, {normals, pts}]}];
p = polyhedron[normals, pts]; vertices = p["vertices"]; faces = p["faces"]; v = Length[vertices];
polyGraphics = Graphics3D[{GraphicsComplex[vertices, 
 {Opacity[0.5], Polygon[faces], 
  PointSize[0.015], Red, Opacity[1], Point[Range[v]]}]}];
Show[dataGraphics, polyGraphics, Boxed -> False]

Answer 2

Thanks to everyone for their helpful comments. I think I have patched together a solution that does what I want (using the package @Daniel Lichtblau linked and the ideas in others' comments). I demonstrate it here with the example from my question:

Get["C:\\VertexEnum.m"]; (*Change path as appropriate.*)
(*If you don't use semicolons, you may want: 
Off[General::spell];  Off[General::spell1]; 
Off[Syntax::com]; Off[Write::noopen]; *)
normals = {{0, 0, -1}, {0, -2, 2}, {2, 0, 2}, {0, 2, 2}, {-2, 0, 2}}; 
(*The normals must be outwardly oriented*)
pts = {{1, 1, 0}, {0, 0, 1}, {0, 0, 1}, {0, 0, 1}, {0, 0, 1}}; 
shift = {10, 10, 10};
(*The shift has to be large enough to push the polyhedron into the first octant 
for the package to work, but its actual values don't matter.*)
ds = MapThread[#1.(#2 + shift) &, {normals, pts}]; 
{shiftedvs,activesets} = VertexEnum`VertexEnumeration[normals, ds]; 
vs = MapThread[#1 - shift &, {shiftedvs}]; 
poly = VertexEnum`Polyhedron3D[vs, activesets];
Show[Graphics3D[poly, 
 ListPointPlot3D[vs, PlotStyle -> PointSize[0.05]], 
 BoxRatios -> Automatic]

vs is a list of the vertices if needed for other things, poly is a list of lists of vertices to represent the faces, and the Show command generates something like this:

Edit: I realized it's better to let the package handle the making of the polyhedron (although possibly the polygons are oriented the wrong way). While something like RegionPlot3D[ And @@ MapThread[#1.({x, y, z} - #2) <= 0 &, {normals, pts}], {x, -2, 2}, {y, -2, 2}, {z, -1, 2}, PlotPoints -> 25] will technically work to roughly show the polyhedron, the newly-edited above code uses the package to extract the faces exactly.