function [fmean,Dfmean,fvar,Dfvar]=Dfstats(mean_x,r_xx,h,w)
% FVAR	Estimate the statistics and the  differentiated statistcs of features
%	The modelled texture feature extraction process constitutes
%	 1) filtering by h (two-dimensional)
%	 2) squaring
%	 3) smoothing by the filter w (two-dimensional)
%	
%	[fmean,Dfmean]=Dfstats(mean_x,r_xx,h,w) estimates the mean and 
%	the derivative of the mean with respect to h
%	
%	[fmean,Dfmean,fvar,Dfvar]=Dfstats(mean_x,r_xx,h,w) 
%	estimates the mean and variance and the derivative of the mean 
%	and variance with respect to h


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

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

if nargout>2
  % Estimate differentiated variance
  [M_w,N_w]=size(w);
  w_=reshape(w',M_w*N_w,1);
  % Estimate r_yy and Dr_yy (computationally intensive)
  r_yy=NaN*zeros(2*M_w-1,2*N_w-1);
  Dr_yy=NaN*zeros(M_h*N_h,(2*M_w-1)*(2*N_w-1));
  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_xx=r2R(r_xx,M_h,N_h,k,l);
	r_yy(M_w+k,N_w+l)=h_'*R_xx*h_;
	Dr_yy(:,(M_w+k-1)*(2*N_w-1)+N_w+l)=R_xx*h_ + R_xx'*h_;
	if isreal(h) & isreal(r_xx)
	  r_yy(M_w-k,N_w-l)=r_yy(M_w+k,N_w+l);
	  Dr_yy(:,(M_w-k-1)*(2*N_w-1)+N_w-l) ...
	    =Dr_yy(:,(M_w+k-1)*(2*N_w-1)+N_w+l);
	end
      end
    end
  end

  % Estimate r_zz and Dr_zz
  r_zz=2*r_yy.^2-2*mean_y^4+r_yy(M_w,N_w).^2;
  for k=-M_w+1:M_w-1
    for l=-N_w+1:N_w-1
      Dr_zz(:,(M_w+k-1)*(2*N_w-1)+N_w+l)=...
	4*r_yy(M_w+k,N_w+l)*Dr_yy(:,(M_w+k-1)*(2*N_w-1)+N_w+l) ...
	- 8*mean_y^3*mean_x*ones(M_h*N_h,1) ...
	+ 2*r_yy(M_w,N_w)*Dr_yy(:,(M_w-1)*(2*N_w-1)+N_w);
    end
  end

  % Variance estimation
  fvar=w_'*r2R(r_zz,M_w,N_w)*w_-fmean.^2;
  for i=1:M_h*N_h
    Dfvar(i) ...
      =w_'*r2R(reshape(Dr_zz(i,:),2*N_w-1,2*M_w-1).',M_w,N_w)*w_ ...
      -2*fmean*Dfmean(i); 
  end
  Dfvar =reshape(Dfvar ,N_h,M_h).';
end

Dfmean=reshape(Dfmean,N_h,M_h).';

