function ClassMap=cp_classify(FeatureFN,ReferenceFN,TrFeatureFN,TrReferenceFN,ClassId,NumCBookVect,RunLength,M_x,N_x,N_tex,ClusterA,ClusterB,TrainingFraction)
% CP_CLASSIFY  Classify a feature map
%	ClassMap=cp_classify(FeatureFN,ReferenceFN,...
%	  TRFeatureFN,TrReferenceFN,ClassId,NumCBookVect,RunLength,...
%	  M_x,N_x,N_tex,ClusterA,ClusterB,TrainingFraction)

% Initialization
LVQ_path=sprintf('/usr/local/lvq/%s',getenv('HOSTTYPE'));
eveninit=['eveninit']; % 		[LVQ_path '/eveninit'];
balance =['balance']; % 		[LVQ_path '/balance'];
olvq1   =['olvq1']; % 			[LVQ_path '/olvq1'];
classify=['classify']; % 		[LVQ_path '/classify'];
Decimals=4;
CBookFN=tempnam;
LvqFN=tempnam;
TrLvqFN=tempnam;
ClLvqFN=tempnam;
ClFN=tempnam;
TmpFN=tempnam;
DummyFN=tempnam;
if isempty(ClusterA)
  ClusterA=2*ones(hsize(FeatureFN,'Frames'),1);
  ClusterB=1*ones(hsize(FeatureFN,'Frames'),1);
end
M_cl=size(ClusterA,1);
N_cl=max(size(ClusterA,2),size(ClusterB,2));
ClusterA=zeropad(ClusterA,M_cl,N_cl);
ClusterB=zeropad(ClusterB,M_cl,N_cl);

if     strcmp(upper(ClassId),'')
  ClassMap=hipsread(FeatureFN,[],[],0)+1;  
elseif strcmp(upper(ClassId),'LVQ')
  fprintf(1,'LVQ training and classification\n');
  fprintf(1,'    Training reference     %s\n',TrReferenceFN);
  fprintf(1,'    Training feature, HIPS %s\n',TrFeatureFN);
  fprintf(1,'    Training feature, LVQ  %s\n',TrLvqFN);
  fprintf(1,'    Feature file, HIPS     %s\n',FeatureFN);
  fprintf(1,'    Feature file, LVQ      %s\n',LvqFN);
  fprintf(1,'    Codebook file          %s\n',CBookFN);
  fprintf(1,'    Class map file, LVQ    %s\n',ClLvqFN);
  fprintf(1,'    Class map file, HIPS   %s\n',ClFN);

  % --------------------------------------------------------------------
  % Training
  % --------------------------------------------------------------------
  fprintf(1,'Training:\n');
  eval(sprintf('!hips2lvq -d %d -p %f -train %s -out %s %s %s',...
    Decimals,TrainingFraction*100,TrLvqFN,DummyFN,TrFeatureFN,TrReferenceFN)); % 	Convert to LVQ format
  eval(sprintf('!ls -l %s %s %s %s',TrLvqFN,DummyFN,TrFeatureFN,TrReferenceFN));
  eval(sprintf('!sed "s/\?/0/g" %s > %s; mv %s %s',...
      TrLvqFN,TmpFN,TmpFN,TrLvqFN));
  eval(sprintf('!ls -l %s',TrLvqFN));
  eval(sprintf('!%s -din %s -cout %s -noc %d',eveninit,...
    TrLvqFN,CBookFN,NumCBookVect));
  eval(sprintf('!ls -l %s',CBookFN));
  eval(sprintf('!%s -din %s -cin %s -cout %s',...
    balance,TrLvqFN,CBookFN,CBookFN));
  eval(sprintf('!ls -l %s %s',CBookFN,TrLvqFN));
  eval(sprintf('!%s -din %s -cin %s -cout %s -rlen %d',olvq1,...
    TrLvqFN,CBookFN,CBookFN,RunLength));

  % --------------------------------------------------------------------
  % Classification
  % --------------------------------------------------------------------
  fprintf(1,'Classification:\n');
  eval(sprintf('!hips2lvq -out %s %s',LvqFN,FeatureFN));
  eval(sprintf('!sed "s/\?/0/g" %s > %s; mv %s %s',LvqFN,TmpFN,TmpFN,LvqFN));
  eval(sprintf('!%s -din %s -cin %s -dout %s',classify,LvqFN,CBookFN,ClLvqFN));
  eval(sprintf('!lvq2hips -s %d %d %s | htob > %s',M_x,N_x,ClLvqFN,ClFN));
  
  ClassMap=hipsread(ClFN);
  eval(sprintf('!rm -f %s %s %s %s %s',...
    CBookFN,LvqFN,TrLvqFN,ClLvqFN,ClFN,DummyFN));
elseif strcmp(lower(ClassId),'thresh')
  ref=hipsread(ReferenceFN);
  ref_tr=hipsread(TrReferenceFN);
  ClassRepr=0;
  Thresh=[];
  OverThresh=[];
  for FeatureNo=1:size(ClusterA,1)
    % Extract information about the current cluster pair
    ThisClusterA=nonzero(ClusterA(FeatureNo,:));
    ThisClusterB=nonzero(ClusterB(FeatureNo,:));
    % Load feature values
    v   =hipsread(FeatureFN  ,[],[],FeatureNo-1);
    v_tr=hipsread(TrFeatureFN,[],[],FeatureNo-1);
    % Extract feature values for current cluster pair
    v_trA=[]; v_trB=[];
    for TexNo=[ThisClusterA ThisClusterB]
      if sum(ThisClusterA==TexNo)>0
	v_trA=[v_trA; v_tr(find(ref_tr==TexNo))];
      elseif sum(ThisClusterB==TexNo)>0
	v_trB=[v_trB; v_tr(find(ref_tr==TexNo))];
      else
	error('Illegal state');
      end
    end
    % Compute mean and standard deviation for features
    m_trA(FeatureNo)=mean(v_trA(:)); std_trA(FeatureNo)=std(v_trA(:));
    m_trB(FeatureNo)=mean(v_trB(:)); std_trB(FeatureNo)=std(v_trB(:));
    % Threshold => class map for this cluster pair
    [NewClassMap,Thresh(FeatureNo),OverThresh(FeatureNo)] ...
      =opthresh(v,v_trA,v_trB);
    fprintf(1,'Thresholding textures%s vs.%s using threshold %f\n',...
      sprintf(' %d',ThisClusterA),sprintf(' %d',ThisClusterB),...
      Thresh(FeatureNo));
    % Store the cluster-pair-class-maps in one image
    ClassRepr=ClassRepr + NewClassMap * 2^FeatureNo;
  end

  % Combine pairwise class maps into one overall class map
  ClassMap=0;
  for TexNo=1:N_tex
    ThisClass=1;
    % Find out what filters are separating texture TexNo
    for f=find(sum([ClusterA ClusterB]'==TexNo))
      fprintf(1,'.');
      if (sum(ClusterA(f,:)==TexNo) & OverThresh(f)==1) ...
          | (sum(ClusterB(f,:)==TexNo) & OverThresh(f)==2)
        % The output of the filter should be above the threshold
        ThisClass=ThisClass & rem(floor(ClassRepr/2^f),2);
      else
        ThisClass=ThisClass & ~rem(floor(ClassRepr/2^f),2);
      end
    end
    ClassMap = ClassMap + TexNo*ThisClass;
  end
  fprintf(1,'\n');
  
  % Handle doubt regions
  fprintf(1,'Doubt regions %4.2f%%',...
    100*sum(ClassMap(:)==0)/(size(ClassMap,1)*size(ClassMap,2)));
  %% while sum(ClassMap(:)==0)>0
  %%   for k=1:N_tex
  %%     ClassMap = ((ClassMap~=0) .* ClassMap) ...
  %% 	+ k * ((ClassMap==0) .* dilate(ClassMap==k,ones(3,3)));
  %%   end
  %% end
  % Compute sum of weighted distance from each threshold for each textures
  % and each nonclassified pixel
  Idx=find(ClassMap==0);
  if length(Idx)>0
    SumDist=zeros(N_tex,length(Idx));
    for FeatureNo=1:size(ClusterA,1)
      fprintf(1,'.');
      ThisClusterA=nonzero(ClusterA(FeatureNo,:));
      ThisClusterB=nonzero(ClusterB(FeatureNo,:));
      v=hipsread(FeatureFN  ,[],[],FeatureNo-1); % Load features
      % Compute weigthed Euclidean distance from nearest class
      for t=ThisClusterA
	SumDist(t,:)=SumDist(t,:)+abs((v(Idx)-m_trA(FeatureNo)).') ...
	  /(std_trA(FeatureNo)*sum(sum([ClusterA ClusterB]==t)));
      end
      for t=ThisClusterB
	SumDist(t,:)=SumDist(t,:)+abs((v(Idx)-m_trB(FeatureNo)).') ...
	  /(std_trB(FeatureNo)*sum(sum([ClusterA ClusterB]==t)));
      end
    end
    [dummy,ClosestClass]=min(SumDist);
    ClassMap(Idx)=ClosestClass;
  end
  fprintf(1,'\n');
elseif strcmp(lower(ClassId),'min')
  % Set the class equal to the feature frame number with minimum feature
  % value, corresponding to prediction error energy or similar features
  ClassMap=zeros(M_x,N_x);
  v_min   =Inf * ones(M_x,N_x);
  N_features=hsize(FeatureFN,'Frames');
  for FeatureNo=1:N_features
    % Load feature image
    v=hipsread(FeatureFN,[],[],FeatureNo-1);
    LowestBM=v<v_min;
    ClassMap=LowestBM*FeatureNo + (~LowestBM).*ClassMap;
    v_min=min(v,v_min);
  end
elseif strcmp(lower(ClassId),'max')
  % Set the class equal to the feature frame number with maximum feature
  % value, corresponding to prediction error energy or similar features
  ClassMap=zeros(M_x,N_x);
  v_max   =-Inf * ones(M_x,N_x);
  N_features=hsize(FeatureFN,'Frames');
  for FeatureNo=0:N_features-1
    % Load feature image
    v=hipsread(FeatureFN,[],[],FeatureNo-1);
    HighestBM=v>v_max;
    sum(HighestBM(:))
    ClassMap=HighestBM*(floor(FeatureNo*N_tex/N_features)+1) ...
	+ (~HighestBM).*ClassMap;
    v_max=max(v,v_max);
  end
else
  error(sprintf('Unsupported classification method "%s"',ClassId));
end
