% 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]
%
% XY can also be a 3xNpts array of homogeneous coordinates.
%
% 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-2008 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.
% November 2008 - Normalising of coordinates added to condition the solution.
function [C, dist] = fitline(XY)
[rows,npts] = size(XY);
if npts < 2
error('Too few points to fit line');
end
if rows ==2 % Add homogeneous scale coordinate of 1
XY = [XY; ones(1,npts)];
end
if npts == 2 % Pad XY with a third column of zeros
XY = [XY zeros(3,1)];
end
% Normalise points so that centroid is at origin and mean distance from
% origin is sqrt(2). This conditions the equations
[XYn, T] = normalise2dpts(XY);
% Set up constraint equations of the form XYn'*C = 0,
% where C is a column vector of the line coefficients
% in the form c(1)*X + c(2)*Y + c(3) = 0.
[u d v] = svd(XYn',0); % Singular value decomposition.
C = v(:,3); % Solution is last column of v.
% Denormalise the solution
C = T'*C;
% 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