9

I'm trying to find the point of intersection between a sphere and a line but honestly, I don't have any idea of how to do so. Could anyone help me on this one ?

1
  • Look up "ray sphere intersection" - the same test is used all of the time in ray-tracing and there's plenty of examples online, and even quite a few here on stackoverflow.
    – Alnitak
    May 4, 2011 at 12:24

8 Answers 8

17

Express the line as an function of t:

{ x(t) = x0*(1-t) + t*x1
{ y(t) = y0*(1-t) + t*y1
{ z(t) = z0*(1-t) + t*z1

When t = 0, it will be at one end-point (x0,y0,z0). When t = 1, it will be at the other end-point (x1,y1,z1).

Write a formula for the distance to the center of the sphere (squared) in t (where (xc,yc,zc) is the center of the sphere):

f(t) = (x(t) - xc)^2 + (y(t) - yc)^2 + (z(t) - zc)^2

Solve for t when f(t) equals R^2 (R being the radius of the sphere):

(x(t) - xc)^2 + (y(t) - yc)^2 + (z(t) - zc)^2 = R^2

A = (x0-xc)^2 + (y0-yc)^2 + (z0-zc)^2 - R^2
B = (x1-xc)^2 + (y1-yc)^2 + (z1-zc)^2 - A - C - R^2
C = (x0-x1)^2 + (y0-y1)^2 + (z0-z1)^2

Solve A + B*t + C*t^2 = 0 for t. This is a normal quadratic equation.

You can get up to two solutions. Any solution where t lies between 0 and 1 are valid.

If you got a valid solution for t, plug it in the first equations to get the point of intersection.

I assumed you meant a line segment (two end-points). If you instead want a full line (infinite length), then you could pick two points along the line (not too close), and use them. Also let t be any real value, not just between 0 and 1.

Edit: I fixed the formula for B. I was mixing up the signs. Thanks M Katz, for mentioning that it didn't work.

10

I believe there is an inaccuracy in the solution by Markus Jarderot. Not sure what the problem is, but I'm pretty sure I translated it faithfully to code, and when I tried to find the intersection of a line segment known to cross into a sphere, I got a negative discriminant (no solutions).

I found this: http://www.codeproject.com/Articles/19799/Simple-Ray-Tracing-in-C-Part-II-Triangles-Intersec, which gives a similar but slightly different derivation.

I turned that into the following C# code and it works for me:

    public static Point3D[] FindLineSphereIntersections( Point3D linePoint0, Point3D linePoint1, Point3D circleCenter, double circleRadius )
    {
        // http://www.codeproject.com/Articles/19799/Simple-Ray-Tracing-in-C-Part-II-Triangles-Intersec

        double cx = circleCenter.X;
        double cy = circleCenter.Y;
        double cz = circleCenter.Z;

        double px = linePoint0.X;
        double py = linePoint0.Y;
        double pz = linePoint0.Z;

        double vx = linePoint1.X - px;
        double vy = linePoint1.Y - py;
        double vz = linePoint1.Z - pz;

        double A = vx * vx + vy * vy + vz * vz;
        double B = 2.0 * (px * vx + py * vy + pz * vz - vx * cx - vy * cy - vz * cz);
        double C = px * px - 2 * px * cx + cx * cx + py * py - 2 * py * cy + cy * cy +
                   pz * pz - 2 * pz * cz + cz * cz - circleRadius * circleRadius;

        // discriminant
        double D = B * B - 4 * A * C;

        if ( D < 0 )
        {
            return new Point3D[ 0 ];
        }

        double t1 = ( -B - Math.Sqrt ( D ) ) / ( 2.0 * A );

        Point3D solution1 = new Point3D( linePoint0.X * ( 1 - t1 ) + t1 * linePoint1.X,
                                         linePoint0.Y * ( 1 - t1 ) + t1 * linePoint1.Y,
                                         linePoint0.Z * ( 1 - t1 ) + t1 * linePoint1.Z );
        if ( D == 0 )
        {
            return new Point3D[] { solution1 };
        }

        double t2 = ( -B + Math.Sqrt( D ) ) / ( 2.0 * A );
        Point3D solution2 = new Point3D( linePoint0.X * ( 1 - t2 ) + t2 * linePoint1.X,
                                         linePoint0.Y * ( 1 - t2 ) + t2 * linePoint1.Y,
                                         linePoint0.Z * ( 1 - t2 ) + t2 * linePoint1.Z );

        // prefer a solution that's on the line segment itself

        if ( Math.Abs( t1 - 0.5 ) < Math.Abs( t2 - 0.5 ) )
        {
            return new Point3D[] { solution1, solution2 };
        }

        return new Point3D[] { solution2, solution1 };
    }
3
  • The formula for B in my answer was wrong. I have fixed it now. Jul 6, 2013 at 11:54
  • This worked perfectly for me, I was able to rapidly convert it to Java too. Apr 22, 2015 at 11:06
  • It would be helpful if you point out that that code is tied to GPLv3.
    – Semimono
    Jul 16, 2019 at 22:12
4

Don't have enough reputation to comment on M Katz answer, but his answer assumes that the line can go on infinitely in each direction. If you need only the line SEGMENT's intersection points, you need t1 and t2 to be less than one (based on the definition of a parameterized equation). Please see my answer in C# below:

        public static Point3D[] FindLineSphereIntersections(Point3D linePoint0, Point3D linePoint1, Point3D circleCenter, double circleRadius)
    {

        double cx = circleCenter.X;
        double cy = circleCenter.Y;
        double cz = circleCenter.Z;

        double px = linePoint0.X;
        double py = linePoint0.Y;
        double pz = linePoint0.Z;

        double vx = linePoint1.X - px;
        double vy = linePoint1.Y - py;
        double vz = linePoint1.Z - pz;

        double A = vx * vx + vy * vy + vz * vz;
        double B = 2.0 * (px * vx + py * vy + pz * vz - vx * cx - vy * cy - vz * cz);
        double C = px * px - 2 * px * cx + cx * cx + py * py - 2 * py * cy + cy * cy +
                   pz * pz - 2 * pz * cz + cz * cz - circleRadius * circleRadius;

        // discriminant
        double D = B * B - 4 * A * C;

        double t1 = (-B - Math.Sqrt(D)) / (2.0 * A);

        Point3D solution1 = new Point3D(linePoint0.X * (1 - t1) + t1 * linePoint1.X,
                                         linePoint0.Y * (1 - t1) + t1 * linePoint1.Y,
                                         linePoint0.Z * (1 - t1) + t1 * linePoint1.Z);

        double t2 = (-B + Math.Sqrt(D)) / (2.0 * A);
        Point3D solution2 = new Point3D(linePoint0.X * (1 - t2) + t2 * linePoint1.X,
                                         linePoint0.Y * (1 - t2) + t2 * linePoint1.Y,
                                         linePoint0.Z * (1 - t2) + t2 * linePoint1.Z);

        if (D < 0 || t1 > 1 || t2 >1)
        {
            return new Point3D[0];
        }
        else if (D == 0)
        {
            return new [] { solution1 };
        }
        else
        {
            return new [] { solution1, solution2 };
        }
    }
2
  • One more comment, my code is not complete - definition of a parametric line is that 0<t<1. As such, code should be updated slightly with (D < 0 || t1 > 1 || t2 >1) being changed to (D < 0 || t1 > 1 || t1 < 0 || t2 >1 || t2 < 0).
    – Ashavsky
    Jul 28, 2015 at 14:43
  • thanks for pasting your code. I might have found a bug in it though. I'm looking at a line segment that starts outside a sphere and enters a sphere but the intersection is returning null (empty array). Please see: paste.ofcode.org/aJSjXFpXvvRGf5hYcn8M9Q. I'd much appreciate any comment on this, thank you!
    – Jason
    Feb 23, 2016 at 19:57
2

You may use Wolfram Alpha to solve it in the coordinate system where the sphere is centered.

In this system, the equations are:

Sphere:

     x^2 + y^2 + z^2 = r^2  

Straight line:

    x = x0 + Cos[x1] t
    y = y0 + Cos[y1] t
    z = z0 + Cos[z1] t 

Then we ask Wolfram Alpha to solve for t: (Try it!)

and after that you may change again to your original coordinate system (a simple translation)

1

Find the solution of the two equations in (x,y,z) describing the line and the sphere.

There may be 0, 1 or 2 solutions.

  • 0 implies they don't intersect
  • 1 implies the line is a tangent to the sphere
  • 2 implies the line passes through the sphere.
1

Here's a more concise formulation using inner products, less than 100 LOCs, and no external links. Also, the question was asked for a line, not a line segment.

Assume that the sphere is centered at C with radius r. The line is described by P+l*D where D*D=1. P and C are points, D is a vector, l is a number.

We set PC = P-C, pd = PC*D and s = pd*pd - PC*PC + r*r. If s < 0 there are no solutions, if s == 0 there is just one, otherwise there are two. For the solutions we set l = -pd +- sqrt(s), then plug into P+l*D.

0

Or you can just find the formula of both:
line: (x-x0)/a=(y-y0)/b=(z-z0)/c, which are symmetric equations of the line segment between the points you can find.
sphere: (x-xc)^2+(y-yc)^2+(z-zc)^2 = R^2.

Use the symmetric equation to find relationship between x and y, and x and z.

Then plug in y and z in terms of x into the equation of the sphere.
Then find x, and then you can find y and z.

If x gives you an imaginary result, that means the line and the sphere doesn't intersect.

0

I don't have the reputation to comment on Ashavsky's solution, but the check at the end needed a bit more tweaking.

if (D < 0)
    return new Point3D[0];
else if ((t1 > 1 || t1 < 0) && (t2 > 1 || t2 < 0))
    return new Point3D[0];
else if (!(t1 > 1 || t1 < 0) && (t2 > 1 || t2 < 0))
    return new [] { solution1 };
else if ((t1 > 1 || t1 < 0) && !(t2 > 1 || t2 < 0))
    return new [] { solution2 };
else if (D == 0)
    return new [] { solution1 };
else
    return new [] { solution1, solution2 };

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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