intersectCircles

INTERSECTCIRCLES Intersection points of two circles.

   POINTS = intersectCircles(CIRCLE1, CIRCLE2)
   Computes the intersetion point of the two circles CIRCLE1 and CIRCLE1.
   Both circles are given with format: [XC YC R], with (XC,YC) being the
   coordinates of the center and R being the radius.
   POINTS is a 2-by-2 array, containing coordinate of an intersection
   point on each row. 
   In the case of tangent circles, the intersection is returned twice. It
   can be simplified by using the 'unique' function.

   Example
     % intersection points of two distant circles
     c1 = [0  0 10];
     c2 = [10 0 10];
     pts = intersectCircles(c1, c2)
     pts =
         5   -8.6603
         5    8.6603

     % intersection points of two tangent circles
     c1 = [0  0 10];
     c2 = [20 0 10];
     pts = intersectCircles(c1, c2)
     pts =
         10    0
         10    0
     pts2 = unique(pts, 'rows')
     pts2 = 
         10    0

   References
   http://local.wasp.uwa.edu.au/~pbourke/geometry/2circle/
   http://mathworld.wolfram.com/Circle-CircleIntersection.html

   See also
   circles2d, intersectLineCircle, radicalAxis

0001 function points = intersectCircles(circle1, circle2)
0002 %INTERSECTCIRCLES Intersection points of two circles.
0003 %
0004 %   POINTS = intersectCircles(CIRCLE1, CIRCLE2)
0005 %   Computes the intersetion point of the two circles CIRCLE1 and CIRCLE1.
0006 %   Both circles are given with format: [XC YC R], with (XC,YC) being the
0007 %   coordinates of the center and R being the radius.
0008 %   POINTS is a 2-by-2 array, containing coordinate of an intersection
0009 %   point on each row.
0010 %   In the case of tangent circles, the intersection is returned twice. It
0011 %   can be simplified by using the 'unique' function.
0012 %
0013 %   Example
0014 %     % intersection points of two distant circles
0015 %     c1 = [0  0 10];
0016 %     c2 = [10 0 10];
0017 %     pts = intersectCircles(c1, c2)
0018 %     pts =
0019 %         5   -8.6603
0020 %         5    8.6603
0021 %
0022 %     % intersection points of two tangent circles
0023 %     c1 = [0  0 10];
0024 %     c2 = [20 0 10];
0025 %     pts = intersectCircles(c1, c2)
0026 %     pts =
0027 %         10    0
0028 %         10    0
0029 %     pts2 = unique(pts, 'rows')
0030 %     pts2 =
0031 %         10    0
0032 %
0033 %   References
0034 %   http://local.wasp.uwa.edu.au/~pbourke/geometry/2circle/
0035 %   http://mathworld.wolfram.com/Circle-CircleIntersection.html
0036 %
0037 %   See also
0038 %   circles2d, intersectLineCircle, radicalAxis
0039 %
0040 
0041 % ------
0042 % Author: David Legland
0043 % E-mail: david.legland@inrae.fr
0044 % Created: 2011-01-20, using Matlab 7.9.0.529 (R2009b)
0045 % Copyright 2011-2024 INRA - Cepia Software Platform
0046 
0047 % adapt sizes of inputs
0048 n1 = size(circle1, 1);
0049 n2 = size(circle2, 1);
0050 if n1 ~= n2
0051     if n1 > 1 && n2 == 1
0052         circle2 = repmat(circle2, n1, 1);
0053     elseif n2 > 1 && n1 == 1
0054         circle1 = repmat(circle1, n2, 1);
0055     else 
0056         error('Both input should have same number of rows');
0057     end
0058 end
0059    
0060 % extract center and radius of each circle
0061 center1 = circle1(:, 1:2);
0062 center2 = circle2(:, 1:2);
0063 r1 = circle1(:,3);
0064 r2 = circle2(:,3);
0065 
0066 % allocate memory for result
0067 nPoints = length(r1);
0068 points = NaN * ones(2*nPoints, 2);
0069 
0070 % distance between circle centers
0071 d12 = distancePoints(center1, center2, 'diag');
0072 
0073 % get indices of circle couples with intersections
0074 inds = d12 >= abs(r1 - r2) & d12 <= (r1 + r2);
0075 
0076 if sum(inds) == 0
0077     return;
0078 end
0079 
0080 % angle of line from center1 to center2
0081 angle = angle2Points(center1(inds,:), center2(inds,:));
0082 
0083 % position of intermediate point, located at the intersection of the
0084 % radical axis with the line joining circle centers
0085 d1m  = d12(inds) / 2 + (r1(inds).^2 - r2(inds).^2) ./ (2 * d12(inds));
0086 tmp = polarPoint(center1(inds, :), d1m, angle);
0087 
0088 % distance between intermediate point and each intersection point
0089 h   = sqrt(r1(inds).^2 - d1m.^2);
0090 
0091 % indices of valid intersections
0092 inds2 = find(inds)*2;
0093 inds1 = inds2 - 1;
0094 
0095 % create intersection points
0096 points(inds1, :) = polarPoint(tmp, h, angle - pi/2);
0097 points(inds2, :) = polarPoint(tmp, h, angle + pi/2);