Thursday, June 25, 2015

java - Sphere-Sphere intersection and Circle-Sphere intersection


I have code for circle-circle intersection. But I need to expand it to 3-D. How do I calculate:



  • Radius and center of the intersection circle of two spheres

  • Points of the intersection of a sphere and a circle?


Given two spheres (sc0,sr0) and (sc1,sr1), I need to calculate a circle of intersection whose center is ci and whose radius is ri.

Moreover, given a sphere (sc0,sr0) and a circle (cc0, cr0), I need to calulate the two intersection points (pi0, pi1)
I have checked this link and this link, but I could not understand the logic behind them and how to code them.


I tried ProGAL library for sphere-sphere-sphere intersection, but the resulting coordinates are rounded. I need precise results.



Answer



Sphere-Sphere Intersection


Let's start with the more obvious one - sphere-sphere. It's almost identical to the circle-circle case in 2D. We can project down on any plane containing the line between the spheres' centers to get an identical 2D picture:


Reducing the problem to 2D


Here the first sphere has center c_1 and radius r_1, the second c_2 and r_2, and their intersection has center c_i and radius r_i. Let d = ||c_2 - c_1||, the distance between the spheres.


Edge Cases:


As was pointed out previously, if r_1 + r_2 < d then there is no intersection at all. The spheres have a separation between them.



If r_1 + r_2 == d, then the intersection is a single point, located a distance of r_1 on the line from c_1 to c_2, or: c_i = c_1 + (c_2 - c_1) * r_1/d


This can also happen from the other side, if d + min(r_1, r_2) = max(r_1, r_2). In this case, the single point of intersection lies on the larger sphere, at c_i = c_larger + (c_smaller - c_larger) * r_larger/d


Finally, if one circle is entirely encompassed within the other, such that d + min(r_1, r_2) < max(r_1, r_2), then once again there is no intersection.


Edge case examples


Typical intersections


Okay, now the meaty bit.


Let's define c_i = c_1 + h * (c_2 - c_1), for some value of h. Note that h can be negative or greater than one in the case where one circle encompasses the center of the other. Fortunately, the signs will "just work," so we don't need to add any special case logic for this condition.


The remaining, interesting cases


We can express r_1 and r_2 in terms of r_i, h, and d using Pythagorean Theorem, then solve the system of equations to find h:


$$\require{cancel} \begin{array}{r|l} (hd)^2 + r_i^2 = r_1^2 & \left((1-h)d\right)^2 + r_i^2 = r_2^2\\ h^2d^2 - r_1^2 = -r_i^2 & \left(1 - 2h + h^2 \right)d^2 - r_2^2 = -r_i^2 \end{array}\\ \begin{align} \cancel{h^2d^2} - r_1^2 &= d^2 - 2hd^2 + \cancel{h^2d^2} - r_2^2\\ 2hd^2 &= d^2 + r_1^2 - r_2^2\\ h &= \frac 1 2 + \frac {r_1^2 - r_2^2} {2 d^2} \end{align}$$



h = 1/2 + (r_1 * r_1 - r_2 * r_2)/(2 * d*d)

We can sub this into our formula for c_i above to find the center of the circle of intersections. Then, reversing one of our earlier Pythagorean relations:


r_i = sqrt(r_1*r_1 - h*h*d*d)

So, now we have the center and radius of our intersection. Now we can revolve this around the separating axis to get our full circle of solutions. The circle lies in a plane perpendicular to the separating axis, so we can take n_i = (c_2 - c_1)/d as the normal of this plane.


Our glorious intersection


After choosing a tangent and bitangent t_i and b_i perpendicular to this normal and each other, you can write any point on this circle as:


p_i(theta) = c_i + r_i * (t_i * cos(theta) + b_i sin(theta));


Because of the Hairy Ball Theorem, there's no one universal way to choose the tangent/bitangent to use. My recommendation would be to pick one of the coordinate axes not parallel to n_i, and set t_i = normalize(cross(axis, n_i)), and b_i = cross(t_i, n_i) or somesuch.




Sphere-Circle Intersection


A circle, being flat/2-dimensional, can only intersect with things within its own plane. So, the first thing is: find out where that plane intersects our sphere.


Slice


Given a circle with center & radius c_c and r_c, respectively, in the plane with (unit) normal n, and a sphere with center & radius c_s and r_s...


The plane of the circle cuts the sphere d = dot(n, c_c - c_s) units from the sphere's center. Again, this can be negative if the normal is pointing toward the sphere's center, but the signs will work out.


If abs(d) > r_s then there is no intersection. Our plane passes above/below the sphere entirely.


If there is an intersection, it will be between our original circle and a new one formed where this plane meets the sphere, with center c_p = c_s + d*n.


If d == r_s then this is the sole point of intersection with the plane; otherwise we have a circle whose radius r_p we can find with Pythagorean Theorem:



r_p = sqrt(r_s*r_s - d*d)

And now we've reduced the problem to circle-circle, which we already know how to solve. Checking the distance/radii involved according to the edge cases above will let us take care of the non-intersecting and single point of intersection cases.


Looks familiar


If the two circles do overlap by more than one point, we can use the formulae from the Sphere-Sphere check to find c_i and r_i as before.


Instead of a full circle of intersection points though, we'll have only the two points lying in the plane of the circle.


To express these points, we can define a tangent vector in the plane,


t = normalize(cross(c_p - c_c, n))

Then our intersection points are:



p_0 = c_i - t * r_i

p_1 = c_i + t * r_i

And we're done. :D


I hope that helps!


No comments:

Post a Comment

Simple past, Present perfect Past perfect

Can you tell me which form of the following sentences is the correct one please? Imagine two friends discussing the gym... I was in a good s...