function [N_features,ClusterA,ClusterB,p1,p2,p3,p4]=cp_feature(x,FeatureFN,ImageId,MethodId,NonlinId,SmoothId,NormId,FeatDec,N_tex,ClusterA,ClusterB,p1,p2,p3,p4)
% CP_FEATURE  Compute features for texture comparison experiment
%	[N_features,ClusterA,ClusterB,p1,p2,p3,p4] ...
%	  =cp_feature(x,FeatureFN,ImageId,...
%	  MethodId,NonlinId,SmoothId,NormId,FeatDec,N_tex,...
%	  ClusterA,ClusterB,p1,p2,p3,p4)
%	
%	PS! It is required from most subroutines that they generate files
%	FeatureFN_001, FeatureFN_002, etc., by means of catframe(...,'Interm').

% Initialization
N_features=0;
fprintf(1,'Feature extr: ');
if exist('ClusterA')~=1, ClusterA=[]; end
if exist('ClusterB')~=1, ClusterB=[]; end
for k=1:4, 
  if exist(sprintf('p%d',k))~=1;
    eval(sprintf('p%d=[];',k));
  end
end

if idcmp(MethodId,'QMF')
  % QMF FB, Wavelet Transforms/Packets/Frames, DCT
  SepPos=findstr(MethodId,':');
  FilterId=MethodId(SepPos(1)+1:SepPos(2)-1);
  Structure=str2num(MethodId(SepPos(2)+1:SepPos(3)-1));
  DoDecFB=strcmp(lower(MethodId(SepPos(3)+1:length(MethodId))),'dec');
  N_features=cp_qmf(x,FeatureFN,NonlinId,SmoothId,NormId,FeatDec,...
    FilterId,Structure,DoDecFB);
elseif  idcmp(MethodId,'cooc')
  % Haralick Co-Occurrence
  if length(MethodId)==4
    BlockSize=16;
  else 
    BlockSize=str2num(MethodId(6:length(MethodId))); 
  end
  N_features=cp_cooc(x,FeatureFN,BlockSize,8,FeatDec);
elseif  idcmp(MethodId,'SAR')
  % Mao's Simultaneous Autoregressive Features
 N_features=cp_sar(x,FeatureFN,NonlinId,SmoothId,NormId,FeatDec);
elseif  idcmp(MethodId,'DCT')
  % Laws' texture masks
  N_features=cp_dct(x,FeatureFN,NonlinId,SmoothId,NormId,FeatDec);
elseif  idcmp(MethodId,'Laws')
  % Laws' texture masks
  N_features=cp_laws(x,FeatureFN,NonlinId,SmoothId,NormId,FeatDec);
elseif  idcmp(MethodId,'Isotropic'),
  % Coggins' ring and wedge shaped filters
  N_features=cp_isotropic(x,FeatureFN,NonlinId,SmoothId,NormId,FeatDec);
elseif  idcmp(MethodId,'Gabor:'),
  % Alternatively decomposed Gabor filters
  SepPos=findstr(MethodId,':');
  Structure=str2num(MethodId(SepPos(1)+1:length(MethodId)));
  N_features=cp_gabor2(x,FeatureFN,NonlinId,SmoothId,NormId,Structure,FeatDec);
elseif  idcmp(MethodId,'Gabor'),
  % Farrokhnias Gabor filters
  N_features=cp_gabor(x,FeatureFN,NonlinId,SmoothId,NormId,FeatDec);
elseif  idcmp(MethodId,'Eigen'),
  % Ades eigenfilter
  N_features=cp_eigen(x,FeatureFN,ImageId,NonlinId,SmoothId,NormId,FeatDec);
elseif  idcmp(MethodId,'Bovik'),
  NumFiltPrTex=2;
  N_features=cp_bovik(x,FeatureFN,ImageId,NonlinId,SmoothId,NormId,...
    FeatDec,N_tex,NumFiltPrTex);
elseif  idcmp(MethodId,'Weldon'),
  % Weldons multi-texture optimal Gabor filter design method
  SepPos=findstr(MethodId,':');
  NumFilt=str2num(MethodId(SepPos(1)+1:length(MethodId)));
  [N_features,p1]=cp_weldon(x,FeatureFN,ImageId,NonlinId,SmoothId,...
    NormId,FeatDec,N_tex,NumFilt,p1);
elseif  idcmp(MethodId,'Dunn'),
  % Dunn/Weldonss two-texture optimal Gabor filter design method
  SepPos=findstr(MethodId,':');
  sigma=str2num(MethodId(SepPos(1)+1:length(MethodId)));
  N_features=cp_dunn(x,FeatureFN,ImageId,NonlinId,SmoothId,...
    NormId,FeatDec,sigma);
elseif idcmp(MethodId,'Singh') | idcmp(MethodId,'Unser') 
  % Optimal filter bank using the Singh or Unser criterion
  SepPos=findstr(MethodId,':');
  if length(SepPos>=2)
    sigma=str2num(MethodId(SepPos(2)+1:length(MethodId)));
    MethodId=MethodId(1:SepPos(2)-1);
  else
    sigma=[];
  end
  if length(SepPos>=1)
    MN=str2num(MethodId(SepPos(1)+1:length(MethodId)));
    MN=[MN MN];
    MethodId=MethodId(1:SepPos(1)-1);
  else 
    MN=[7 7];
  end
  [N_features,ClusterA,ClusterB,p1]=cp_moptxf(x,FeatureFN,MethodId,...
    ImageId,NonlinId,SmoothId,NormId,FeatDec,N_tex,ClusterA,ClusterB,...
    MN,sigma,p1);
elseif idcmp(MethodId,'Fisher')
  % Optimal filter methods
  N_features=cp_optxf(x,FeatureFN,MethodId,...
    ImageId,NonlinId,SmoothId,NormId,FeatDec);
elseif  idcmp(MethodId,'NNFilt')
  if length(MethodId)==length('NNFilt')
    MaskSize=11;
  else 
    MaskSize=str2num(MethodId(length('NNFilt')+2:length(MethodId))); 
  end
  if idcmp(ImageId,'Nat-16')
    Nodes=[30 30];
  else
    Nodes=[20 20];
  end
  N_features=cp_nnfilt(x,FeatureFN,ImageId,NonlinId,...
    SmoothId,NormId,FeatDec,N_tex,MaskSize,Nodes);
elseif  idcmp(MethodId,'LinPred')
  if length(MethodId)==length('LinPred')
    r=4;
  else 
    r=str2num(MethodId(length('LinPred')+2:length(MethodId))); 
  end
  [N_features,p1]=cp_linpred(x,FeatureFN,...
    ImageId,NonlinId,SmoothId,NormId,FeatDec,N_tex,r,8,p1);
else
  error(sprintf('Method "%s" not implemented',MethodId));
end



