function [DJ,DG]=multifilt_gopt_DJ(h_,M_h,N_h,Criterion,m_xA,m_xB,r_xxA,r_xxB,w,dummy)

h=reshape(h_,N_h,M_h)';
[M_r,N_r]=size(r_xxA); M_r=M_r/length(m_xA);

if strcmp(lower(Criterion),'fisher') | strcmp(lower(Criterion),'fisherg')
  % Estimate the means and variances
  m_A=[]; var_A=[]; m_B=[]; var_B=[];
  for i=1:length(m_xA)
    [m_A(i),Dm,var_A(i),Dvar]=Dfstats(m_xA(i),r_xxA((i-1)*M_r+1:i*M_r,:),h,w);
    Dm_A(1:M_h*N_h,i)=reshape(Dm.',M_h*N_h,1);
    Dvar_A(1:M_h*N_h,i)=reshape(Dvar.',M_h*N_h,1);
  end
  for i=1:length(m_xB)
    [m_B(i),Dm,var_B(i),Dvar]=Dfstats(m_xB(i),r_xxB((i-1)*M_r+1:i*M_r,:),h,w);
    Dm_B(1:M_h*N_h,i)=reshape(Dm.',M_h*N_h,1);
    Dvar_B(1:M_h*N_h,i)=reshape(Dvar.',M_h*N_h,1);
  end
  % Compute key elements of the criterion function and its derivative
  m_oA=mean(m_A); Dm_oA=mean(Dm_A.').'; 
  m_oB=mean(m_B); Dm_oB=mean(Dm_B.').'; 
  sc_A=(m_A-m_oA).'; sc_A=sc_A(:,ones(1,M_h*N_h));
  sc_B=(m_B-m_oB).'; sc_B=sc_B(:,ones(1,M_h*N_h));
  var_oA=mean(var_A) + mean((m_A-m_oA).^2);
  var_oB=mean(var_B) + mean((m_B-m_oB).^2);
  Dvar_oA=mean(Dvar_A.').' ...
    + mean(2*sc_A.*(Dm_A-Dm_oA(:,ones(1,length(m_xA)))).' ).';
  Dvar_oB=mean(Dvar_B.').' ...
    + mean(2*sc_B.*(Dm_B-Dm_oB(:,ones(1,length(m_xB)))).' ).';
  % Compute the criterion function
  DJ=-2*(m_oA - m_oB) * (Dm_oA - Dm_oB)  / (var_oA + var_oB) ...
    + (m_oA - m_oB)^2 * (Dvar_oA + Dvar_oB)  / (var_oA + var_oB)^2;
elseif strcmp(lower(Criterion),'meanscatter')
  % Estimate the means and their derivatives
  m_A=[]; m_B=[]; Dm_A=[]; Dm_B=[];
  for i=1:length(m_xA)
    [m_A(i),Dm]=Dfstats(m_xA(i),r_xxA((i-1)*M_r+1:i*M_r,:),h,w);
    Dm_A(1:M_h*N_h,i)=reshape(Dm.',M_h*N_h,1);
  end
  for i=1:length(m_xB)
    [m_B(i),Dm]=Dfstats(m_xB(i),r_xxB((i-1)*M_r+1:i*M_r,:),h,w);
    Dm_B(1:M_h*N_h,i)=reshape(Dm.',M_h*N_h,1);
  end
  % Compute key elements of the criterion function and its derivative
  m_oA=mean(m_A); Dm_oA=mean(Dm_A.').'; 
  m_oB=mean(m_B); Dm_oB=mean(Dm_B.').'; 
  sc_A=(m_A-m_oA).'; sc_A=sc_A(:,ones(1,M_h*N_h));
  sc_B=(m_B-m_oB).'; sc_B=sc_B(:,ones(1,M_h*N_h));
  % Compute the derivative of the criterion function
  DJ=-2*(m_oA - m_oB) * (Dm_oA - Dm_oB) ...
    / (mean((m_A-m_oA).^2) + mean((m_B-m_oB).^2)) ...
    + (m_oA - m_oB)^2 ...
    * (mean(2*sc_A.*(Dm_A-Dm_oA(:,ones(1,length(m_xA)))).' ).' ...
    +  mean(2*sc_B.*(Dm_B-Dm_oB(:,ones(1,length(m_xB)))).' ).') ...
    / (mean((m_A-m_oA).^2) + mean((m_B-m_oB).^2))^2;
else
  error(sprintf('Unsupported criterion "%s"',Criterion));  
end

DJ=reshape(DJ.',M_h*N_h,1);

DG=-2*h_;

