intersectLineSphere

INTERSECTLINESPHERE Return intersection points between a line and a sphere.

   PTS = intersectLineSphere(LINE, SPHERE);
   Returns the two points which are the intersection of the given line and
   sphere. 
   LINE   : [x0 y0 z0  dx dy dz]
   SPHERE : [xc yc zc  R]
   PTS     : [x1 y1 z1 ; x2 y2 z2]
   If there is no intersection between the line and the sphere, return a
   2-by-3 array containing only NaN.

   Example
     % draw the intersection between a sphere and a collection of parallel
     % lines 
     sphere = [50.12 50.23 50.34 40];
     [x, y] = meshgrid(10:10:90, 10:10:90);
     n = numel(x);
     lines = [x(:) y(:) zeros(n,1) zeros(n,2) ones(n,1)];
     figure; hold on; axis equal;
     axis([0 100 0 100 0 100]); view(3);
     drawSphere(sphere);
     drawLine3d(lines);
     pts = intersectLineSphere(lines, sphere);
     drawPoint3d(pts, 'rx');

     % apply rotation on set of lines to check with non vertical lines
     rot = eulerAnglesToRotation3d(20, 30, 10);
     rot2 = recenterTransform3d(rot, [50 50 50]);
     lines2 = transformLine3d(lines, rot2);
     figure; hold on; axis equal;
     axis([0 100 0 100 0 100]); view(3);
     drawSphere(sphere);
     drawLine3d(lines2);
     pts2 = intersectLineSphere(lines2, sphere);
     drawPoint3d(pts, 'rx');

   See also
   spheres, circles3d, intersectPlaneSphere

0001 function points = intersectLineSphere(line, sphere, varargin)
0002 %INTERSECTLINESPHERE Return intersection points between a line and a sphere.
0003 %
0004 %   PTS = intersectLineSphere(LINE, SPHERE);
0005 %   Returns the two points which are the intersection of the given line and
0006 %   sphere.
0007 %   LINE   : [x0 y0 z0  dx dy dz]
0008 %   SPHERE : [xc yc zc  R]
0009 %   PTS     : [x1 y1 z1 ; x2 y2 z2]
0010 %   If there is no intersection between the line and the sphere, return a
0011 %   2-by-3 array containing only NaN.
0012 %
0013 %   Example
0014 %     % draw the intersection between a sphere and a collection of parallel
0015 %     % lines
0016 %     sphere = [50.12 50.23 50.34 40];
0017 %     [x, y] = meshgrid(10:10:90, 10:10:90);
0018 %     n = numel(x);
0019 %     lines = [x(:) y(:) zeros(n,1) zeros(n,2) ones(n,1)];
0020 %     figure; hold on; axis equal;
0021 %     axis([0 100 0 100 0 100]); view(3);
0022 %     drawSphere(sphere);
0023 %     drawLine3d(lines);
0024 %     pts = intersectLineSphere(lines, sphere);
0025 %     drawPoint3d(pts, 'rx');
0026 %
0027 %     % apply rotation on set of lines to check with non vertical lines
0028 %     rot = eulerAnglesToRotation3d(20, 30, 10);
0029 %     rot2 = recenterTransform3d(rot, [50 50 50]);
0030 %     lines2 = transformLine3d(lines, rot2);
0031 %     figure; hold on; axis equal;
0032 %     axis([0 100 0 100 0 100]); view(3);
0033 %     drawSphere(sphere);
0034 %     drawLine3d(lines2);
0035 %     pts2 = intersectLineSphere(lines2, sphere);
0036 %     drawPoint3d(pts, 'rx');
0037 %
0038 %   See also
0039 %   spheres, circles3d, intersectPlaneSphere
0040 %
0041 
0042 %   ---------
0043 %   author : David Legland
0044 %   INRA - TPV URPOI - BIA IMASTE
0045 %   created the 18/02/2005.
0046 %
0047 
0048 %   HISTORY
0049 %   2011-06-21 bug for tangent lines, add tolerance
0050 
0051 
0052 %% Process input arguments
0053 
0054 % check if user-defined tolerance is given
0055 tol = 1e-14;
0056 if ~isempty(varargin)
0057     tol = varargin{1};
0058 end
0059 
0060 % difference between centers
0061 dc = bsxfun(@minus, line(:, 1:3), sphere(:, 1:3));
0062 
0063 % equation coefficients
0064 a = sum(line(:, 4:6) .* line(:, 4:6), 2);
0065 b = 2 * sum(bsxfun(@times, dc, line(:, 4:6)), 2);
0066 c = sum(dc.*dc, 2) - sphere(:,4).*sphere(:,4);
0067 
0068 % solve equation
0069 delta = b.*b - 4*a.*c;
0070 
0071 % initialize empty results
0072 points = NaN * ones(2 * size(delta, 1), 3);
0073 
0074 
0075 %% process couples with two intersection points
0076 
0077 % process couples with two intersection points
0078 inds = find(delta > tol);
0079 if ~isempty(inds)
0080     % delta positive: find two roots of second order equation
0081     u1 = (-b(inds) -sqrt(delta(inds))) / 2 ./ a(inds);
0082     u2 = (-b(inds) +sqrt(delta(inds))) / 2 ./ a(inds);
0083     
0084     % convert into 3D coordinate
0085     points(inds, :) = line(inds, 1:3) + bsxfun(@times, u1, line(inds, 4:6));
0086     points(inds+length(delta),:) = line(inds, 1:3) + bsxfun(@times, u2, line(inds, 4:6));
0087 end
0088 
0089 
0090 %% process couples with one intersection point
0091 
0092 % proces couples with two intersection points
0093 inds = find(abs(delta) < tol);
0094 if ~isempty(inds)
0095     % delta around zero: find unique root, and convert to 3D coord.
0096     u = -b(inds) / 2 ./ a(inds);
0097     
0098     % convert into 3D coordinate
0099     pts = line(inds, 1:3) + bsxfun(@times, u, line(inds, 4:6));
0100     points(inds, :) = pts;
0101     points(inds+length(delta),:) = pts;
0102 end