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 ------ Author: David Legland e-mail: david.legland@grignon.inra.fr Created: 2011-01-20, using Matlab 7.9.0.529 (R2009b) Copyright 2011 INRA - Cepia Software Platform.
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@grignon.inra.fr 0044 % Created: 2011-01-20, using Matlab 7.9.0.529 (R2009b) 0045 % Copyright 2011 INRA - Cepia Software Platform. 0046 0047 % HISTORY 0048 % 2011-06-06 add support for multiple inputs 0049 0050 0051 % adapt sizes of inputs 0052 n1 = size(circle1, 1); 0053 n2 = size(circle2, 1); 0054 if n1 ~= n2 0055 if n1 > 1 && n2 == 1 0056 circle2 = repmat(circle2, n1, 1); 0057 elseif n2 > 1 && n1 == 1 0058 circle1 = repmat(circle1, n2, 1); 0059 else 0060 error('Both input should have same number of rows'); 0061 end 0062 end 0063 0064 % extract center and radius of each circle 0065 center1 = circle1(:, 1:2); 0066 center2 = circle2(:, 1:2); 0067 r1 = circle1(:,3); 0068 r2 = circle2(:,3); 0069 0070 % allocate memory for result 0071 nPoints = length(r1); 0072 points = NaN * ones(2*nPoints, 2); 0073 0074 % distance between circle centers 0075 d12 = distancePoints(center1, center2, 'diag'); 0076 0077 % get indices of circle couples with intersections 0078 inds = d12 >= abs(r1 - r2) & d12 <= (r1 + r2); 0079 0080 if sum(inds) == 0 0081 return; 0082 end 0083 0084 % angle of line from center1 to center2 0085 angle = angle2Points(center1(inds,:), center2(inds,:)); 0086 0087 % position of intermediate point, located at the intersection of the 0088 % radical axis with the line joining circle centers 0089 d1m = d12(inds) / 2 + (r1(inds).^2 - r2(inds).^2) ./ (2 * d12(inds)); 0090 tmp = polarPoint(center1(inds, :), d1m, angle); 0091 0092 % distance between intermediate point and each intersection point 0093 h = sqrt(r1(inds).^2 - d1m.^2); 0094 0095 % indices of valid intersections 0096 inds2 = find(inds)*2; 0097 inds1 = inds2 - 1; 0098 0099 % create intersection points 0100 points(inds1, :) = polarPoint(tmp, h, angle - pi/2); 0101 points(inds2, :) = polarPoint(tmp, h, angle + pi/2);