×
INTELLIGENT WORK FORUMS
FOR ENGINEERING PROFESSIONALS

Are you an
Engineering professional?
Join Eng-Tips Forums!
• Talk With Other Members
• Be Notified Of Responses
• Keyword Search
Favorite Forums
• Automated Signatures
• Best Of All, It's Free!

*Eng-Tips's functionality depends on members receiving e-mail. By joining you are opting in to receive e-mail.

#### Posting Guidelines

Promoting, selling, recruiting, coursework and thesis posting is forbidden.

# Intersection of a line and the surface of a sphere2

## Intersection of a line and the surface of a sphere

(OP)
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

SIG:Please see FAQ731-376: Eng-Tips.com Forum Policies for tips on how to make the best use of Eng-Tips.

### RE: Intersection of a line and the surface of a sphere

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.

### RE: Intersection of a line and the surface of a sphere

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.

### RE: Intersection of a line and the surface of a sphere

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.

### RE: Intersection of a line and the surface of a sphere

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.

### RE: Intersection of a line and the surface of a sphere

Error alert. Let me try again.

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

### RE: Intersection of a line and the surface of a sphere

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

Cheers

Greg Locock

SIG:Please see FAQ731-376: Eng-Tips.com Forum Policies for tips on how to make the best use of Eng-Tips.

### RE: Intersection of a line and the surface of a sphere

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

M

--
Dr Michael F Platten

### RE: Intersection of a line and the surface of a sphere

(OP)
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

SIG:Please see FAQ731-376: Eng-Tips.com Forum Policies for tips on how to make the best use of Eng-Tips.

### RE: Intersection of a line and the surface of a sphere

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

#### CODE

% 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: http://en.wikipedia.org/wiki/Line%E2%80%93sphere_intersection
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

### RE: Intersection of a line and the surface of a sphere

delete the line that begins "norms"

M

--
Dr Michael F Platten

### RE: Intersection of a line and the surface of a sphere

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

#### CODE

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

### RE: Intersection of a line and the surface of a sphere

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

Cheers

Greg Locock

SIG:Please see FAQ731-376: Eng-Tips.com Forum Policies for tips on how to make the best use of Eng-Tips.

#### Red Flag This Post

Please let us know here why this post is inappropriate. Reasons such as off-topic, duplicates, flames, illegal, vulgar, or students posting their homework.

#### Red Flag Submitted

Thank you for helping keep Eng-Tips Forums free from inappropriate posts.
The Eng-Tips staff will check this out and take appropriate action.

Close Box

# Join Eng-Tips® Today!

Join your peers on the Internet's largest technical engineering professional community.
It's easy to join and it's free.

Here's Why Members Love Eng-Tips Forums:

• Talk To Other Members
• Notification Of Responses To Questions
• Favorite Forums One Click Access
• Keyword Search Of All Posts, And More...

Register now while it's still free!