Finding point of intersection between a line and a sphere

Hello all,

I have a MATLAB code that plots a 3D sphere using:

[xs,ys,zs] = sphere(10);

surface = surf(350*zs+1000,350*ys,350*xs);

I also have a line that represents the normal between three points (labeled P0, P1, P2) on a plane, which is plotted from the middle-point between all three points:

P0 = [tz1,ty1,tx1]; P1 = [tz2,ty2,tx2]; P2 = [tz3,ty3,tx3]; %represent a triangle

Pm = mean([P0;P1;P2]); %represents the midpoint between P0, P1 and P2

normal = cross(P0-P1,P0-P2);

cn = normal + Pm;

normal_vector = plot3([Pm(1) cn(1)],[Pm(2) cn(2)],[Pm(3) cn(3)],'k--'); %normal

What I am trying to do is find the coordinates of the point of intersection between the line "normal_vector" and the sphere "surface".

This is what the plot looks like:

The points P0, P1 and P2 are shown as coloured circles and are always inside the sphere, so their normal is always showing 'outwards' through the surface of the sphere.

Thank you in advance!

Answer 1

[Torsten]

Sphere:

(x-xs)^2 + (y-ys)^2 + (z-zs)^2 = R^2

Line:

[x y z] = Pm + l*normal

Thus

(Pm(1)+l*normal(1)-xs)^2 + (Pm(2)+l*normal(2)-ys)^2 + (Pm(3)+l*normal(3)-zs)^2 = R^2

Solve for (the positive) l.

[André C. D.]

Thank you for the response. I tried that by solving for "l" using the solve function, but unfortunately, I get an 11x11 complex double as a result (when I use the values I have for x, y, z, etc.), which can't be plotted.

Is there something I am missing?

[Torsten]

Executable code available ?

[André C. D.]

The important part of the plotting code looks like this:

[xs,ys,zs] = sphere(10);
   
R = 350; %sphere radius
surface = surf(R*zs+1000,R*ys,R*xs);

hold on

surface.EdgeColor = 'w'; surface.FaceColor = 'k'; surface.FaceAlpha = 0.05; surface.FaceLighting = 'none';

P0 = [tz1,ty1,tx1]; P1 = [tz2,ty2,tx2]; P2 = [tz3,ty3,tx3];
Pm = mean([P0;P1;P2]);

normal = cross(P0-P1,P0-P2);

cn = normal + Pm;

normal1 = normal(1); normal2 = normal(2); normal3 = normal(3); %not actually necessary
Pm1 = Pm(1); Pm2 = Pm(2); Pm3 = Pm(3);
    
normal_vector = plot3([Pm(1) cn(1)],[Pm(2) cn(2)],[Pm(3) cn(3)],'k--');
                                                    
plot3(tz1,ty1,tx1,"o") %tz, ty and tx: coordinates of P0, P1, P2 (1x1 doubles)
plot3(tz2,ty2,tx2,"o")
plot3(tz3,ty3,tx3,"o")

The way I solved for "l" using the equation you gave me:

sol = solve((Pm1+(l*normal1)+xs)^2 + (Pm2+(l*normal2)+ys)^2 + (Pm3+(l*normal3)+zs)^2 == R^2,l)
 
sol = 

    -(Pm1*normal1 + Pm2*normal2 + Pm3*normal3 + normal1*xs + normal2*ys + normal3*zs - (- Pm1^2*normal2^2 - Pm1^2*normal3^2 + 2*Pm1*Pm2*normal1*normal2 + 2*Pm1*Pm3*normal1*normal3 + 2*Pm1*normal1*normal2*ys + 2*Pm1*normal1*normal3*zs - 2*Pm1*normal2^2*xs - 2*Pm1*normal3^2*xs - Pm2^2*normal1^2 - Pm2^2*normal3^2 + 2*Pm2*Pm3*normal2*normal3 - 2*Pm2*normal1^2*ys + 2*Pm2*normal1*normal2*xs + 2*Pm2*normal2*normal3*zs - 2*Pm2*normal3^2*ys - Pm3^2*normal1^2 - Pm3^2*normal2^2 - 2*Pm3*normal1^2*zs + 2*Pm3*normal1*normal3*xs - 2*Pm3*normal2^2*zs + 2*Pm3*normal2*normal3*ys + R^2*normal1^2 + R^2*normal2^2 + R^2*normal3^2 - normal1^2*ys^2 - normal1^2*zs^2 + 2*normal1*normal2*xs*ys + 2*normal1*normal3*xs*zs - normal2^2*xs^2 - normal2^2*zs^2 + 2*normal2*normal3*ys*zs - normal3^2*xs^2 - normal3^2*ys^2)^(1/2))/(normal1^2 + normal2^2 + normal3^2)
    -(Pm1*normal1 + Pm2*normal2 + Pm3*normal3 + normal1*xs + normal2*ys + normal3*zs + (- Pm1^2*normal2^2 - Pm1^2*normal3^2 + 2*Pm1*Pm2*normal1*normal2 + 2*Pm1*Pm3*normal1*normal3 + 2*Pm1*normal1*normal2*ys + 2*Pm1*normal1*normal3*zs - 2*Pm1*normal2^2*xs - 2*Pm1*normal3^2*xs - Pm2^2*normal1^2 - Pm2^2*normal3^2 + 2*Pm2*Pm3*normal2*normal3 - 2*Pm2*normal1^2*ys + 2*Pm2*normal1*normal2*xs + 2*Pm2*normal2*normal3*zs - 2*Pm2*normal3^2*ys - Pm3^2*normal1^2 - Pm3^2*normal2^2 - 2*Pm3*normal1^2*zs + 2*Pm3*normal1*normal3*xs - 2*Pm3*normal2^2*zs + 2*Pm3*normal2*normal3*ys + R^2*normal1^2 + R^2*normal2^2 + R^2*normal3^2 - normal1^2*ys^2 - normal1^2*zs^2 + 2*normal1*normal2*xs*ys + 2*normal1*normal3*xs*zs - normal2^2*xs^2 - normal2^2*zs^2 + 2*normal2*normal3*ys*zs - normal3^2*xs^2 - normal3^2*ys^2)^(1/2))/(normal1^2 + normal2^2 + normal3^2)

And the two solutions I got are 11x11 complex doubles.

[Torsten]

sol = solve((Pm1+(l*normal1) - xs)^2 + (Pm2+(l*normal2) - ys)^2 + (Pm3+(l*normal3) - zs)^2 == R^2,l)

And if Pm1, Pm2, Pm3, normal1, normal2, normal3, xs, ys, zs are scalars, I don't know why you get a 11x11 matrix result.

[André C. D.]

xs, ys and zs are not scalars. The sphere(10) function returns three 11x11 doubles in order to 'simulate' the surface of the sphere.

That is probably where the issue is at; however, I am not sure how to get from the xs, ys, zs matrices to ordinary coordinate vectors that represent every point in the surface.

[Torsten]

Ah, I thought (xs,ys,zs) is the center of the sphere.

You have to solve

sol = solve((Pm1+(l*normal1) - xc)^2 + (Pm2+(l*normal2) - yc)^2 + (Pm3+(l*normal3) - zc)^2 == R^2,l)

where (xc,yc,zc) is the center of the sphere.

[André C. D.]

Ahh thank you so much, now it works perfectly!

Answer 2

[EVELYN ROSSANA PARRA LOPEZ]

The last line of code is summarized in replacing the terms x, y and z of the parametric equation of a line in space, in the equation that describes a sphere, and the variable to be found is the parameter, in this case l. I apply the same with a sphere and a known line, but the answer is as follows:

CODE LINES:

syms t

sol=solve((0.2118+t*0.8473-1).^2+(0.06883+t*0.2754-0.5).^2+(0.1135+t*0.4541-0.5).^2==((0.25).^2),t)

RESULT:

sol=

240523932/249992315 - (10^(1/2)*3160661400392057^(1/2))/999969260

(10^(1/2)*3160661400392057^(1/2))/999969260 + 240523932/249992315

Answer 3

[Jason Der]

The accepted answer works, but I wasn't able to find a way to vectorize it for a large dataset.

However, I did read an elegant solution on eng-tips at the following url:

https://www.eng-tips.com/viewthread.cfm?qid=248900

Hope this helps someone one day.