Code covered by the BSD License  

Highlights from
geom2d

image thumbnail
from geom2d by David Legland
Geometry library for matlab. Performs geometric computations on points, lines, circles, polygons...

intersectLineCircle(line, circle) 
function points = intersectLineCircle(line, circle)
%INTERSECTLINECIRCLE Intersection point(s) of a line and a circle
%
%   INTERS = intersectLineCircle(LINE, CIRCLE);
%   Returns a 2-by-2 array, containing on each row the coordinates of an
%   intersection point. If the line and circle do not intersect, the result
%   is filled with 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
%
% ------
% Author: David Legland
% e-mail: david.legland@grignon.inra.fr
% Created: 2011-01-14,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2011 INRA - Cepia Software Platform.

%   HISTORY
%   2011-06-06 fix bug in delta test


% local precision
eps = 1e-14;

% center parameters
center = circle(:, 1:2);
radius = circle(:, 3);

% line parameters
dp = line(:, 1:2) - center;
vl = line(:, 3:4);

% coefficient of second order equation
a = sum(line(:, 3:4).^2, 2);
b = 2*sum(dp .* vl, 2);
c =  sum(dp.^2, 2) - radius.^2;

% discriminant
delta = b .^ 2 - 4 * a .* c;

if delta > eps
    % find two roots of second order equation
    u1 = (-b - sqrt(delta)) / 2 ./ a;
    u2 = (-b + sqrt(delta)) / 2 ./ a;
    
    % convert into 2D coordinate
    points = [line(1:2) + u1 * line(3:4) ; line(1:2) + u2 * line(3:4)];

elseif abs(delta) < eps
    % find unique root, and convert to 2D coord.
    u = -b / 2 ./ a;    
    points = line(1:2) + u*line(3:4);
    
else
    % fill with NaN
    points = NaN * ones(2, 2);
    return;
end

Contact us