Generating convex polyhedron from face planes?
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]
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.