4
$\begingroup$

I am trying to obtain a list of coordinates at which two spheres intersect.

Take for example the spheres: sp1=Sphere[{0, 0, 0}, 1] and sp2=Sphere[{1, 1, 1}, 1.5] When I Plot them it is clear they intersect, but i cannot retrieve the coordinates.

So far I Tried finding them by looking for the nearest or RegionNearest between RandomPoint[Sphere[{0, 0, 0}, 1], 1000] and RandomPoint[Sphere[{1, 1, 1}, 1.5],1000]. Nor did it work by looking for Intersection of the two list (since no coordinates actually matched) I also tried to look for the minimal euclidean distance between all points and the positions of those minimal distances, but in that case my computer crashes. Can somebody help? Thanks in advance

$\endgroup$
  • $\begingroup$ This seems simple enough to compute analytically (pen and paper) that it might not be worth bothering with the numerical approach ... $\endgroup$ – Szabolcs Jul 20 '16 at 12:05
  • 1
    $\begingroup$ I suggest this interesting read: Sphere-Sphere Intersection $\endgroup$ – DPF Jul 20 '16 at 12:07
  • $\begingroup$ it takes just a little tweaking to generalize that mathworld solution to the case where the spheres are not along a coordinate axis.. $\endgroup$ – george2079 Jul 20 '16 at 12:20
  • $\begingroup$ In my example it is simple enough for pen and paper, but i want to do it for 180 different and more complicated spheres, which makes it useful to be able to let mathematica take care of it. And I had seen the sphere sphere intersection explanation and notebook, but dont see how i can apply it to my spheres. $\endgroup$ – Seppen Jul 20 '16 at 12:24
  • $\begingroup$ What do you mean by "a list of coordinates" for the intersection? Just sampling the region? $\endgroup$ – kirma Jul 20 '16 at 14:47
6
$\begingroup$
sp1 = Sphere[{0, 0, 0}, 1];
sp2 = Sphere[{1, 1, 1}, 1.5];
ri1 = RegionIntersection[sp1, sp2];
l = MeshCoordinates[DiscretizeRegion[ri1]];
Show[Graphics3D[{Opacity[0.5], sp1, sp2, Thick, Red, 
   Line@l[[Last[FindShortestTour[l]]]]}]]

intersecting spheres

Note that you can find the center and radius of the circle with:

o = Mean[l];
r = Mean[Table[Norm[l[[i]] - Mean[l]], {i, Length[l]}]];

The normal of the circle is naturally in the direction of either of the centers of the spheres.

$\endgroup$
  • $\begingroup$ +1 - Now the question is, can you reorder the points so that they can be used in a Line, via Graphics3D[{Opacity[0.5], sp1, sp2, Red, Thick, Line@l}]? $\endgroup$ – Jason B. Jul 20 '16 at 14:06
  • $\begingroup$ @JasonB added the center and radius of the circle. Unfortunately Circle[] doesn't work in 3D. $\endgroup$ – Feyre Jul 20 '16 at 14:48
  • 1
    $\begingroup$ @Feyre - it seems in this example you can use FindShortestTour: This works: Graphics3D[{Opacity[0.5], sp1, sp2, Red, Thick, Line@l[[Last@FindShortestTour@l]]}] $\endgroup$ – Jason B. Jul 20 '16 at 14:50
  • $\begingroup$ @JasonB Yes, updated the answer $\endgroup$ – Feyre Jul 20 '16 at 15:02
4
$\begingroup$

This answer uses the circle3D function written by Taiki, found here. But it specifically does not use any of the Region functionality to find the intersection. Instead, I just grabbed some formulas from this MathWorld page

circle3D[centre_: {0, 0, 0}, radius_: 1, normal_: {0, 0, 1}, angle_: {0, 2 Pi}] :=
  Composition[
    Line,
    Map[RotationTransform[{{0, 0, 1}, normal}, centre], #] &,
    Map[Append[#, Last@centre] &, #] &,
    Append[DeleteDuplicates[Most@#], Last@#] &,
    Level[#, {-2}] &,
    MeshPrimitives[#, 1] &,
    DiscretizeRegion,
    If
  ][
    First@Differences@angle >= 2 Pi,
    Circle[Most@centre, radius],
    Circle[Most@centre, radius, angle]
  ]

sphereIntersection[{x1_, r1_}, {x2_, r2_}] := 
  Module[{dvec, d, d1, r},
   dvec = x2 - x1;
   d = Norm@dvec;
   d1 = (d^2 + r1^2 - r2^2)/(2 d);
   circle3D[x1 + d1 Normalize[x2 - x1], Sqrt[r1^2 - d1^2], 
    Normalize[x2 - x1]]
   ];
sphereIntersection[s1_Sphere, s2_Sphere] := 
 sphereIntersection[List @@ s1, List @@ s2]

Using the OP's spheres, you can visualize it via

Graphics3D[{Opacity[.3], sp1, sp2, Opacity[1], Blue, Thick, 
  sphereIntersection[sp2, sp1]}]

Mathematica graphics

You can grab the points via

Cases[Graphics3D[sphereIntersection[sp2, sp1]], 
 Line[{a__}] :> a, Infinity]

or just get random points via

RandomPoint[sphereIntersection[sp2, sp1]]
$\endgroup$
4
$\begingroup$

Using the results from this MathWorld page:

BlockRandom[SeedRandom["spherical"];
            c1 = RandomReal[1, 3]; c2 = RandomReal[1, 3];
            r2 = RandomReal[3/2]; r1 = RandomReal[3/2];
            d = EuclideanDistance[c1, c2];
            u = (d^2 + r1^2 - r2^2)/(2 d^2); cc = {1 - u, u}.{c1, c2}; 
            rc = Sqrt[r1^2 - d^2 u^2];
            Graphics3D[{Opacity[1/2], Sphere[c1, r1], Sphere[c2, r2],
                        {Red, Tube[BSplineCurve[
                                   AffineTransform[{RotationMatrix[{{0, 0, 1},
                                                    Normalize[c2 - c1]}], cc}][
                                   rc PadRight[{{1, 0}, {1, 1}, {0, 1}, {-1, 1},
                                                {-1, 0}, {-1, -1}, {0, -1}, {1, -1},
                                                {1, 0}}, {Automatic, 3}]], 
                                   SplineClosed -> True, SplineDegree -> 2, 
                                   SplineKnots -> {0, 0, 0, 1/4, 1/4, 1/2,
                                                   1/2, 3/4, 3/4, 1, 1, 1},
                                   SplineWeights -> {1, 1/Sqrt[2], 1, 1/Sqrt[2], 1,
                                                     1/Sqrt[2], 1, 1/Sqrt[2], 1}]]}},
                       Boxed -> False]]

two intersecting spheres


Here's a variation using CirclePoints[] to evenly sample the circle:

BlockRandom[SeedRandom["spherical"];
            c1 = RandomReal[1, 3]; c2 = RandomReal[1, 3];
            r2 = RandomReal[3/2]; r1 = RandomReal[3/2];
            d = EuclideanDistance[c1, c2];
            u = (d^2 + r1^2 - r2^2)/(2 d^2); cc = {1 - u, u}.{c1, c2};
            rc = Sqrt[r1^2 - d^2 u^2];
            n = 30;
            Graphics3D[{Opacity[1/2], Sphere[c1, r1], Sphere[c2, r2],
                        {Red, Sphere[AffineTransform[{RotationMatrix[{{0, 0, 1},
                                                      Normalize[c2 - c1]}], cc}][
                                     PadRight[ArrayPad[CirclePoints[rc, n],
                                                       {{0, 1}, {0, 0}}, "Periodic"],
                                              {Automatic, 3}]], 1/100]}}, Boxed -> False]]

crown of spheres

$\endgroup$
  • $\begingroup$ This is the corresponding picture for the OP's case. $\endgroup$ – J. M. is away Jul 20 '16 at 18:10
  • $\begingroup$ +1 Showing off your Spline skills I see $\endgroup$ – Jason B. Jul 20 '16 at 21:47
  • $\begingroup$ @Jason, it was just a means to an end. ;) Let me edit this to include a variation... $\endgroup$ – J. M. is away Jul 21 '16 at 0:52
2
$\begingroup$

If you want to simply sample the intersection, you can do this:

RandomPoint[
   ImplicitRegion[
    RegionMember[
     RegionIntersection[Sphere[{0, 0, 0}, 1], 
      Sphere[{1, 1, 1}, 1.5]], {x, y, z}], {x, y, z}], 100] // 
  Point // Graphics3D

enter image description here

This should work, but maybe it doesn't because of the single-dimensional nature of the sampling region:

RandomPoint[
     RegionIntersection[Sphere[{0, 0, 0}, 1], 
      Sphere[{1, 1, 1}, 3/2]], 100] // Point // Graphics3D
$\endgroup$
  • $\begingroup$ Also, you can construct actual solutions for equations involved around RegionMember, but that needs more specific information of what you're really after. $\endgroup$ – kirma Jul 20 '16 at 15:15
1
$\begingroup$

Why not just generate the parametric equations for the intersection of the spheres? You know that the intersection points form a circle in a plane perpendicular to the vector between the centers of the spheres.

sphereIntersection[{x1_, r1_}, {x2_, r2_}, th_] :=

  Module[{dvec, d, d1, r, rr, ss, nx, ny, nz},
   dvec = x2 - x1;
   d = Norm@dvec;
   dvec = {nx, ny, nz} = Normalize@dvec;
   d1 = (d^2 + r1^2 - r2^2)/(2 d);
   If[r1^2 < d1^2, Return[$Failed]];
   r = Sqrt[r1^2 - d1^2];
   rr = {1 - nx^2/(1 + nz), -nx ny/(1 + nz), -nx};
   ss = {-nx ny/(1 + nz), 1 - ny^2/(1 + nz), -ny};
   x1 + d1 dvec + r (Cos[th] rr + Sin[th] ss)];
sphereIntersection[s1_Sphere, s2_Sphere, th_] := 
 sphereIntersection[List @@ s1, List @@ s2, th]

You get your equation via

sphereIntersection[sp1, sp2, theta] // FullSimplify
(* {(
 7 (3 + Sqrt[3]) + Sqrt[143] (3 + 2 Sqrt[3]) Cos[theta] - 
  Sqrt[429] Sin[theta])/(24 (3 + Sqrt[3])), (
 7 (3 + Sqrt[3]) - Sqrt[429] Cos[theta] + 
  Sqrt[143] (3 + 2 Sqrt[3]) Sin[theta])/(24 (3 + Sqrt[3])), 
 1/24 (7 - Sqrt[143] Cos[theta] - Sqrt[143] Sin[theta])} *)

And generate points via

sphereIntersection[sp1, sp2, #] & /@ Range[0, 2 π, π/2] // N
(* {{0.972304, 0.109291, -0.206594}, {0.109291, 
  0.972304, -0.206594}, {-0.38897, 0.474043, 
  0.789928}, {0.474043, -0.38897, 0.789928}, {0.972304, 
  0.109291, -0.206594}} *)

And plot it via

Show[Graphics3D[{Opacity[0.2], sp1, sp2}],
 ParametricPlot3D[
  sphereIntersection[sp1, sp2, theta], {theta, 0, 2 π}, 
  Evaluated -> True]
 ]

Mathematica graphics

$\endgroup$
  • $\begingroup$ You might want to see this paper. $\endgroup$ – J. M. is away Jul 21 '16 at 8:56
  • $\begingroup$ @J.M. Thanks, it's better without using a random number to generate the basis vectors $\endgroup$ – Jason B. Jul 21 '16 at 20:49
  • $\begingroup$ I know; I was actually trying to look for this paper when I saw your solution. Apart from that little typo I've now fixed, it's pretty nice! $\endgroup$ – J. M. is away Jul 22 '16 at 1:30

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.