points2circle(P,P2,P3)

function xyr = points2circle(P,P2,P3)
% POINTS2CIRCLE - calculate center and radius of circle through three points 
%
%   XYR = POINTS2CIRCLE(P1,P2,P3) or POINTS2CIRCLE(P) calculates the center
%   [XYR(1),XYR(2)] and the radius XYR(3) of the circle passing through
%   three points specified by [Pi(1), Pi(2)]. P1, P2 and P3 are vectors
%   with two elements, or P is a 3-by-2 matrix where each row specifies a
%   point. If points are on a straight line (collinear), a warning is given
%   and XYR will have three NaNs.  
%
%   Example:
%
%   % DATA
%       px = [1 -8 3] ; py = [12 3 2] ;
%   % ENGINE
%       xyr = points2circle([px(:) py(:)]) ;
%       % or points2circle([px(1) py(1)], [px(2) py(2)], [px(3) py(3)]) ;
%   % PLOT RESULTS
%     % create circle
%       t = linspace(0,2*pi,100) ; 
%       xc = xyr(1) + xyr(3) * cos(t) ;  yc = xyr(2) + xyr(3) * sin(t) ;
%     % make lines connecting center and points
%       x00 = [xyr([1 1 1]) ; px ; xyr([1 1 1])] ;
%       y00 = [xyr([2 2 2]) ; py ; xyr([2 2 2])] ;
%     % plot it nicely
%       h = plot(x00(:),y00(:),'k-', ...       % connect center to points
%          px([1:3 1]),py([1:3 1]),'b:', ...   % triangle
%          xc,yc,'r-', ...                     % circle
%          xyr(1),xyr(2),'r.', ...             % center
%          px,py,'bo' ...                      % the three points
%          ) ;
%       legend(h([5 3 4]),{'Points','Circle','Center'}) ;
%       axis square
%
%   See also DET

% for Matlab R13 and up
% version: 2.0 (mar 2008)
% (c) Jos van der Geest
% email: [e-mail]

% History
% Created: dec 2006, based on http://home.att.net/~srschmitt/circle3pts.html
% Revisions
% 1.1 (mar 2008) : modified for the file exchange (added help, changed
%                  variable names, small coding changes)
% 1.2 (mar 2008) : fixed some spelling errors in the help, add comments
%                  improved test for collinearity using rcond
% 2.0 (mar 2008) : shortened algorithm without local minor function

if nargin==1,
    if size(P,1) ~=3 | size(P,2) ~=2,
        error ('A single input should be a 3-by-2 matrix') ;
    end
elseif ~isequal(numel(P),numel(P2),numel(P3),2),
    error('The three points should all have 2 elements.')
else
    P = [P(:) P2(:) P3(:)] .' ; 
end

% build matrix
M = [sum(P.^2,2) P ones(3,1)] ;

if rcond(M(:,2:end)) < eps, 
    % test for collinearity
    xyr = [NaN NaN NaN] ;
    warning([mfilename ':CollinearPoints'], ...
        'No solution! Points may be on a straight line (collinear).') ;
else
    M11 = det(M(:,2:4)) ; % minor
    % allocate three element vector using NaN as a placeholder
    % calculate center
    xyr = [([det(M(:,[1 3 4])) -det(M(:,[1 2 4]))] ./ (2*M11)) NaN];
    % calculate radius
    xyr(3) = sqrt(sum(xyr(1:2).^2) + (det(M(:,[1 2 3])) ./ M11)) ; 
end

% Algorithm. By creating the matrix:
%
%    M = [ 1              1   1  1 
%          (x1^2 + y1^2) x1  y1  1  
%          (x2^2 + y2^2) x2  y2  1  
%          (x3^2 + y3^2) x3  y3  1  ]
%
% the following formulae can be calculated:
%
%    x0 =  M12 / (2*M11) 
%    y0 = -M13 / (2*M11) 
%    R  = sqrt(x0^2 + y0^2 + M14 / M11)) 
%
% in which Mij is the minor of M, i.e., the determinant after removal of
% the i-th row and j-th column of M. Note that the first row is always
% removed in the algorithm.