% 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