Intersecting ellipsoids

Question

First of all I know that finding whether two ellipsoids intersect is complicated and there are specific algorithms created for this purpose (see, for instance, [1]. Actually, I am trying to follow this documentation in order to implement something in Mathematica. Of course there are other approaches).

I have "googled" for a while, searching for Mathematica codes but I could not find anything relevant. Maybe this has to do with the fact that functions like RegionIntersection were introduced rather recently (2014; version 10.0).

Anyway, here is a simple code. It's not panacea of course. (Found here)

RegionIntersectQ[e1_, e2_] := 
 RegionQ[DiscretizeRegion[RegionIntersection[e1, e2]] // Quiet]

It works in many cases. E.g.

Graphics3D[{g1 = 
   Ellipsoid[{0, 0, 0}, {{5, 2, 3}, {2, 3, 2}, {3, 2, 5}}], 
  g2 = Ellipsoid[{1, 1, 1}, {{10, 2, 3}, {2, 3, 2}, {3, 2, 5}}]}, 
 Axes -> True]
RegionIntersectQ[g1, g2] // Timing
(*2.62...,True*)

My first question is about the time efficiency. Why 2.6 sec are needed?

We can find (in version 10.3 at least!) other cases that the code finds correctly whether the two ellipsoids intersect (the code tries harder when there is a tiny region of intersection).

But there are cases where it fails strangely and I do not understand why.

Graphics3D[{{Opacity[0.6], 
   g1 = Ellipsoid[{0, 0, 0}, {{5, 2, 3}, {2, 3, 2}, {3, 2, 5}}]}, 
  g2 = Ellipsoid[{4.4, 1, 1}, {{10, 2, 3}, {2, 3, 2}, {3, 2, 7}}]}, 
 Axes -> True]
RegionIntersectQ[g1, g2] // Timing

I am not expecting other people to spend their time to implement algorithms of others but I am wondering how we can use efficient the current functionality of Mathematica (built-in functions like RegionIntersection) to check if two ellipsoids intersect or not.

EDIT After writing the question I found this which has connection with my previous question here. @JasonB addresses the issue very thoroughly.

But the query remains the same: Can I expect Mathematica's relevant functionality to do the job (even with a bit help of us) or should I forget the whole idea?

Just checking, but you do know RegionIntersection is not implemented for MeshRegion objects embedded in 3D? — M.R., Dec 04 '15 at 09:05
@M.R. I don't see your point. Can you elaborate more on your comment? Thanks. — Dimitris, Dec 04 '15 at 09:13
Check the docs for ref/RegionIntersection under Possible issues, and it says RegionIntersection functionality in 3d is not complete. — M.R., Dec 04 '15 at 09:18
@M.R. Thanks again. Very interesting. I read a lot of times this documentation page but I didn't pay attention to this! But an Ellipsoid seems less complicated than DelaunayMesh[RandomReal[1, {20, 3}]] :-)! — Dimitris, Dec 04 '15 at 09:22
I agree that the reason for the false negative result is the RegionIntersection, but the reason for the timing problem is definitely DiscretizeRegion, evaluate: First /@ {step1 = RegionIntersection[g1, g2]; // AbsoluteTiming, step2 = DiscretizeRegion[step1]; // AbsoluteTiming, step3 = RegionQ[step2]; // AbsoluteTiming} — Jason B., Dec 04 '15 at 10:12
BTW, @dimitris, what are you working on with all these graphics problems? They are really quite fun and I'm tempted to blow off a days work to implement the algorithms you linked to. — Jason B., Dec 04 '15 at 10:13
@JasonB. You are right. But as I said I was wondering if we can use the current functionality of Mathematica for a simple test. — Dimitris, Dec 04 '15 at 10:16
This is an active research in computational geometry. Do you want to just visualize the ellipsoids and the intersections, or do you want to obtain the equation representing the surface of intersection. Also, do you want to determine if there is an intersection, or the input is always intersecting ellipsoids? You can also post this question in the Mathematics community. — Marvin, Dec 05 '15 at 09:41
@Saurav: I want to determine if there is intersection. No need for analytic expressions. — Dimitris, Dec 07 '15 at 14:38

score 3 · Accepted Answer · answered Jan 03 '17 at 01:12

FindInstance is faster:

FindInstance[RegionMember[g1, {x, y, z}] && RegionMember[g2, {x, y, z}],
 {x, y, z}] // AbsoluteTiming
(*  {0.002422, {{x -> 0, y -> 0, z -> 0}}}  *)

FindInstance[
  RegionMember[
    Ellipsoid[{5, 0, 0}, {{5, 2, 3}, {2, 3, 2}, {3, 2, 5}}], {x, y, z}] && RegionMember[
    Ellipsoid[{1, 1, 1}, {{10, 2, 3}, {2, 3, 2}, {3, 2, 5}}], {x, y, z}],
  {x, y, z}] // AbsoluteTiming
(*  {0.045233, {}}  *)

Also, Reduce is faster than the OP's approximate discretizing, but a bit slower than the above in the positive (first) case.

Reduce[Exists[{x, y, z}, RegionMember[g1, {x, y, z}] && RegionMember[g2, {x, y, z}]],
 {x, y, z}] // AbsoluteTiming
(*  {0.011632, True}  *)

Reduce[Exists[{x, y, z},
   RegionMember[Ellipsoid[{5, 0, 0}, {{5, 2, 3}, {2, 3, 2}, {3, 2, 5}}], {x, y, z}] &&
   RegionMember[Ellipsoid[{1, 1, 1}, {{10, 2, 3}, {2, 3, 2}, {3, 2, 5}}], {x, y, z}]],
 {x, y, z}] // AbsoluteTiming
(*  {0.044602, False}  *)

I saw the answer right now. Thank you very much! – Dimitris Jan 14 '17 at 12:15 — Dimitris, Jan 14 '17 at 12:15

score 3 · Answer 2 · answered Oct 05 '18 at 03:29

Starting in M11.1 you can use RegionDisjoint. Here are the two OP examples:

g1 = Ellipsoid[{0, 0, 0}, {{5, 2, 3}, {2, 3, 2}, {3, 2, 5}}];
g2 = Ellipsoid[{1, 1, 1}, {{10, 2, 3}, {2, 3, 2}, {3, 2, 5}}];

RegionDisjoint[g1, g2] //AbsoluteTiming

{0.002053, False}

g1 = Ellipsoid[{0, 0, 0}, {{5, 2, 3}, {2, 3, 2}, {3, 2, 5}}];
g2 = Ellipsoid[{4.4, 1, 1}, {{10, 2, 3}, {2, 3, 2}, {3, 2, 7}}];

RegionDisjoint[g1, g2] //AbsoluteTiming

{0.001416, False}

And here is @Michael's disjoint example:

r1 = Ellipsoid[{5, 0, 0}, {{5, 2, 3}, {2, 3, 2}, {3, 2, 5}}];
r2 = Ellipsoid[{1, 1, 1}, {{10, 2, 3}, {2, 3, 2}, {3, 2, 5}}];

RegionDisjoint[r1, r2] //AbsoluteTiming

{0.041539, True}

Intersecting ellipsoids

2 Answers2