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=fstats(mean_x,r_xx,h) estimates the mean feature value
%	[fmean,fvar]=fstats(mean_x,r_xx,h,w) estimates the mean and variance
%	[fmean,fvar]=fstats(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=sum(w(:)) * h_'*r2R(r_xx,M_h,N_h)*h_;

% Estimate variance (if requested)
if nargout>1
  [M_w,N_w]=size(w);
  w_=reshape(w',M_w*N_w,1);
  % Estimate r_yy (computationally intensive)
  
  r_yy=NaN*zeros(2*M_w-1,2*N_w-1);
  r_yy_IsSymmetric=isreal(h) & isreal(r_xx);
  for k=-M_w+1:M_w-1
    for l=-N_w+1:N_w-1
      if isnan(r_yy(M_w+k,N_w+l))
	r_yy(M_w+k,N_w+l)=h_'*r2R(r_xx,M_h,N_h,k,l)*h_;
	if r_yy_IsSymmetric
	  r_yy(M_w-k,N_w-l)=r_yy(M_w+k,N_w+l);
	end
      end
    end
  end

  % Estimate r_zz
  mean_y=mean_x*sum(h_);
  r_zz=2*r_yy.^2-2*mean_y^4+r_yy(M_w,N_w).^2;
  
  if AR1
    % Estimate AR(1) parameters
    rho(1)=mean(mean(r_zz(M_w+1:2*M_w-1,:)./r_zz(M_w:2*M_w-2,:)));
    rho(2)=mean(mean(r_zz(:,N_w+1:2*N_w-1)./r_zz(:,N_w:2*N_w-2)));
    [l,k]=meshgrid(-N_w+1:N_w-1,-M_w+1:M_w-1); l=abs(l); k=abs(k);
    r_zz_ar=r_zz(M_w,N_w).*rho(1).^k.*rho(2).^l;
  
    % Normalize
    weight=zeros(size(r_zz_ar));
    for k=-M_w+1:M_w-1
      for l=-N_w+1:N_w-1
	sift=zeros(size(r_zz_ar)); sift(M_w+k,N_w+l)=1;
	weight(M_w+k,N_w+l)=w_'*r2R(sift,M_w,N_w)*w_;
      end
    end
    weight=weight/mean(weight(:));
    r_zz=r_zz_ar/mean(mean(weight.*r_zz_ar./r_zz)); clear r_zz_ar
  end
    
  % Variance estimation
  if M_w*N_w>41^2
    % Slow technique, can be necessary for very large smoothing
    % filters. Size limit where it is necessary depends on memory
    % configuration of the computer. The speed problems are due to Matlab
    % command interpretation.
    fvar=-fmean.^2;
    for k_1=0:M_w-1
      for l_1=0:N_w-1
	for k_2=0:M_w-1
	  for l_2=0:N_w-1
	    fvar=fvar+...
	      w(k_1+1,l_1+1)*w(k_2+1,l_2+1)*r_zz(M_w+k_1-k_2,N_w+l_1-l_2);
	  end
	end
      end
    end
  else
    % Significantly faster technique (Matlab)
    fvar=w_'*r2R(r_zz,M_w,N_w)*w_-fmean.^2;
  end

end
