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 % 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 % Author: David Legland, david.legland@inrae.fr
0041 % Author: JuanPi Carbajal <ajuanpi+dev@gmail.com>
0042 % Created: 2011-01-14,    using Matlab 7.9.0.529 (R2009b)
0043 % Copyright 2011 INRA - Cepia Software Platform.
0044 
0045 % HISTORY
0046 % 2011-06-06 fix bug in delta test
0047 % 2017-05-05 included some suggestions from code by JuanPi Carbajal <ajuanpi+dev@gmail.com>
0048 % 2017-08-08 update doc
0049             
0050 % check size of inputs
0051 nLines = size(line, 1);
0052 nCircles = size(circle, 1);
0053 if nLines ~= nCircles
0054   error ('matGeom:geom3d:invalidArguments', ...
0055       'Requires same number of lines and circles');
0056 end
0057             
0058 % center parameters
0059 center = circle(:, 1:2);
0060 radius = circle(:, 3);
0061 
0062 % line parameters
0063 dp = line(:, 1:2) - center;
0064 vl = line(:, 3:4);
0065 
0066 % coefficients of second order equation
0067 a = sum(line(:, 3:4).^2, 2);
0068 b = 2 * sum(dp .* vl, 2);
0069 c = sum(dp.^2, 2) - radius.^2;
0070 
0071 % discriminant
0072 delta = b .^ 2 - 4 * a .* c;
0073 
0074 points = nan(2, 2, nCircles);
0075 
0076 valid = delta >= 0;
0077 
0078 if any(valid)
0079     % compute roots (as a N-by-N-by-2 array)
0080     u = bsxfun(@plus, -b(valid), bsxfun(@times, [-1 1], sqrt(delta(valid))));
0081     u = bsxfun(@rdivide, u, a(valid)) / 2;
0082 
0083     if sum(valid) == 1
0084         points = [...
0085             line(1:2) + u(:,1) .* line(3:4); ...
0086             line(1:2) + u(:,2) .* line(3:4)];
0087     else
0088         tmp = [...
0089             line(valid, 1:2) + u(:,1) .* line(valid, 3:4) ...
0090             line(valid, 1:2) + u(:,2) .* line(valid, 3:4)].';
0091         points(:, :, valid) = permute(reshape(tmp, [2, 2, nCircles]), [2 1 3]);
0092     end
0093 end