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

% 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 Y
Y=vshift(y,N(1,1)*hs+N(2,1)) + vshift(y,-N(1,1)*hs-N(2,1));
for i=2:M,
  Y=[Y vshift(y,N(1,i)*hs+N(2,i))+vshift(y,-N(1,i)*hs-N(2,i))];
end

% Solve the normal equations: Y'*Y*theta=Y'*y
theta=(Y'*Y)\(Y'*y);

% Calculate mean square error
mse=mean((y-Y*theta).^2);

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

