Intersection of a line and the surface of a sphere 2

[GregLocock]

Hi I need an efficient Matlab algorithm for the intersection of a line in 3d space (or its extension) defined by 2 points A and B , and a sphere of radius r and known centre location C.

In this case the line always intersects the sphere, at two points D and E, ie I don't need error checking for the tangential or missing case, and I am only interested in the solution for whichever of D or E is nearest A.

This calculation needs repeating millions of times, hence a fast solution would be preferred.

Any hints?

Cheers

Greg Locock

SIGlease see FAQ731-376 for tips on how to make the best use of Eng-Tips.

 

[electricpete]

Let C be center of circle and R be radius of circle
Let L1 and L2 be points at beginning and end of the line.
Let t be a variable for parameterization of the line:
Points on the line satisfy: L1 + t*(L2-L1)

The intersections will be on the line and be a distance R from C.

Using the distance equation:
R^2 = [Cx -<L1x + t(L2x-L1x)>]^2 +R^2 + [Cy -<L1y + t(L2y-L1y)>]^2 +[Cz -<L1z + t(L2z-L1z)>]^2

Simplify it and you have a quadratic equation in t. Two solutions, one for each intersection.


=====================================
Eng-tips forums: The best place on the web for engineering discussions.

 

[electricpete]

Whoops - Had an extra R^2 floating around in the distance equation for some reason. Should've been:
R^2 = [Cx -<L1x + t(L2x-L1x)>]^2 [Cy -<L1y + t(L2y-L1y)>]^2 +[Cz -<L1z + t(L2z-L1z)>]^2


=====================================
Eng-tips forums: The best place on the web for engineering discussions.

 

[electricpete]

Try one more time:
R^2 = [Cx -<L1x + t(L2x-L1x)>]^2 + [Cy -<L1y + t(L2y-L1y)>]^2 +[Cz -<L1z + t(L2z-L1z)>]^2

=====================================
Eng-tips forums: The best place on the web for engineering discussions.

 

[electricpete]

I presume with some work you could solve the quadratic equation by hand to give the desired solution t in terms of all the coordinates.

Attached I let the computer do the work. Defining M = L2-L1 simplifies the algebra just a little but still a long formula.

=====================================
Eng-tips forums: The best place on the web for engineering discussions.

 

[electricpete]

Error alert. Let me try again.

=====================================
Eng-tips forums: The best place on the web for engineering discussions.

 

[electricpete]

Attached I have corrected the error in eq1. Even a longer solution ;-(

=====================================
Eng-tips forums: The best place on the web for engineering discussions.

 

[GregLocock]

Thanks, I see where you are going, using paramaterization on the line is clever.



Cheers

Greg Locock

SIGlease see FAQ731-376 for tips on how to make the best use of Eng-Tips.

 

[GregLocock]

Here's my human readable solution



Cheers

Greg Locock

SIGlease see FAQ731-376 for tips on how to make the best use of Eng-Tips.

 

[MikeyP]

Greg, can the repeats of your calculation run concurrently or are they iterative? i.e can you vectorise?

M

--
Dr Michael F Platten

 

[GregLocock]

Turns out that the above solution works fine, in the form shown it runs what I'd call vectorised, ie for all the values of Ax(etc) at once.

Cheers

Greg Locock

SIGlease see FAQ731-376 for tips on how to make the best use of Eng-Tips.

 

[MikeyP]

OK here's my attempt.
It runs 1,000,000 examples in a smidge over 1 second.

% generate a load of lines and spheres that we know will intersect
clear all
n = 1e6; % number of examples to generate

% A and B are 2 points on the line
A = 10.*(rand(3,n)-0.5); % values between -5 and 5
B = 10.*(rand(3,n)-0.5); % values between -5 and 5

% C is the centre of the sphere
C = 2.*(rand(3,n)-0.5); % values between -1 and 1

% r is the radius of the sphere
r = 20.*(rand(1,n));

% CALCULATE INTERSECTION

% profile on

% move the origin to A
B = B-A;
C = C-A;

% calculate unit direction vector of line

norms = sqrt(dot(B,B));
U = B.*repmat(1./sqrt(dot(B,B)),[3 1]);

% find distance along line where intersects occur
% See: [URL unfurl="true"]http://en.wikipedia.org/wiki/Line%E2%80%93sphere_intersection[/URL]
D = dot(U,C);
E = sqrt( dot(U,C).^2 - dot(C,C) + r.^2 );
F = [D+E; D-E];

% find which solution is closest to origin
[dummy, idx] = min(abs(F));
G = F(idx);

% find the intersect point
intersect = U.*repmat(G, [3 1]) + A;

% profile off

--
Dr Michael F Platten

 

[MikeyP]

delete the line that begins "norms"

M

--
Dr Michael F Platten

 

[MikeyP]

Slight improvement: down to 0.67 secs for 1e6 evaluations (assumes A, B, C and r are set up as before)

B = B-A;
C = C-A;
U = B./repmat(sqrt(dot(B,B)),[3 1]);
D = dot(U,C);
E = sqrt( D.^2 - dot(C,C) + r.^2 );
F = [D+E; D-E];
[dummy, idx] = min(abs(F));
G = F(idx);
intersect = U.*repmat(G, [3 1]) + A;

--
Dr Michael F Platten

 

[GregLocock]

Well, thanks for that, very elegant, and it certainly uses matlab in ways I haven't so far.



Cheers

Greg Locock

SIGlease see FAQ731-376 for tips on how to make the best use of Eng-Tips.