function [fmean,fvar]=fstats(mean_x,r_xx,h,w,P1)
% FVAR	Estimate the statistcs of extracted texture features
%	The modelled feature extraction process constitutes
%	 1) filtering by h (two-dimensional)
%	 2) squaring
%	 3) smoothing by the filter w (two-dimensional)
%	
%	fmean=ftstats2(mean_x,r_xx,h) estimates the mean feature value
%	[fmean,fvar]=ftstats2(mean_x,r_xx,h,w) estimates the mean and variance
%	[fmean,fvar]=ftstats2(mean_x,r_xx,h,w,'AR(1)') estimates the mean and 
%	variance using the separable AR(1) assumption

AR1=0;
if exist('P1')==1
  if strcmp(upper(P1),'AR(1)')
    AR1=1;
  end
end


% Some initial assignments etc.
[M_h,N_h]=size(h);
h_=reshape(h',M_h*N_h,1);

% Estimate mean
fmean=h_'*r2R(r_xx,M_h,N_h)*h_;

% Estimate variance (if requested)
if nargout>1
  [M_w,N_w]=size(w);
  fvar=0;
  for beta=-N_w+1:N_w-1
    for alpha=-M_w+1:M_w-1
      fvar=fvar ...
	+(M_w-abs(alpha))*(N_w-abs(beta)) ...
	*(h_'*r2R(r_xx,M_h,N_h,alpha,beta)*h_)^2;
    end
  end
  fvar=fvar * 2/(M_w^2*N_w^2) - 2*(mean_x*sum(h_))^4;
end
