function theta = LS_est (y,N)
% theta = LS_est (y,N)
%
% Generate the non-symmetric least squares noncausal parameter estimate 
% based on the input image y and the neighbor set N.
%
%
% N is given by a matrix with two lines, first line is relative model position
% in vertical direction, second is corresponding column.
% For example
%
%   N= [-1	-1	0
%	-1	0	1]
%
% means a neighbor set of {(-1,-1),(-1,0),(0,1)}, which in turn means
% the GMRF model
%   y(i,j)=theta(-1,-1) y(i-1,j-1) + theta(-1,0) y(i-1,j) + theta(0,1) y(i,j+1)

% Compute some basic model parameters
M=size(N,2);				% Model order
vs=size(y,1);				% Image vertical size
hs=size(y,2);				% Image horizontal size

% Convert the input image y to a vector, line by line from y
y=reshape(y',size(y,1)*size(y,2),1);

% Generate the matrix YY=Y'*Y without generating the huge matrix Y.
for i=1:M,
  for j=1:M,
    YY(i,j) = vshift(y,N(1,i)*hs+N(2,i))' ...
	* vshift(y,N(1,j)*hs+N(2,j));
  end
end

% Generate the matrix Yy=Y'*y without generating the huge matrix Y.
for i=1:M,
  Yy(i,1) = vshift(y,N(1,i)*hs+N(2,i))' * y;
end

% Solve the normal equations: Y'*Y*theta=Y'*y
theta=YY\Yy;

% Export the parameters in packed parameter format:
theta=[N
       theta'];
