intersectLineCircle

INTERSECTLINECIRCLE Intersection point(s) of a line and a circle.

   INTERS = intersectLineCircle(LINE, CIRCLE);
   Returns a 2-by-2-by-N array, containing on each row the coordinates of
   an intersection point for each line-circle pair, i.e. INTERS(:,:,k)
   contains the intersections between LINE(k,:) and CIRCLE(k,:).

   If a line-circle pair does not intersect, the corresponding results are
   set to NaN. 

   Example
     % base point
     center = [10 0];
     % create vertical line
     l1 = [center 0 1];
     % circle
     c1 = [center 5];
     pts = intersectLineCircle(l1, c1)
     pts =
     10   -5
     10    5
     % draw the result
     figure; clf; hold on;
     axis([0 20 -10 10]);
     drawLine(l1);
     drawCircle(c1);
     drawPoint(pts, 'rx');
     axis equal;

   See also
   lines2d, circles2d, intersectLines, intersectCircles

   References
   http://local.wasp.uwa.edu.au/~pbourke/geometry/sphereline/
   http://mathworld.wolfram.com/Circle-LineIntersection.html

0001 function points = intersectLineCircle(line, circle)
0002 %INTERSECTLINECIRCLE Intersection point(s) of a line and a circle.
0003 %
0004 %   INTERS = intersectLineCircle(LINE, CIRCLE);
0005 %   Returns a 2-by-2-by-N array, containing on each row the coordinates of
0006 %   an intersection point for each line-circle pair, i.e. INTERS(:,:,k)
0007 %   contains the intersections between LINE(k,:) and CIRCLE(k,:).
0008 %
0009 %   If a line-circle pair does not intersect, the corresponding results are
0010 %   set to NaN.
0011 %
0012 %   Example
0013 %     % base point
0014 %     center = [10 0];
0015 %     % create vertical line
0016 %     l1 = [center 0 1];
0017 %     % circle
0018 %     c1 = [center 5];
0019 %     pts = intersectLineCircle(l1, c1)
0020 %     pts =
0021 %     10   -5
0022 %     10    5
0023 %     % draw the result
0024 %     figure; clf; hold on;
0025 %     axis([0 20 -10 10]);
0026 %     drawLine(l1);
0027 %     drawCircle(c1);
0028 %     drawPoint(pts, 'rx');
0029 %     axis equal;
0030 %
0031 %   See also
0032 %   lines2d, circles2d, intersectLines, intersectCircles
0033 %
0034 %   References
0035 %   http://local.wasp.uwa.edu.au/~pbourke/geometry/sphereline/
0036 %   http://mathworld.wolfram.com/Circle-LineIntersection.html
0037 %
0038 
0039 % ------
0040 % Authors: David Legland, JuanPi Carbajal
0041 % E-mail: david.legland@inrae.fr, ajuanpi+dev@gmail.com
0042 % Created: 2011-01-14, using Matlab 7.9.0.529 (R2009b)
0043 % Copyright 2011-2024 INRA - Cepia Software Platform
0044 
0045             
0046 % check size of inputs
0047 nLines = size(line, 1);
0048 nCircles = size(circle, 1);
0049 if nLines ~= nCircles
0050   error ('matGeom:geom3d:invalidArguments', ...
0051       'Requires same number of lines and circles');
0052 end
0053             
0054 % center parameters
0055 center = circle(:, 1:2);
0056 radius = circle(:, 3);
0057 
0058 % line parameters
0059 dp = line(:, 1:2) - center;
0060 vl = line(:, 3:4);
0061 
0062 % coefficients of second order equation
0063 a = sum(line(:, 3:4).^2, 2);
0064 b = 2 * sum(dp .* vl, 2);
0065 c = sum(dp.^2, 2) - radius.^2;
0066 
0067 % discriminant
0068 delta = b .^ 2 - 4 * a .* c;
0069 
0070 points = nan(2, 2, nCircles);
0071 
0072 valid = delta >= 0;
0073 
0074 if any(valid)
0075     nValids = sum(valid);
0076     % compute roots (as a N-by-N-by-2 array)
0077     u = bsxfun(@plus, -b(valid), bsxfun(@times, [-1 1], sqrt(delta(valid))));
0078     u = bsxfun(@rdivide, u, a(valid)) / 2;
0079 
0080     if nValids == 1
0081         points = [...
0082             line(1:2) + u(:,1) .* line(3:4); ...
0083             line(1:2) + u(:,2) .* line(3:4)];
0084     else
0085         tmp = [...
0086             line(valid, 1:2) + u(:,1) .* line(valid, 3:4) ...
0087             line(valid, 1:2) + u(:,2) .* line(valid, 3:4)].';
0088         points(:, :, valid) = permute(reshape(tmp, [2, 2, nValids]), [2 1 3]);
0089     end
0090 end