Intersection of two spheres and a plane

Maybe you could draw a picture? I'm not sure what you mean by the plane lying on the sphere of the first equation which is a sphere centered at the origin. The plane equation given will also pass thru the origin. Thus the plane simply bisects rather than lies on the first sphere equation using the normal (a,b,c).

Did you mean ax + by + cz = 1 instead?

It seems to me like you already have most of this problem solved. If you have the equation of the plane the circle is on, then intersect that with the plane ax + by + cz = 0 to get a line. Then you can intersect the line and the circle to get your intersection points.

This problem also lends itself naturally to a 2D approach. You've already calculated the circle plane, and you have the intersection plane. It's possible to then work entirely in the 2D sub-space defined by the two vectors u = (a, b, c) and v = (x_t, y_t, z_t). In this subspace, the circle appears as a line segment (finite length), and the intersection plane appears as a line (infinite length). All we do is see if the line intersects the line segment. One major simplification you should make is if dot(u,v) < 0, let u = -u to make things easier later on. The radius of the circle is sqrt(r² - r⁴/4), which is half the length of the line segment, which is now perpendicular and centered on the v-vector at a distance of 1 - r²/2 from the origin. The line segment will intersect the line if the angle between (a, b, c) and (x_t, y_t, z_t) is large enough. If we let dot(u,v) be the dot product between these vectors, and let t = sqrt(1 - (1 - dot(u,v))²) * (1 - r²/2) / (1 - dot(u,v)), then there's an intersection if t² ≤ r² - r⁴/4. When the intersection plane and the circle plane are parallel, then dot(u,v) is 1, and 1 - dot(u,v) is 0, which gives a divide by 0 and t approaches ∞, which is clearly not less than some other quantity. In the case that the two planes are perpendicular then dot(u,v) will be 0, 1 - dot(u,v) will be 1, and t will be 0 which is always less than or equal to the right-hand quantity (hence a guaranteed intersection).

From this point we can geometrically construct the intersection points. Take the cross product of the two basis vectors and normalize it, cross(u,v)/|cross(u,v)|, and give it a name like n. Then take the cross product between this vector and the v-vector, cross(n,v)/|cross(n,v)|, and give it a name like s. Here comes some heavy vector math, required to go from the sub-space back to full 3D. The intersection points you're looking for will be:

p_i = (1 - r²/2)*(x_t, y_t, z_t) + t*s +/- sqrt(r² - r⁴/4 - t²)*n

The +/- of the square-root gives you two solutions. If t² = r² - r⁴/4, then there is only one solution, as expected. It seems long and complicated, but it's really only a few cross products and square roots.

Hopefully I didn't make any mistakes either, I tend to lose my train of thought sometimes when dealing with a lot of math [smile]

First, find the intersection between your plane ax+by+cz=0 and the sphere-sphere plane x_tx + y_ty + z_tz = 1 - (r²) / 2

n₁ = [a b c]
n₂ = normalize([x_t y_t z_t])

a = n₁ cross n₂

If a is 0, the planes are parallel, so there is no solution.

Now find a point p₀ that is on both planes by saying x₀ = 0 and solving the plane equations as a system. (You'll have to code a special case using y₀ = 0 or z₀ = 0 for when that doesn't work.)

Now you have an equation for a line the solution must be on:

s = p₀ + t*a

where t is a scalar value.

Now just use a line-sphere intersection with the unit sphere to find the solutions.

s_x = p_0x + a_x*t
s_y = p_0y + a_y*t
s_z = p_0z + a_z*t
x² + y² + z² = 1

Substitute and solve the resulting quadratic.
(p_0x + a_x*t)² + (p_0y + a_y*t)² + (p_0z + a_z*t)² = 1

(a_x²+a_y²+a_z²)t² + 2(p_0xa_x + p_0ya_y + p_0za_z)t + (p_0x + p_0y + p_0z) - 1 = 0

Just use the quadratic formula on this to get your answer.

[edit] Oops, I was ninja'd by an hour. That's what I get for starting a post before dinner then doing laundry afterwards. Also, Zipster's solution doesn't require a corner case, so it's a cleaner solution.

Actually I think I might have made a mistake in the equation for t, and it should be t = dot(u,v) * (1 - r²/2) / sqrt(1 - dot(u,v)²). The idea is that you're not interested in the angle between u and v, but the angle between perp(u) and v. So if the dot product of two normalized vectors is cos(θ), then if you use the perpendicular of one of those vectors instead it's cos(θ +/- PI/2) (+/- depending on which produces an absolute angle of less than PI/2), which gives us sin(θ), or sqrt(1 - dot(u,v)²). I was instead just using 1 - dot(u,v).

Thanks for the help, although I did sort of resolve this on my own. What I did was I converted the ax + by + cz = 0 plane into parametric form, where one of the vectors was (a, b, c) x (1, 0, 0) and the other was that product x(a, b, c). I then used x, y, and z (which were then in terms of two parametric variables) in the equation for the plane of intersection of the two spheres. Then I solved the parametric equation for one of the parameters in terms of the other one. This gave me a linear equation: a vector plus a vector times a scalar. Then, I plugged this into my amazing Line-Sphere Intersector Of Science, to get the place where this line intersected the unit sphere. Now, it works (I think).

As for InetRoadKill's question, I said (a, b, c) lies on the sphere, not "the plane ax + by + cz = 0 lies on the sphere." The expression ax + by + cz is indeed equal to 0, not 1.





If you would like to know the Line-Sphere Intersector Of Science^TM algorithm, it is as follows:

For a sphere (x - x₀)² + (y - y₀)² + (z - z₀)² = r², and a line defined as x = x₀ + tv,

let o be the vector from the origin to center of the circle, minus x₀, and
let s be the unit vector parallel to v.

Then the value of t that satisfies the line equation is s dot o +or- sqrt((s dot o)² + r² - |o|²),

and the value of x can be found by evaluating the line equation with the calculated value of t.