function [h,ClusterA,ClusterB]=moptxf(x,r_xx,M_h,N_h,Criterion,w,ClusterA,ClusterB,MaxNumIter,h,M_h_cluster,N_h_cluster,N_T)
% MOPTXF  Multiple optimal texture filters design 
%       [h,ac,cb]=moptxf([x_1;...;x_m],[r_xx1;...;r_xxm],M_h,N_h,Criterion,w,m)
%       clusters the m textures and designs filters for each cluster pair.
%       x_1...x_m are the texture images and r_xx1...r_xxm are their
%       autocorrelation functions (both textures and correlations must be
%       of identical size). The filters, h, are of size M_h by N_h, and 
%       designed according to Criterion (see below). All filters are stacked
%       in the matrix h. w is the smoothing filter, and m is the number of 
%       textures. The cluster pairs are returned in ca and cb.
%       
%       h=moptxf([x_1;...;x_m],[r_xx1;...;r_xxm],M_h,N_h,Criterion,w,ca,cb)
%       uses the supplied clusters ca and cb, and determines the number
%       of textures from the largest element in ca and cb.
%       
%       ...=moptxf(...cb,MaxNumIter,h) will, if a gradient search criterion
%       is selected, use the h input as initial filters, and will not iterate 
%       more than MaxNumIter.
%       
%       Supported criteria are:
%         'Unser'       Closed form 
%         'Singh'       Closed form
%         'MeanScatter' Relatively fast gradient search solution
%         'FisherG'     Complex gradient search solution

if exist('ClusterB')~=1,    ClusterB=[];     end
if exist('M_h_cluster')~=1, M_h_cluster=[];  end
if exist('N_h_cluster')~=1, N_h_cluster=[];  end
if exist('N_T')~=1,         N_T=[];          end
if isempty(M_h_cluster),    M_h_cluster=M_h; end
if isempty(N_h_cluster),    N_h_cluster=N_h; end

if isempty(N_T)
  if isempty(ClusterB)
    N_T=ClusterA;
  else
    N_T=max(max([ClusterA ClusterB]));    % Number of textures
  end
end
[M_r,N_r]=size(r_xx); M_r=M_r/N_T;      % Size of r_xxi
[M_x,N_x]=size(x);    M_x=M_x/N_T;      % Size of x_i  
if isempty(ClusterB)
  N_T=ClusterA;
  Distance=txtursim([],r_xx,N_T,Criterion,M_h,N_h);
  [ClusterA,ClusterB]=txturgroup(Distance,mean(reshape(x.',M_x*N_x,N_T))',...
    r_xx,M_h_cluster,N_h_cluster,w,'MeanSep',1.2,[],2,1);
end
N_H=size(ClusterA,1); % Number of filters to be designed

for i=1:N_H
  ClusterAA=nonzero(ClusterA(i,:));
  ClusterBB=nonzero(ClusterB(i,:));
  % Filter optimization 
  x_A=x((1:M_x)'*ones(1,length(ClusterAA)) ...
    +M_x*(ClusterAA(ones(1,M_x),:)-1),:); 
  x_B=x((1:M_x)'*ones(1,length(ClusterBB)) ....
    +M_x*(ClusterBB(ones(1,M_x),:)-1),:); 
  if  strcmp(lower(Criterion),'unser') | strcmp(lower(Criterion),'singh')
    r_xxA=eval(sprintf('+r_xx((1:M_r)+%d,:)',(ClusterAA-1)*M_r))...
      /length(ClusterAA);
    r_xxB=eval(sprintf('+r_xx((1:M_r)+%d,:)',(ClusterBB-1)*M_r))...
      /length(ClusterBB);
  else
    r_xxA=eval(sprintf('[%s]',...
      sprintf(';r_xx((1:M_r)+%d,:)',(ClusterAA-1)*M_r)));
    r_xxB=eval(sprintf('[%s]',...
      sprintf(';r_xx((1:M_r)+%d,:)',(ClusterBB-1)*M_r)));
  end
  if  strcmp(lower(Criterion),'unser') ...
      | strcmp(lower(Criterion),'singh')
    h((1:M_h)+(i-1)*M_h,1:N_h) ...
      =optxf(Criterion,x_A,x_B,[M_h N_h],w,r_xxA,r_xxB);
  elseif strcmp(lower(Criterion),'fisher') ...
      |    strcmp(lower(Criterion),'fisherg')
    % Set gradient search options
    opt=[];
    w_1=ones(5,5)/5^2;    opt(1,1)=1; opt(1,3)=1e-3; opt(1,14)=300;
    w_2=ones(9,9)/9^2;    opt(2,1)=1; opt(2,3)=1e-3; opt(2,14)=150;
    w_3=ones(13,13)/13^2; opt(3,1)=1; opt(3,3)=1e-3; opt(3,14)=75;
    opt(4,1)=1; opt(4,3)=1e-3; opt(4,14)=MaxNumIter;
    h((1:M_h)+(i-1)*M_h,:)=multifilt_gopt(h((1:M_h)+(i-1)*M_h,:),Criterion,...
      mean(reshape(x_A.',M_x*N_x,length(ClusterAA))),...
      mean(reshape(x_B.',M_x*N_x,length(ClusterBB))),...
      r_xxA,r_xxB,w,0, ...
      opt(size(opt,1),:),w_1,opt(1,:),w_2,opt(2,:),w_3,opt(3,:));
  elseif strcmp(lower(Criterion),'meanscatter') 
    if length(ClusterAA)>1 | length(ClusterBB)>1
      opt=1; opt(3)=1e-3; opt(14)=MaxNumIter;
      h((1:M_h)+(i-1)*M_h,:)=multifilt_gopt(h((1:M_h)+(i-1)*M_h,:),...
        Criterion,...
        mean(reshape(x_A.',M_x*N_x,length(ClusterAA))),...
        mean(reshape(x_B.',M_x*N_x,length(ClusterBB))),...
        r_xxA,r_xxB,w,0,opt);
    else
      % Closed form optimal two-texture filter already found
    end
  else
    error(sprintf('Criterion "%s" not supported',Criterion));
  end
  eval(sprintf('h_%d=h;',i-1));
end