Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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}}}]]
share|improve this question
1  
As always, some example data is very helpful. It's probably going to take me a few minutes to come up with a simple and meaningful example to even begin to attempt possible solutions and I'm already bored by that point and move on to something else... The onus is on you, the OP, to put forth a minimal example that saves us the trouble. –  rm -rf Nov 14 '12 at 15:23
2  
Since you can build the Hessian normal form of your planes, turn them into the corresponding inequalities and feed them into RegionPlot3D[]. –  J. M. Nov 14 '12 at 15:26
2  
@rm-rf, thank you for the comment. I hope the example I've edited in clears things up. Let me know if there's anything else I can do to make this question more accessible. J.M., are you saying to use the inequalities for "on the side of the plane towards the center" in each case? If so, is there an easy way to do that automatically if I have too many planes to handle manually? –  Mark S. Nov 14 '12 at 15:50
2  
Are your normals oriented, that is, do they consistently point outward? They're not in your example data, but maybe they are in the real data you want to work with. If so, you can do the RegionPlot3D thing. If not, dealing with unoriented planes can be quite problematic (see my comment on @whuber's answer). –  Rahul Narain Nov 14 '12 at 23:32
2  
There is a package for polytope vertex enumeration, coauthored by Komei Fukuda. –  Daniel Lichtblau Nov 15 '12 at 16:21
show 3 more comments

2 Answers

up vote 8 down vote accepted

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.

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

Points and normals

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

Region plot

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] &];

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:

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, vertexGraphics, polyGraphics, Boxed -> False]

Final

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

Random polytope


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]
share|improve this answer
    
You attempted to tackle the issue of poor plane orientation and you have a better presentation of the final result. Thanks for your effort. –  Mark S. Nov 16 '12 at 22:37
    
Your answers are terrific! –  cormullion Nov 16 '12 at 23:07
add comment

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: pyramid

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.

share|improve this answer
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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