function N_features=cp_bovik(x,FeatureFN,ImageId,NonlinId,SmoothId,NormId,FeatDec,N_tex,NumFiltPrTex)
% CP_BOVIK  Comparative texture measure study: Optimal Gabor filterbank feat.


% Set initial values
MaxNumFiltPrTex=6;
B=1.0;
alpha=1;
N_features=0;

% Load database of normalization factors (variables: Ids and NormFacts)
if exist('BovikCF.mat')==2
  load BovikCF.mat
else
  ImageIds=[];
  UVs=[];
end

[dummy,N_tex,dummy,Components]=cp_image(ImageId,0);
for TexNo=1:N_tex
  % Select/retrieve filter center frequencies
  CF_idx=find(ImageIds==Components(TexNo));
  if isempty(CF_idx)
    x_i=cp_image(ImageId,TexNo);
    [M_x,N_x]=size(x_i);
    X_i=sqrt(abs(fftshift(fft2(x_i))));
    gimage(X_i,'colorbar','off');
    title(sprintf(...
      'Pick %d highest peaks for image D%d in descending order',...
      MaxNumFiltPrTex,Components(TexNo)));
    hold on; 
    plot([1/2 1/2]*M_x,[  0   1]*N_x,'c');
    plot([  0   1]*M_x,[1/2 1/2]*N_x,'c');
    hold off;
    [V,U]=ginput(MaxNumFiltPrTex);
    hold on; plot(V,U,'oy'); hold off; drawnow; pause(1);
    U=(U-1)/M_x - 0.5;
    V=(V-1)/N_x - 0.5;
    UV=U'+sqrt(-1)*V';
    fprintf(1,'Selected %d filters for texture D%d\n',...
      length(UV),Components(TexNo));
    ImageIds=[ImageIds; Components(TexNo)];
    UVs=[UVs; UV];
    save BovikCF.mat ImageIds UVs
    NewCFs=1;
  else
    UV=UVs(CF_idx,:);
    NewCFs=0;
  end
  
  % Design filter and extract features
  fprintf(1,'\t(D%d:',Components(TexNo));
  for fn=1:NumFiltPrTex
    N_features=N_features+1;
    CtrFreq=abs(UV(fn));
    Angle=angle(UV(fn));
    fprintf(1,' %4.2f/%4.0f',CtrFreq,Angle*180/pi);
    SigmaX=alpha/(pi*CtrFreq) * (2^B+1)/(2^B-1);
    SigmaY=20*atan((2^B-1)/(2^B+1));
    [y,h]=sepgabor(x,CtrFreq,Angle,SigmaX,SigmaY);
    if NewCFs
      % H=abs(fftshift(fft2(zeropad(h,M_x,N_x))));
      % gimage(X_i .* H,'colorbar','off'); colormap('jet'); 
      % title(sprintf('Filtering with filter #%d optimized for D%d',...
      %   fn,Components(TexNo)));
      % drawnow; pause(1);
    end
    v=cp_energy(y,NonlinId,SmoothId,NormId,[1 1]*CtrFreq); % Local energy func.
    catframe(v,FeatureFN,[],'Interm',FeatDec); % Write feature file
  end
  fprintf(1,')');
  if NewCFs
    clf;
  end
end




