function [N_features,ClusterA,ClusterB,h]=cp_moptxf(x,FeatureFN,Criterion,ImageId,NonlinId,SmoothId,NormId,FeatDec,N_tex,ClusterA,ClusterB,MN,sigma,h)
% CP_MOPTXF  Comparative texture measure study: Optimal filterbank features


% Initialization
if strcmp(lower(SmoothId),'gauss')
  if isempty(sigma)
    sigma=8;
  end
  f_0=1./(2*sqrt(2)*sigma);
  [dummy,w]=gausslp([],sigma);
elseif strcmp(lower(SmoothId),'average')
  M_w=31;
  N_w=31;
  f_0=[M_w N_w];
  w=ones(M_w,N_w)/(M_w*N_w);
else
  error(fprintf('Smoothing option "%s" not supported',SmoothId));
end
if ~strcmp(lower(NonlinId),'sqr')
  error(fprintf('Nonlinearity option "%s" not supported',NonlinId));
end
M_h=MN(1);
N_h=MN(2);

% Design filter
if isempty(h)
  % Clustering of the textures
  r_xx=[];
  x_design=[];
  for TexNo=1:N_tex
    % Filter design
    [x_i,N_tex,r_xxi]=cp_image(ImageId,TexNo);
    r_xx=[r_xx ; r_xxi];
    x_design=[x_design ; x_i(1:128,1:128)];
  end
  % Filter design
  [h,ClusterA,ClusterB]=moptxf(x_design,r_xx,M_h,N_h,Criterion,w,...
    N_tex,[],[],[],7,7);
end
N_features=size(ClusterA,1);
fprintf(1,'(%d filter(s))',N_features);

% Feature extraction
for n_h=1:N_features
  fprintf(1,'.');
  h_n=h((1:M_h)+(n_h-1)*M_h,1:N_h);
  y=filter2(h_n,x);   % Filtering
  z=cp_energy(y,NonlinId,SmoothId,NormId,f_0);  % Local energy function
  catframe(z,FeatureFN,[],'Interm',FeatDec);  % Write feature file
end



