lineFit

LINEFIT Fit a straight line to a set of points.

   L = lineFit(X, Y)
   Computes parametric line minimizing square error of all points (X,Y).
   Result is a 4*1 array, containing coordinates of a point of the line,
   and the direction vector of the line, that is  L=[x0 y0 dx dy];

   L = lineFit(PTS) 
   Gives coordinats of points in a single array.

   L = lineFit(PT0, PTS);
   L = lineFit(PT0, X, Y);
   with PT0 = [x0 y0], imposes the line to contain point PT0.

   Requires:
   Optimiaztion toolbox

   See also:
   lines2d, polyfit, polyfit2, lsqlin


   -----
   author : David Legland
   INRA - TPV URPOI - BIA IMASTE
   created the 30/04/2004.

0001 function line = lineFit(varargin)
0002 %LINEFIT Fit a straight line to a set of points.
0003 %
0004 %   L = lineFit(X, Y)
0005 %   Computes parametric line minimizing square error of all points (X,Y).
0006 %   Result is a 4*1 array, containing coordinates of a point of the line,
0007 %   and the direction vector of the line, that is  L=[x0 y0 dx dy];
0008 %
0009 %   L = lineFit(PTS)
0010 %   Gives coordinats of points in a single array.
0011 %
0012 %   L = lineFit(PT0, PTS);
0013 %   L = lineFit(PT0, X, Y);
0014 %   with PT0 = [x0 y0], imposes the line to contain point PT0.
0015 %
0016 %   Requires:
0017 %   Optimiaztion toolbox
0018 %
0019 %   See also:
0020 %   lines2d, polyfit, polyfit2, lsqlin
0021 %
0022 %
0023 %   -----
0024 %   author : David Legland
0025 %   INRA - TPV URPOI - BIA IMASTE
0026 %   created the 30/04/2004.
0027 %
0028 
0029 %   HISTORY
0030 %   09/12/2004 : update implementation
0031 
0032 
0033 
0034 % ---------------------------------------------
0035 % extract input arguments
0036 
0037 if length(varargin)==1
0038     % argument is an array of points
0039     var = varargin{1};
0040     x = var(:,1);
0041     y = var(:,2);   
0042 elseif length(varargin)==2
0043     var = varargin{1};
0044     if size(var, 1)==1
0045         var = varargin{2};
0046         x = var(:,1);
0047         y = var(:,2);
0048     else
0049         % two arguments : x and y
0050         x = var;
0051         y = varargin{2};
0052     end
0053 elseif length(varargin)==3
0054     % three arguments : ref point, x and y
0055     x = varargin{2};
0056     y = varargin{3};
0057 end
0058     
0059 % ---------------------------------------------
0060 % Initializations :
0061 
0062 
0063 N = size(x, 1);
0064 
0065 % ---------------------------------------------
0066 % Main algorithm :
0067 
0068 
0069 % main matrix of the problem
0070 X = [x y ones(N,1)];
0071 
0072 % conditions initialisations
0073 A = zeros(0, 3);
0074 b = [];
0075 Aeq1 = [1 1 0];
0076 beq1 = 1;
0077 Aeq2 = [1 -1 0];
0078 beq2 = 1;
0079 
0080 % disable verbosity of optimisation
0081 opt = optimset('lsqlin');
0082 opt.LargeScale = 'off';
0083 opt.Display = 'off';
0084 
0085 % compute line coefficients [a;b;c] , in the form a*x + b*y + c = 0
0086 % using linear regression
0087 % Not very clean : I could not impose a*a+b*b=1, so I checked for both a=1
0088 % and b=1, and I kept the result with lowest residual error....
0089 [coef1, res1] = lsqlin(X, zeros(N, 1), A, b, Aeq1, beq1, [], [], [], opt);
0090 [coef2, res2] = lsqlin(X, zeros(N, 1), A, b, Aeq2, beq2, [], [], [], opt);
0091 
0092 % choose the best line
0093 if res1<res2
0094     coef = coef1;
0095 else
0096     coef = coef2;
0097 end
0098 
0099 line = cartesianLine(coef');