% CIRCLEINTERSECT - Finds intersection of two circles.%% Function to return the intersection points between two circles% given their centres and radi.%% Usage: [i1, i2] = circleintersect(c1, r1, c2, r2, lr)%% Where:%         c1 and c2 are 2-vectors specifying the centres of the two circles.%         r1 and r2 are the radii of the two circles.%         i1 and i2 are the two 2D intersection points (if they exist)%         lr is an optional string specifying what solution is wanted:%            'l' for the solution to the left of the line from c1 to c2%            'r' for the solution to the right of the line from c1 to c2%            'lr' if both solutions are wanted (default).%  Peter Kovesi  %  School of Computer Science & Software Engineering%  The University of Western Australia%  Peter Kovesi%  pk @ csse uwa edu au%  http://www.csse.uwa.edu.au/~pk%%  April   2000   - original version%  October 2003   - mods to allow selection of left/right solutions%                   and catching of degenerate trianglesfunction [i1, i2] = circleintersect(c1, r1, c2, r2, lr)if nargin == 4  lr = 'lr';endmaxmag = max([max(r1, r2), max(c1), max(c2)]);  % maximum value in data inputtol = 100*(maxmag+1)*eps;  % Tolerance used for deciding whether values are equal                           % scaling by 100 is a bit arbitrary...bv = (c2-c1);         % Baseline vector from c1 to c2b = norm(bv);         % The distance between the centres% Trap case of  baseline of zero length.  If r1 == r2% we have a valid geometric situation, but there are an infinite number of% solutions.  Here we simply choose to add r1 in the x direction to c1.if b < eps & abs(r1-r2) < tol    i1 = c1 + [r1 0];  i2 = i1;  returnendbv = bv/b;            % Normalise baseline vector.bvp = [-bv(2) bv(1)]; % Vector perpendicular to baseline% Trap the degenerate cases where one of the radii are zero, or nearly zeroif r1 < tol & abs(b-r2) < tol  i1 = c1;  i2 = c1;  return;elseif r2 < tol & abs(b-r1) < tol  i1 = c2;  i2 = c2;  return;end% Check triangle inequalityif b > (r1+r2) | r1 > (b+r2) | r2 > (b+r1)  c1,  c2, r1, r2, b  error('No solution to circle intersection');end% Normal solutioncosR2 = (b^2 + r1^2 - r2^2)/(2*b*r1);sinR2 = sqrt(1-cosR2^2);if strcmpi(lr,'lr')                      i1 = c1 + r1*cosR2*bv + r1*sinR2*bvp;   % 'left' solution  i2 = c1 + r1*cosR2*bv - r1*sinR2*bvp;   % and 'right solutionelseif strcmpi(lr,'l')  i1 = c1 + r1*cosR2*bv + r1*sinR2*bvp;   % 'left' solution  i2 = [];elseif strcmpi(lr,'r')  i1 = c1 + r1*cosR2*bv - r1*sinR2*bvp;   % 'right solution  i2 = [];else  error('illegal left/right solution request');end