function c=cp_main(ImageId,MethodId,NonlinId,SmoothId,NormId,ClassId,Conditional,SepTrainFract)
% CP_MAIN  Texture method comparison main function
%	c=cp_main(ImageId,MethodId,NonlinId,SmoothId,NormId,ClassId)
%	segments the specified image by the specified method. Only one
%	classifier is currently supported: ClassId='LVQ'. 
%	
%	Note: not all possible parameter combinations are meaningful!
%	
%	ImageId:
%	- 'Nat-5c'	Complex composite image Nat-5c
%	- 'Nat-16c'	Complex composite image Nat-16c
%	- [4 84]	Brodatz texture pair (example)
%	- [MSU_GMRF1 MSU_GMRF2]	See "help msu_texture_id" for a list of options
%	
%	NonlinId:	SmoothId:	NormId:		ClassId
%	- ''		- ''		- ''		- ''
%	- 'tanh'	- 'Gauss'	- 'Sqrt'	- 'LVQ'
%	- 'Sqr'  	- 'Average'	- 'Log'		- 'Thresh'
%	- 'Abs'				
%	
%	MethodId:
%	- 'cooc##'	Co-occurrence (Haralick et al., 1973) with 8 gray
%			levels and ## by ## windows.
%	- 'cooc'	Identical to cooc16
%	- 'SAR'		Simultaneous auto regressive (Mao & Jain, 1992)
%	- 'Laws'	Laws masks (Laws, 1980)
%	- 'Isotropic'	Isotropic filters (Coggins & Jain, 1985)
%	- 'Gabor'	Dyadic Gabor filter bank (Jain & Farrokhnia, 1991)
%	- 'QMF:$$:##:$'	Quadrature mirror filter (Randen & Husøy, 1994)
%			Wavelet transform (Mallat, 1989)
%			Wavelet packet (Chang & Kuo, 1993)
%			Wavelet frames (Unser, 1995)
%			Discrete cosine transform (Ng et al., 1992)
%	- 'Eigen'	Eigenfilter (Ade, 1983)
%	- 'Bovik'	Optimal Gabor filter bank (Bovik, 1990)
%	- 'Dunn'	Optimal Gabor filter (Dunn & Higgins, 1995)
%	- 'Weldon'	Optimal Gabor filter bank (Weldon & Higgins, 1996)
%	- 'Singh'	Singh correlation filters (Mahalanobis & Singh, 1994)
%	- 'Unser'	Unser-criterion optimal filters (Randen & Husøy)
%	- 'Fisher'	Fisher-criterion optimal filter (Randen & Husøy)
%	- 'NNFilt'	Back propagation method (Jain & Karu, 1996)
%	
%	Examples:
%	cp_main('Nat-5c','Laws','Sqr','Gauss','Log','LVQ');
%	cp_main('Nat-5c','Laws','Sqr','Gauss','Log','LVQ',...
%	  'Conditional','NonSepTrain');

% Necessary initializations to run as script:
% ImageId='GMRF';MethodId='Weldon:4';NonlinId='Sqr';SmoothId='Gauss';
% NormId='Log';ClassId='LVQ';

% --------------------------------------------------------------------
% INITIALIZATION
% --------------------------------------------------------------------

if ~exist('Conditional')
  Conditional=0;
else
  if idcmp(Conditional,'Conditional')
    Conditional=1;
  else
    Conditional=0;
  end
end

if ~exist('SepTrainFract')
  % Train on separate training feature vectors
  SepTrainFract=1;
else
  if isstr(SepTrainFract)
    if idcmp(SepTrainFract,'NonSepTrain')
      SepTrainFract=0;
    else
      SepTrainFract=1;
    end
  end
end

fprintf(1,'Removing files:\n');
TmpDir=tempnam; TmpDir=TmpDir(1:findstr(TmpDir,'.tmp.')-2);
eval(sprintf('!ls -lc %s/.tmp.*',TmpDir));
eval(sprintf('!rm -f %s/.tmp.*',TmpDir));

TrainDec=4;
TrainEdgeCut=32;
TestEdgeCut=20;
NumCBookVect=800;
RunLength=20000;

TmpFirstName=[]; % 			Empty => temporary files local
TmpFN=tempnam(TmpFirstName); % 		Temporary file
Tmp2FN=tempnam(TmpFirstName); % 	Temporary file
FeatureFN=tempnam(TmpFirstName); % 	Feature file name
TrFeatureFN=tempnam(TmpFirstName); % 	Train Feature file name
ReferenceFN=tempnam(TmpFirstName); % 	Reference file name
TrReferenceFN=tempnam(TmpFirstName); % 	Training reference file name
ClassMapFN=tempnam(TmpFirstName); % 	Classification result file name

% Result identifier (first part of file name for result files)
if isstr(ImageId)
  ResultId=ImageId;
else
  ResultId=sprintf('D%04d',ImageId);
end
ResultId=sprintf('%s-%s-%s-%s-%s-%s-%03d',...
  ResultId,MethodId,NonlinId,SmoothId,NormId,ClassId,...
  round(100*SepTrainFract));
if Conditional & (exist(['results/' ResultId '.gif'])==2)
  fprintf(1,'+----------------------------------------------------------------------\n');
  fprintf(1,'EXPERIMENT: %s EXISTING RESULTS, EXITING\n',ResultId);
  fprintf(1,'+----------------------------------------------------------------------\n');
  return;
end


fprintf(1,'+----------------------------------------------------------------------\n');
fprintf(1,'EXPERIMENT: %s\n',ResultId);
fprintf(1,'Feature file:  %s\n',FeatureFN);
fprintf(1,'Training file: %s\n',TrFeatureFN);
fprintf(1,'+----------------------------------------------------------------------\n');

[x,N_tex]=cp_image(ImageId); % 		Load image
[M_x,N_x]=size(x);


% --------------------------------------------------------------------
% FEATURE EXTRACTION (AND IN SOME CASES ALSO SEGMENTATION)
% Write output to the HIPS sequence FeatureFN
% --------------------------------------------------------------------
fprintf(1,'F E A T U R E   E X T R A C T I O N -----------------------\n');
[N_features,ClusterA,ClusterB,p1]=cp_feature(x,FeatureFN,...
  ImageId,MethodId,NonlinId,SmoothId,NormId,1,N_tex);
cp_normalize(FeatureFN,N_features,ImageId,MethodId,NonlinId,SmoothId,NormId);
if strcmp(ResultId,'Nat-5c-Gabor-Abs-Gauss--LVQ')
  eval(sprintf('!cp %s features.hips',FeatureFN));
end

% --------------------------------------------------------------------
% TRAINING DATA FEATURE EXTRACTION
% Write output to the HIPS sequence TrFeatureFN
% --------------------------------------------------------------------
TrReference=[];
if ~idcmp(MethodId,'NNFilt')
  if SepTrainFract>0
    for n=1:N_tex
      x_i=cp_image(ImageId,n);
      [M_xi,N_xi]=size(x_i);
      % Extract features
      cp_feature(x_i,TmpFN,ImageId,MethodId,NonlinId,...
	SmoothId,NormId,TrainDec,N_tex,ClusterA,ClusterB,p1);
      M_vi=hsize(sprintf('%s_000',TmpFN),'Rows');
      N_vi=hsize(sprintf('%s_000',TmpFN),'Cols');
      % Variance normalization
      cp_normalize(TmpFN,N_features,ImageId,MethodId,NonlinId,SmoothId,NormId);
      % Cut borders 
      M_vi=M_vi-2*TrainEdgeCut/TrainDec;
      N_vi=N_vi-2*TrainEdgeCut/TrainDec;
      eval(sprintf(...
	'!extract -s %d %d -p %d %d %s > %s',...
	M_vi,N_vi,[1 1]*TrainEdgeCut/TrainDec-1,TmpFN,Tmp2FN));
      % Concatenate 
      if n==1
	eval(sprintf('!mv %s %s',Tmp2FN,TrFeatureFN));
      else
	eval(sprintf('!%s %s %s | %s -s %d %d -p 0 0 > %s',...
	  'joinframes',TrFeatureFN,Tmp2FN,...
	  'extract',M_vi,n*N_vi,TmpFN));
	eval(sprintf('!mv %s %s',TmpFN,TrFeatureFN));
      end
      TrReference=[TrReference n*ones(M_vi,N_vi)];
    end
    if strcmp(ResultId,'Nat-5c-Gabor-Abs-Gauss--LVQ')
      eval(sprintf('!cp %s trfeatures.hips',TrFeatureFN));
    end
  else % SepTrainFract<=0
    eval(sprintf('!decimate -N 4 %s > %s',FeatureFN,TrFeatureFN));
    TrReference=cp_image(ImageId,0);
    TrReference=TrReference(1:4:end,1:4:end);
    SepTrainFract=1;
  end % if SepTrainFract
end

% --------------------------------------------------------------------
% CLASSIFICATION
% --------------------------------------------------------------------
fprintf(1,'C L A S S I F I C A T I O N -------------------------------\n');
Reference=cp_image(ImageId,0); % 	Get reference image
hipswrite(Reference,TmpFN);   eval(sprintf('!htob %s>%s',TmpFN,ReferenceFN));
hipswrite(TrReference,TmpFN); eval(sprintf('!htob %s>%s',TmpFN,TrReferenceFN));
ClassMap=cp_classify(FeatureFN,ReferenceFN,...
  TrFeatureFN,TrReferenceFN,ClassId,NumCBookVect,RunLength,...
  hsize(FeatureFN,'Rows'),hsize(FeatureFN,'Cols'),N_tex,...
  ClusterA,ClusterB,SepTrainFract);
Decimation=M_x/hsize(FeatureFN,'Rows');
% Decimate reference similar to the feature image
[M_R,N_R]=size(Reference);
Reference=Reference(1:Decimation:M_R,1:Decimation:N_R);
% Cut edges
MIdx=TestEdgeCut/Decimation+1:M_x/Decimation-TestEdgeCut/Decimation;
NIdx=TestEdgeCut/Decimation+1:N_x/Decimation-TestEdgeCut/Decimation;
Reference=Reference(MIdx,NIdx);
ClassMap=ClassMap(MIdx,NIdx);
[M_C,N_C]=size(ClassMap);


% --------------------------------------------------------------------
% SAVE AND ASSESS RESULTS
% --------------------------------------------------------------------
% Save class map
ColorMap=1:N_tex;
ColorMap=[rem(floor(ColorMap/2^0),2) 
  rem(floor(ColorMap/2^1),2) 
  rem(floor(ColorMap/2^2),2)]';
for n=9:8:N_tex
  idx=n:min(n+8-1,N_tex);
  ColorMap(idx,:)=(n-1)/N_tex*ColorMap(idx,:);
end
if min(ClassMap(:)) <= 0
  eval(sprintf('save results/%s.mat ClassMap',ResultId));
else
  gifwrite(ClassMap,ColorMap,['results/' ResultId '.gif']);
end
% Compute and save confusion matrix
Confusion=zeros(N_tex,N_tex);
for k=1:N_tex
  for l=1:N_tex
    Confusion(k,l)=100 * sum(sum((Reference==k) & (ClassMap==l))) ...
      / sum(sum(Reference==k));
  end
end
eval(sprintf('save -ascii results/%s-ConfMatr.txt Confusion',ResultId));
% Compute and save total classification error
Error=100 * sum(sum(Reference~=ClassMap)) / (M_C*N_C);
eval(sprintf('save -ascii results/%s-Error.txt Error',ResultId));
% Store number of features and clustering (if any)
eval(sprintf('save -ascii results/%s-NFeat.txt N_features',ResultId));
if ~isempty(ClusterA)
  Clustering=[ClusterA NaN*ones(N_features,1) ClusterB];
  eval(sprintf('save -ascii results/%s-Clust.txt Clustering',ResultId));
end




% --------------------------------------------------------------------
% CLEAN UP
% --------------------------------------------------------------------
eval(sprintf('!rm -f %s %s',TmpFN,Tmp2FN));
eval(sprintf('!rm -f %s %s',FeatureFN,TrFeatureFN));
eval(sprintf('!rm -f %s %s %s',ReferenceFN,TrReferenceFN,ClassMapFN));


fprintf(1,'FINISH: %s, %d features, error %5.2f%%\n',...
  ResultId,N_features,Error);

