21
$\begingroup$

I have two convex polyhedra stored in the following form: a set of vertices vertices = {{x1,y1,z1},...}, a set of faces, where each face is a convex polygon specified by the ordered list of the numbers of its vertices. For example,

vertices1 = {{-1, -1, -1}, {-1, -1, 1}, {-1, 1, -1}, {-1, 1, 1}, {1, -1, -1}, {1, -1, 1}, 
             {1, 1, -1}, {1, 1, 1}};
vertices2 = {{-(1/2), -(1/2), -(1/2)}, {5/6, -(7/6), 5/6}, {-(7/6), 5/6, 5/6}, 
             {1/6, 1/6, 13/6}, {5/6, 5/6, -(7/6)}, {13/6, 1/6, 1/6}, {1/6, 13/6, 1/6}, 
             {3/2, 3/2, 3/2}};
faces = {{5, 6, 8, 7}, {1, 2, 4, 3}, {3, 4, 8, 7}, {1, 2, 6, 5}, 
         {2, 4, 8, 6}, {1, 3, 7, 5}};

These are just two cubes, one of them rotated and translated. They can be visualized by GraphicsComplex:

Show[{Graphics3D@GraphicsComplex[vertices1, Polygon /@ faces], Graphics3D@GraphicsComplex[vertices2, Polygon /@ faces]}]

Mathematica graphics

I need to find a way to calculate the exact coordinates of the vertices for the concave polyhedron that is the difference of these two ($P_1\setminus P_2$ or $P_2\setminus P_1$), which, I guess, is almost equivalent to finding their intersection.

Obviously, in my example the difference consists of several polyhedra, not one, but the idea is still the same -- this is just a set of faces stored as lists of numbers of vertices. The output has to have the same form as the input. The algorithm has to be applicable to any pair of convex polyhedra.

Edit: one interesting case is when the polyhedra do not even contain any of each other's vertices:

`vertices1 = {{-1, -1, -1}, {-1, -1, 1}, {-1, 1, -1}, {-1, 1, 1}, {1, -1, -1}, {1, -1, 1}, {1, 1, -1}, {1, 1, 1}};`
`vertices2 = {{-(1/2) - 1/Sqrt[2], -(1/2) + 1/Sqrt[2],0}, {-(1/2) + 1/Sqrt[2], -(1/2) - 1/Sqrt[2], 0}, {1/2 - 1/Sqrt[2], 1/2 + 1/Sqrt[2], Sqrt[2]}, {1/2 + 1/Sqrt[2], 1/2 - 1/Sqrt[2], Sqrt[2]}, {1/2 - 1/Sqrt[2], 1/2 + 1/Sqrt[2], -Sqrt[2]}, {1/2 + 1/Sqrt[2], 1/2 - 1/Sqrt[2], -Sqrt[2]}, {3/2 - 1/Sqrt[2], 3/2 + 1/Sqrt[2], 0}, {3/2 + 1/Sqrt[2], 3/2 - 1/Sqrt[2], 0}};`
`Show[{Graphics3D@GraphicsComplex[vertices1, Polygon /@ faces], Graphics3D@GraphicsComplex[vertices3, Polygon /@ faces]}]`

enter image description here

$\endgroup$
  • 1
    $\begingroup$ I haven't looked at it myself, but there is some stuff here. $\endgroup$ – wxffles Jul 2 '14 at 21:12
  • $\begingroup$ Just some thoughts: if the convex hulls of both polyhedra P1 and P2 are known, then a first step would be to determine which vertices of P1 are outside of P2 and vice versa. This is easy using the sign of the circulation of vertices face by face. Next, for each rib having a vertex inside and the other outside of the other polyhedron, determine the intersection of the line (rib of P1) with the face of P2. With "the face", I mean the one separating the two vertices at the extremes of "the rib". Now, take the "outside vertex" and the "rib-face" intersections. Forgive my handwaving... $\endgroup$ – Wouter Jul 2 '14 at 22:28
  • $\begingroup$ @Wouter, thanks, sounds like a good plan! The biggest problem still is figuring out which face is "the face" for the given "rib". Especially if it's not actually a face, but an edge or even a vertex, like in my example. $\endgroup$ – level1807 Jul 3 '14 at 5:23
  • $\begingroup$ @Wouter moreover, I came up with an example when the polyhedra do not contain any of each other's vertices. See the Edit in the question. $\endgroup$ – level1807 Jul 3 '14 at 5:34
18
+50
$\begingroup$

Update for v10

I used MeshRegion and MeshPrimitives for intersected points.

linePoly[v1_, v2_, f_] := Module[{fC = Append[#, #[[1]]] & /@ f},
  {x, y, z} /. NSolve[
    Or@@ ({x, y, z}\[Element]# & /@ 
        MeshPrimitives[MeshRegion[v1, Line /@ fC], 1]) &&
     Or@@ ({x, y, z}\[Element]# & /@ 
        MeshPrimitives[MeshRegion[v2, Polygon /@ f], 2])]]
polysIntersected[v1_, v2_, f_] := Union[linePoly[v1, v2, f], linePoly[v2, v1, f]]

There are your data.

vertices3 = {{-1, -1, -1}, {-1, -1, 1}, {-1, 1, -1}, {-1, 1,1}, {1, -1, -1}, {1, -1, 1}, {1, 1, -1}, {1, 1, 1}};
vertices4 = {{-(1/2) - 1/Sqrt[2], -(1/2) + 1/Sqrt[2], 0}, {-(1/2) + 1/Sqrt[2], -(1/2) - 1/Sqrt[2], 0}, {1/2 - 1/Sqrt[2], 1/2 + 1/Sqrt[2], Sqrt[2]}, {1/2 + 1/Sqrt[2], 1/2 - 1/Sqrt[2], Sqrt[2]}, {1/2 - 1/Sqrt[2],  1/2 + 1/Sqrt[2], -Sqrt[2]}, {1/2 + 1/Sqrt[2], 1/2 - 1/Sqrt[2], -Sqrt[2]}, {3/2 - 1/Sqrt[2], 3/2 + 1/Sqrt[2],  0}, {3/2 + 1/Sqrt[2], 3/2 - 1/Sqrt[2], 0}};
faces = {{5, 6, 8, 7}, {1, 2, 4, 3}, {3, 4, 8, 7}, {1, 2, 6, 5}, {2,4, 8, 6}, {1, 3, 7, 5}};

You can test for second your data.

points = polysIntersected[vertices3, vertices4, faces] // Chop;
Show[Graphics3D[GraphicsComplex[#, {Opacity[0.5],  Polygon /@ faces}]] & /@ {vertices3, vertices4}];
Graphics3D[{PointSize[Large], Red, Point /@ points}];
Show[%, %%]

enter image description here

And you might want to know intersected points like this.

points
{{-1., 0, 0}, {-1., 0, 0}, {-1., 0.414214, -0.292893}, {-1., 0.414214,
   0.292893}, {-0.5, 0.914214, -1.}, {-0.5, 0.914214, 1.}, {-0.414214,
   1., 1.}, {-0.414214, 1., -1.}, {0, -1., 
  0}, {0.414214, -1., -0.292893}, {0.414214, -1., 
  0.292893}, {0.585786, 1., -1.}, {0.585786, 1., 
  1.}, {0.914214, -0.5, -1.}, {0.914214, -0.5, 1.}, {1., -0.414214, 
  1.}, {1., -0.414214, -1.}, {1., 0.585786, -1.}, {1., 0.585786, 
  1.}, {1., 1., -0.707107}, {1., 1., 0.707107}}

Old

Have try this following code.

vertices1 = {{-1, -1, -1}, {-1, -1, 1}, {-1, 1, -1}, {-1, 1, 
    1}, {1, -1, -1}, {1, -1, 1}, {1, 1, -1}, {1, 1, 1}};
vertices2 = {{-(1/2), -(1/2), -(1/2)}, {5/6, -(7/6), 5/6}, {-(7/6), 
    5/6, 5/6}, {1/6, 1/6, 13/6}, {5/6, 5/6, -(7/6)}, {13/6, 1/6, 
    1/6}, {1/6, 13/6, 1/6}, {3/2, 3/2, 3/2}};
faces = {{5, 6, 8, 7}, {1, 2, 4, 3}, {3, 4, 8, 7}, {1, 2, 6, 5}, {2, 
    4, 8, 6}, {1, 3, 7, 5}};

makeSurfaceEq[v_List, {x_, y_, z_}] :=
 Module[{a, b, c, d, eq},
   eq = a x + b y + c z + d == 0;
   eq /. Solve[(a #1 + b #2 + c #3 + d == 0 &) @@@ v, {a, b, c}][[
      1]] /. d -> 1
   ] /; Length[v] >= 3

polygonToSurfaceInEq[v_, f_, {x_, y_, z_}] := Module[{cOfG, eqs, ineq},
cOfG = (Plus @@ v)/Length[v];
eqs = makeSurfaceEq[v[[#]], {x, y, z}] & /@ f;
Table[
  ineq = eqs[[i]] /. Equal -> Greater;
  If[ineq /. Thread[{x, y, z} -> cOfG], ineq, ! ineq],
  {i, Length[eqs]}] /.
 {Greater -> GreaterEqual, Less -> LessEqual}
]

makeLineInEq[v_List, {x_, y_, z_}] :=
 Module[{eq, t},
   eq = (v[[1]] - (Subtract @@ v) t) == {x, y, z};
   Reduce[Exists[{t}, eq && 0 <= t <= 1], {x, y, z}, Reals]
   ] /; Length[v] == 2


polygonToLineInEq[v_, f_, {x_, y_, z_}] := Module[{eqs},
  eqs = v[[#]] & /@ f;
  eqs = Partition[#, 2, 1, 1] & /@ eqs;
  eqs = Union[Flatten[eqs, {1, 2}], 
    SameTest -> (Sort[#1] == Sort[#2] &)];
  eqs = makeLineInEq[#, {x, y, z}] & /@ eqs
  ]

polyInterSected[v1_, v2_, f_] := Module[
  {x, y, z, ieq1, ieq2, ieq3, ieqs, eqs1, eqs2},
  ieq1 = polygonToLineInEq[v1, f, {x, y, z}];
  ieq2 = polygonToSurfaceInEq[v2, f, {x, y, z}];
  ieq3 = ieq2 /. {Greater | GreaterEqual | Less | LessEqual -> Equal};
  Off[NSolve::svars];
  p1 = Cases[{x, y, z} /. 
     NSolve[(Or @@ ieq1) && (And @@ ieq2) && (Or @@ ieq3), {x, y, 
       z}], {_?NumberQ, _?NumberQ, _?NumberQ}]
  ]

polyInterSected[vertices1, vertices2, faces]

{{-1., 0.75, 1.}, {-1., 1., 0.75}, {-0.5, 1., 1.}, {0.75, -1., 1.}, {0.75, 1., -1.}, {1., -1., 0.75}, {1., -0.5, 1.}, {1., 0.75, -1.}, {1., 1., -0.5}}

g1 = Show[{
    Graphics3D[GraphicsComplex[vertices1,
      {Opacity[0.5], Polygon /@ faces}]],
    Graphics3D[GraphicsComplex[vertices2,
      {Opacity[0.5], Polygon /@ faces}]]}];
g2 = Graphics3D[{PointSize[Large], Red, 
    Point /@ polyInterSected[vertices1, vertices2, faces]}];
Show[g1, g2]

enter image description here

$\endgroup$
  • $\begingroup$ Very nice! Something's off with the second picture though :) And can you find all the vertices of the intersection/difference, i.e. the internal ones as well? $\endgroup$ – level1807 Jul 3 '14 at 7:33
  • $\begingroup$ I replaced the second picture, so it's OK. Still, constructing the final intersection/difference is interesting. $\endgroup$ – level1807 Jul 3 '14 at 7:45
  • $\begingroup$ @Juhno Lee...thank you for this excellent answer and the version 10 update...I learned a lot. I had already up voted original answer.:) $\endgroup$ – ubpdqn Jun 11 '15 at 0:13
  • $\begingroup$ Thanks, this looks very good. And I'm very pleased with how v10 can handle these tasks with regions. I gave you the bounty and I hope that either I or someone else will be able to work on this and construct the faces of the final intersection (which is simply the convex hull of points) or difference (probably more tricky). $\endgroup$ – level1807 Jun 14 '15 at 18:34
4
$\begingroup$

This takes some time for pre-processing:

h1 = Hexahedron[{{-1.1666666666666667`, 0.8333333333333334`, 
    0.8333333333333334`}, {0.16666666666666666`, 2.1666666666666665`, 
    0.16666666666666666`}, {1.5`, 1.5`, 1.5`}, {0.16666666666666666`, 
    0.16666666666666666`, 
    2.1666666666666665`}, {-0.5`, -0.5`, -0.5`}, {0.8333333333333334`,
     0.8333333333333334`, -1.1666666666666667`}, {2.1666666666666665`,
     0.16666666666666666`, 
    0.16666666666666666`}, {0.8333333333333334`, -1.1666666666666667`,
     0.8333333333333334`}}]; h2 = 
 Hexahedron[{{-1.`, 1.`, -1.`}, {1.`, 1.`, -1.`}, {1.`, 1.`, 
    1.`}, {-1.`, 1.`, 
    1.`}, {-1.`, -1.`, -1.`}, {1.`, -1.`, -1.`}, {1.`, -1.`, 
    1.`}, {-1.`, -1.`, 1.`}}];
ri = RegionIntersection[h1, h2];
rd1 = RegionDifference[h1, ri];
rd2 = RegionDifference[h2, ri];

Visualizing:

With[{rp1 = RegionPlot3D[ri, PlotPoints -> 100, PlotStyle -> Red], 
  rp2 = RegionPlot3D[rd1, PlotPoints -> 100, PlotStyle -> Green], 
  rp3 = RegionPlot3D[rd2, PlotPoints -> 100, PlotStyle -> Blue]}, 
 Manipulate[
  Show[If[p == 1, rp1, Graphics3D[]],
   If[q == 1, rp2, Graphics3D[]],
   If[r == 1, rp3, Graphics3D[]], Boxed -> False, Axes -> False, 
   Background -> Black, 
   PlotRange -> Table[{-3, 3}, {3}]], {p, {0, 1}}, {q, {0, 
    1}}, {r, {0, 1}}]]

enter image description here

An estimate of volume of intersection: Volume[DiscretizeRegion@ri] is 3.61699

Update

To find points of intersection:

pts = Part[vertices1, #] & /@ faces;
pts2 = Part[vertices2, #] & /@ faces;
ip1 = InfinitePlane[#[[1 ;; 3]]] & /@ pts;
ip2 = InfinitePlane[#[[1 ;; 3]]] & /@ pts2;
ans = Cases[
   RegionIntersection @@@ 
    Tuples[RegionIntersection @@@ Tuples[{ip1, ip2}], 2], Point[x_]];
rmfun1[x_] := Or @@ Through[(RegionMember /@ (Polygon /@ pts))[x]]
rmfun2[x_] := Or @@ Through[(RegionMember /@ (Polygon /@ pts2))[x]]
rmf[x_] := And[rmfun1[x], rmfun2[x]]
pck = Union[Pick[ans, rmf /@ ans[[All, 1]]]];
Graphics3D[{Red, PointSize[0.02], pck, , Blue, Opacity[0.5], 
  Polygon /@ pts, Yellow, Polygon /@ pts2}]

enter image description here

Also works for second set of vertices (but no efficient):

enter image description here

with points:

{Point[{-1, 0, 0}], Point[{-(1/2), -(1/2) + Sqrt[2], -1}], 
 Point[{-(1/2), -(1/2) + Sqrt[2], 1}], Point[{0, -1, 0}], 
 Point[{1, 1 - Sqrt[2], -1}], 
 Point[{1, 1 - Sqrt[2], 
   1 - Sqrt[2] - (3 - 2 Sqrt[2])/Sqrt[2] + 1/2 (-4 + 5 Sqrt[2])}], 
 Point[{1, 2 - Sqrt[2], -1}], Point[{1, 2 - Sqrt[2], 1}], 
 Point[{1 - Sqrt[2], 1, 
   1 - 3/Sqrt[2] - Sqrt[2] + 1/2 (-4 + 5 Sqrt[2])}], 
 Point[{-(1/2) + Sqrt[2], -(1/2), -1}], 
 Point[{-(1/2) + 1/2 (3 - 2 Sqrt[2]), -(1/2) + Sqrt[2] + 
    1/2 (3 - 2 Sqrt[2]), 1}], 
 Point[{1/2 + (2 - Sqrt[2])/Sqrt[2], -(3/2) + Sqrt[2] - (2 - Sqrt[2])/
    Sqrt[2], 1}], 
 Point[{1/2 (1 + Sqrt[2] + (-2 + Sqrt[2])/Sqrt[2]), 
   1 + 1/2 (1 - Sqrt[2] - (-2 + Sqrt[2])/Sqrt[2]), 
   1/Sqrt[2] - Sqrt[2]}], 
 Point[{1/2 (1 + Sqrt[2] + (-2 + Sqrt[2])/Sqrt[2]), 
   1 + 1/2 (1 - Sqrt[2] - (-2 + Sqrt[2])/Sqrt[2]), -(1/Sqrt[2]) + 
    Sqrt[2]}], 
 Point[{1/2 (2 - Sqrt[2]) (1 + Sqrt[2] + (-2 + Sqrt[2])/Sqrt[2]), 
   1 + 1/2 (2 - Sqrt[2]) (1 - Sqrt[2] - (-2 + Sqrt[2])/Sqrt[2]), 
   Sqrt[2] - (2 - Sqrt[2])/Sqrt[2]}], 
 Point[{1/2 (2 - Sqrt[2]) (1 + Sqrt[2] + (-2 + Sqrt[2])/Sqrt[2]), 
   1 + 1/2 (2 - Sqrt[2]) (1 - Sqrt[2] - (-2 + Sqrt[2])/Sqrt[
       2]), -Sqrt[2] + (2 - Sqrt[2])/Sqrt[2]}], 
 Point[{1/2 (-1 + Sqrt[2]) (1 + Sqrt[2] + (-2 + Sqrt[2])/Sqrt[
      2]), -1 + 
    1/2 (-1 + Sqrt[2]) (1 - Sqrt[2] - (-2 + Sqrt[2])/Sqrt[2]), -((-1 +
      Sqrt[2])/Sqrt[2])}], 
 Point[{1/2 (-1 + Sqrt[2]) (1 + Sqrt[2] + (-2 + Sqrt[2])/Sqrt[
      2]), -1 + 
    1/2 (-1 + Sqrt[2]) (1 - Sqrt[2] - (-2 + Sqrt[2])/Sqrt[2]), (-1 + 
    Sqrt[2])/Sqrt[2]}], 
 Point[{-1 + 
    1/2 (-1 + Sqrt[2]) (1 + Sqrt[2] + (-2 - Sqrt[2])/Sqrt[2]), 
   1/2 (-1 + Sqrt[2]) (1 - Sqrt[2] - (-2 - Sqrt[2])/Sqrt[2]), -((-1 + 
     Sqrt[2])/Sqrt[2])}], 
 Point[{-1 + 
    1/2 (-1 + Sqrt[2]) (1 + Sqrt[2] + (-2 - Sqrt[2])/Sqrt[2]), 
   1/2 (-1 + Sqrt[2]) (1 - Sqrt[2] - (-2 - Sqrt[2])/Sqrt[2]), (-1 + 
    Sqrt[2])/Sqrt[2]}]}
$\endgroup$
  • $\begingroup$ This may be enough for some simple visualization. However, I specifically stated that the result needs to be analytic, i. e. a list of vertices and faces. Can this information be extracted from RegionDifference? $\endgroup$ – level1807 Jun 7 '15 at 13:02
  • $\begingroup$ @level1807 yes apologies...I only foolishly read title and looked at graphics...and post fact realized the aim was intersection points. I have upvoted Juhno Lee's answer... $\endgroup$ – ubpdqn Jun 7 '15 at 13:05

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.