% FITLINE - Least squares fit of a line to a set of points%% Usage:   [C, dist] = fitline(XY)%% Where:   XY  - 2xNpts array of xy coordinates to fit line to data of%                the form %                [x1 x2 x3 ... xN%                 y1 y2 y3 ... yN]%% Returns: C    - 3x1 array of line coefficients in the form%                 c(1)*X + c(2)*Y + c(3) = 0%          dist - Array of distances from the fitted line to the supplied%                 data points.  Note that dist is only calculated if the%                 function is called with two output arguments.%% The magnitude of C is scaled so that line equation corresponds to%   sin(theta)*X + (-cos(theta))*Y + rho = 0% where theta is the angle between the line and the x axis and rho is the% perpendicular distance from the origin to the line.  Rescaling the% coefficients in this manner allows the perpendicular distance from any% point (x,y) to the line to be simply calculated as%   r = abs(c(1)*X + c(2)*Y + c(3))%%% If you want to convert this line representation to the classical form %     Y = a*X + b% use%  a = -c(1)/c(2)%  b = -c(3)/c(2)%% Note, however, that this assumes c(2) is not zero% Copyright (c) 2003-2005 Peter Kovesi% School of Computer Science & Software Engineering% The University of Western Australia% http://www.csse.uwa.edu.au/% % Permission is hereby granted, free of charge, to any person obtaining a copy% of this software and associated documentation files (the "Software"), to deal% in the Software without restriction, subject to the following conditions:% % The above copyright notice and this permission notice shall be included in % all copies or substantial portions of the Software.%% The Software is provided "as is", without warranty of any kind.% June      2003 - Original version% September 2004 - Rescaling to allow simple distance calculation.function [C, dist] = fitline(XY)    [rows,npts] = size(XY);        if rows ~=2    error('Data is not 2D');  end    if npts < 2    error('Too few points to fit line');  end    % Set up constraint equations of the form  AC = 0,  % where C is a column vector of the line coefficients  % in the form   c(1)*X + c(2)*Y + c(3) = 0.    A = [XY' ones(npts,1)];  % Build constraint matrix    if npts == 2             % Pad A with zeros    A = [A; zeros(1,3)];   end  [u d v] = svd(A);        % Singular value decomposition.  C = v(:,3);              % Solution is last column of v.  % Rescale coefficients so that line equation corresponds to  %   sin(theta)*X + (-cos(theta))*Y + rho = 0  % so that the perpendicular distance from any point (x,y) to the line  % to be simply calculated as   %   r = abs(c(1)*X + c(2)*Y + c(3))  theta = atan2(C(1), -C(2));  % Find the scaling (but avoid dividing by zero)  if abs(sin(theta)) > abs(cos(theta))      k = C(1)/sin(theta);  else      k = -C(2)/cos(theta);  end    C = C/k;    % If requested, calculate the distances from the fitted line to  % the supplied data points   if nargout==2         dist = abs(C(1)*XY(1,:) + C(2)*XY(2,:) + C(3));  end