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$
4
  • 1
    $\begingroup$ I haven't looked at it myself, but there is some stuff here. $\endgroup$
    – wxffles
    Jul 2, 2014 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, 2014 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, 2014 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, 2014 at 5:34

2 Answers 2

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$
4
  • $\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, 2014 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, 2014 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, 2015 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, 2015 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$
2
  • $\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, 2015 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, 2015 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.