Intersection between a line and a sphere

Question

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 ?

Cowardly

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. – Viki

Answer 1

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.

Fleer

Answer 2

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 };
    }

Etom

The formula for B in my answer was wrong. I have fixed it now. – Fleer

This worked perfectly for me, I was able to rapidly convert it to Java too. – Wadley

It would be helpful if you point out that that code is tied to GPLv3. – Sorel

Answer 3

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 };
        }
    }

Saxena

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). – Saxena

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! – Beauvais

Answer 4

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)

Entity

Answer 5

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.

Polynices

Answer 6

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.

Proverbial

Answer 7

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.

Ashla

Answer 8

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 };

Cussedness

Answer 9

User Ashavsky almost got it right, however his code doesn't handle case when one of the segment points is inside, and the other outside of the sphere. Here is code that handles also those cases (implemented in Unreal Engine 5):

int UGGMathFunctionLibrary::SegmentSphereIntersection_Internal2(FVector linePoint0, FVector linePoint1, FVector sphereCenter, double sphereRadius, FVector* OutIntersections)
{
    /* my handling of a case when one is inside of the sphere and second not. original algorithm doesn't handle this*/ 
    OutIntersections[0] = FVector::ZeroVector;
    OutIntersections[1] = FVector::ZeroVector;

    float Len0 = (sphereCenter - linePoint0).Length();
    float Len1 = (sphereCenter - linePoint1).Length();
    
    if (Len0 < sphereRadius && Len1 < sphereRadius)
    {
        return 0;
    }
    if (Len0 < sphereRadius)
    {
        int NumIntersections = RaySphereIntersection(UGGCoreFunctionLibrary::GetGameWorld(), linePoint0, (linePoint1 - linePoint0).GetSafeNormal(),
            sphereCenter, sphereRadius, OutIntersections[0], OutIntersections[1]);

        check(NumIntersections == 1);
        return 1;
    }
    if (Len1 < sphereRadius)
    {
        int NumIntersections = RaySphereIntersection(UGGCoreFunctionLibrary::GetGameWorld(), linePoint1, (linePoint0 - linePoint1).GetSafeNormal(),
    sphereCenter, sphereRadius, OutIntersections[0], OutIntersections[1]);

        check(NumIntersections == 1);
        return 1;
    }
    // end of custom code
        
    // taken from https://mcmap.net/q/1039536/-intersection-between-a-line-and-a-sphere (user Ashavsky answer)
    double cx = sphereCenter.X;
    double cy = sphereCenter.Y;
    double cz = sphereCenter.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 - sphereRadius * sphereRadius;

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

    double t1 = (-B - FMath::Sqrt(D)) / (2.0 * A);

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

    double t2 = (-B + FMath::Sqrt(D)) / (2.0 * A);
    FVector solution2 = FVector(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 || t1 < 0 || t2 >1 || t2 < 0)
    {
        return 0;
    }
    if (D == 0)
    {
        OutIntersections[0] = solution1;
        return 1;
    }
    
    OutIntersections[0] = solution1;
    OutIntersections[1] = solution2;
    return 2;
}

int UGGMathFunctionLibrary::RaySphereIntersection(UObject* WorldContext, const FVector& RayPos, const FVector& RayDir, const FVector& SpherePos, float SphereRadius, FVector& OutCloserIntersection, FVector& OutFurtherIntersection, float& T1, float& T2, float& Discr, float& b, float& c, float DebugTime /*= -1.0f*/)
{
    OutCloserIntersection = OutFurtherIntersection = FVector::ZeroVector;

    FVector RayDirection = RayDir.GetUnsafeNormal();


    FVector FromSphereToRay = RayPos - SpherePos;
    b = FVector::DotProduct(FromSphereToRay, RayDirection);
    c = FVector::DotProduct(FromSphereToRay, FromSphereToRay) - SphereRadius * SphereRadius; // dot product btwn same vector is length squared

    // No intersection if ray origin is outside of sphere (c > 0.0f) and  ray is pointing away from sphere (b > 0)
    if (c > 0.0f && b > 0.0f)
    {
        return 0;
    }


    Discr = b * b - c;
    // A negative discriminant corresponds to ray missing sphere
    if (Discr < 0.0f)
    {
        return 0;
    }


    T1 = -b - FMath::Sqrt(Discr);
    T2 = -b + FMath::Sqrt(Discr);

    // Ray started inside sphere - one intersection
    if (T1 < 0.0f)
    {
        OutCloserIntersection = RayPos + T2 * RayDirection;
        
        return 1;
    }

    // Two intersections
    OutCloserIntersection = RayPos + T1 * RayDirection;
    OutFurtherIntersection = RayPos + T2 * RayDirection;


    return 2;
}

Epilimnion