function curve = linearregress(X,Y,Z)
% LINEARREGRESS Fit a line by regression to data.
%   CURVE = LINEARREGRESS(X,Y) returns a vector of coefficients [c b]
%   which correspond to the linear function y = bx + c.  For N data
%   points, X is a 1-by-N vector of predictors and Y is a 1-by-N vector 
%   of observations.
%
%   CURVE = LINEARREGRESS(X,Y,Z) where Z is a vector of weights,
%   performs a weighted least squares regression.
%
% The coefficients in COEFF are with respect to the image coordinate
% system. 
%
%   See also PARABOLICREGRESS, REGRESS

% Copyright (c) 2011
% Jerod Weinman
% jerod@acm.org
  
% $Id: linearregress.m 595 2014-10-24 18:15:54Z weinman $
 
% Jerod Weinman, Zach Butler, Dugan Knoll, Alexander Turner, Yirong Jing
% jerod@acm.org
% Copyright 2011, 2012, 2013, 2014
 
%    This file is part of PIONEER WordNorm.
%
%    PIONEER WordNorm is free software: you can redistribute it and/or modify
%    it under the terms of the GNU General Public License as published by
%    the Free Software Foundation, either version 3 of the License, or
%    (at your option) any later version.
%
%    PIONEER WordNorm is distributed in the hope that it will be useful,
%    but WITHOUT ANY WARRANTY; without even the implied warranty of
%    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%    GNU General Public License for more details.
%
%    You should have received a copy of the GNU General Public License
%    along with PIONEER WordNorm.  If not, see <http://www.gnu.org/licenses/>.
  
% Subtract off the averages of the values
    x_mu = mean(X);
    X = X - x_mu;
    y_mu = mean(Y);
    Y = Y - y_mu;
    
    % Set up the X values for the regression  Y = bX + c
    Xfunc = [X ones([length(X) 1])];
    
    if nargin<3
      bc = regress(Y, Xfunc);
    else
      bc = lscov(Xfunc,Y,Z);
    end;

    % We found the regression with the origin being centered around the
    % middle of the data points.  Now we want to readjust the origin to be
    % where it was initially.
    
    curve(2) = bc(1);
    curve(1) = -bc(1)*x_mu + bc(2) + y_mu;
    
