Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

Intersection between a line and a sphere

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 ?

like image 427
Tanguy Avatar asked May 04 '11 12:05

Tanguy


People also ask

How do you find the intersection of a sphere and a line?

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.

Is the line of intersection of a sphere and a plane through its center?

It is a circle. Such a circle formed by the intersection of a plane and a sphere is called the circle of a sphere.

What is the intersection of any plane and a sphere?

When a spherical surface and a plane intersect, the intersection is a point or a circle.

How do you see if a vector intersects a sphere?

You can work out the intersection point of the vectors using linear algebra, and hence the length of the line (or more efficiently the square of the length of the line) and test if this is less than the radius (or the square of the radius) of your sphere. This is the single correct answer on this post.


2 Answers

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.

like image 53
Markus Jarderot Avatar answered Sep 19 '22 17:09

Markus Jarderot


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 };
    }
like image 37
M Katz Avatar answered Sep 23 '22 17:09

M Katz