Difference (or intersection) of two convex polyhedra

Question

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

I haven't looked at it myself, but there is some stuff here. — wxffles, Jul 2, 2014 at 21:12
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... — Wouter, Jul 2, 2014 at 22:28
@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. — level1807, Jul 3, 2014 at 5:23
@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. — level1807, Jul 3, 2014 at 5:34

Junho Lee · Accepted Answer · 2015-06-10 22:40:16Z

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

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? — level1807, Jul 3, 2014 at 7:33
I replaced the second picture, so it's OK. Still, constructing the final intersection/difference is interesting. — level1807, Jul 3, 2014 at 7:45
@Juhno Lee...thank you for this excellent answer and the version 10 update...I learned a lot. I had already up voted original answer.:) — ubpdqn, Jun 11, 2015 at 0:13
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). — level1807, Jun 14, 2015 at 18:34

ubpdqn · Accepted Answer · 2015-06-08 01:26:45Z

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

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? — level1807, Jun 7, 2015 at 13:02
@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... — ubpdqn, Jun 7, 2015 at 13:05

Stack Exchange Network

Difference (or intersection) of two convex polyhedra

2 Answers 2

Your Answer

Not the answer you're looking for? Browse other questions tagged
computational-geometry
intersection
polyhedra
or ask your own question.

Linked

Hot Network Questions

Difference (or intersection) of two convex polyhedra

2 Answers 2

Your Answer

Sign up or log in

Post as a guest

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

Linked

Related

Hot Network Questions

Not the answer you're looking for? Browse other questions tagged
computational-geometry
intersection
polyhedra
or ask your own question.