function [J,G]=multifilt_gopt_J(h_,M_h,N_h,Criterion,m_xA,m_xB,r_xxA,r_xxB,w,P_save)

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),var_A(i)]=fstats(m_xA(i),r_xxA((i-1)*M_r+1:i*M_r,:),h,w);
  end
  for i=1:length(m_xB)
    [m_B(i),var_B(i)]=fstats(m_xB(i),r_xxB((i-1)*M_r+1:i*M_r,:),h,w);
  end
  % Compute the overall means and variances
  m_oA=mean(m_A);
  m_oB=mean(m_B);
  var_oA=mean(var_A) + mean((m_A-m_oA).^2);
  var_oB=mean(var_B) + mean((m_B-m_oB).^2);
  % Compute the criterion function
  J=-(m_oA - m_oB)^2 / (var_oA + var_oB);
elseif strcmp(lower(Criterion),'meanscatter')
  % Estimate the means
  m_A=[]; m_B=[];
  for i=1:length(m_xA)
    m_A(i)=fstats(m_xA(i),r_xxA((i-1)*M_r+1:i*M_r,:),h,w);
  end
  for i=1:length(m_xB)
    m_B(i)=fstats(m_xB(i),r_xxB((i-1)*M_r+1:i*M_r,:),h,w);
  end
  % Compute key elements of the criterion function
  m_oA=mean(m_A);
  m_oB=mean(m_B);
  % Compute the criterion function
  J=-(m_oA - m_oB)^2 / (mean((m_A-m_oA).^2) + mean((m_B-m_oB).^2));
else
  error(sprintf('Unsupported criterion "%s"',Criterion));  
end

% Save intermediate results with probability P_save
if exist('P_save')~=1
  P_save=0;
end
if rand<=P_save
  save multifilt_gopt_h.mat h J
end

% Compute the constraint function
G=(1-norm(h_)^2);
