mbindata2pca.cc

 const char *help = "\
 progname: mbindata2pca.cc\n\
 code2html: This program reads multiple bindata files and a PCA model, projects patterns into it.\n\
 version: Torch3 vision2.0, 2003-2005\n\
 (c) Sebastien Marcel (marcel@idiap.ch)\n";
 
 #include "string_utils.h"
 #include "ImageGray.h"
 #include "PCAMachine.h"
 #include "MyMeanVarNorm.h"
 #include "HistoEqualSmoothImagePreProcessing.h"
 #include "FileListCmdOption.h"
 #include "str_utils.h"
 #include "Random.h"
 #include "CmdLine.h"
 
 using namespace Torch;
 
 real DFS(int n_input, real *realinput, real *mean, int n_output, real *projection, real *eigenvalues, real *dffs_, real *difs_)
 {
 	real alpha = 0.5;
 	real alpha2 = alpha * alpha;
 	
 	real dffs = 0.0;
 	for(int i = 0 ; i < n_input ; i++)
 	{
 		real z = realinput[i] - mean[i];
 			
 		dffs += z * z;
 	}
 
 	real difs = 0.0;
 
 	for(int i = 0 ; i < n_output ; i++)
 	{
 		real z = projection[i] * projection[i];
 
 		dffs -= z;
 		difs += z / (eigenvalues[i] * alpha2);
 	}
 
 	dffs /= (alpha * eigenvalues[n_output]);
 
 	*dffs_ = dffs;
 	*difs_ = difs;
 
 	real dfs = difs + dffs;
 	
 	return dfs;
 }
 
 int main(int argc, char *argv[])
 {
 	char *dir_name;
 	char *model_filename;
 	int n_input;
 	real variance;
 	int n_output;
 	bool verbose;
 	int verbose_level;
 	bool normalize;
 	bool image_normalize;
 	int width, height;
 	bool noProjection;
 	bool doDFFS;
 	bool doDIFS;
 	int the_seed;
 	int n_perturbations;
 	int dim;
 	real nmin, nmax;
 
 	//
 	FileListCmdOption filelist("file name", "the list files or one data file");
         filelist.isArgument(true);
 
 	//
 	CmdLine cmd;
 	cmd.setBOption("write log", false);
 	
 	cmd.info(help);
 	cmd.addText("\nArguments:");
   	cmd.addCmdOption(&filelist);
 	cmd.addSCmdArg("model filename", &model_filename, "PCA model filename");
 	cmd.addICmdArg("n_input", &n_input, "number of inputs");
 	cmd.addText("\nOptions:");
 	cmd.addSCmdOption("-dir", &dir_name, ".", "dir name");
 	cmd.addRCmdOption("-variance", &variance, 0.95, "variance (-1 => 100\%)");
 	cmd.addICmdOption("-noutput", &n_output, -1, "number of outputs");
 	cmd.addBCmdOption("-verbose", &verbose, false, "verbose");
 	cmd.addICmdOption("-verbose_level", &verbose_level, 1, "level of verbose");
 	cmd.addBCmdOption("-normalize", &normalize, false, "normalize");
 	cmd.addBCmdOption("-imagenorm", &image_normalize, false, "considers the input pattern as an image and performs a photometric normalization");
 	cmd.addICmdOption("-width", &width, -1, "width");
 	cmd.addICmdOption("-height", &height, -1, "height");
 	cmd.addBCmdOption("-noProjection", &noProjection, false, "do not output projection");
 	cmd.addBCmdOption("-dffs", &doDFFS, false, "output DFFS");
 	cmd.addBCmdOption("-difs", &doDIFS, false, "output DIFS");
 	cmd.addICmdOption("-p", &n_perturbations, 0, "number of perturbations to generate");
 	cmd.addICmdOption("-dim", &dim, -1, "dimension to perturbate (all if -1)");
 	cmd.addICmdOption("-seed", &the_seed, 9503, "seed");
 	cmd.addRCmdOption("-nmin", &nmin, -0.1, "min random noise");
 	cmd.addRCmdOption("-nmax", &nmax, +0.1, "max random noise");
 	cmd.read(argc, argv);
 
 	if(verbose == false) verbose_level = 0;
 	
 	Random::manualSeed(the_seed);
 
 	if(image_normalize)
 	{
 		print("Perform photometric normalization ...\n");
 
 		print("The input pattern is an %dx%d image.\n", width, height);
 	}
 
 	if(n_perturbations < 0) n_perturbations = 0;
 
 	//
 	PCAMachine *pca_machine = NULL;
 	pca_machine = new PCAMachine(n_input);
 
 	//
 	PreProcessing *preprocess = NULL;
 	if(image_normalize)
 	{
 		print("Pre-processing image %dx%d ...\n", width, height);
 		preprocess = new HistoEqualSmoothImagePreProcessing(width, height);
 	}
 
 	//
 	MyMeanVarNorm *mv_norm = NULL;
 
 	if(normalize)
 	{
 	   	if(verbose_level >= 1)
 			print("Loading PCA model and its normalisation: %s ...\n", model_filename);
 
 		mv_norm = new MyMeanVarNorm(n_input, 1);
 
 	   	pca_machine->load(model_filename, mv_norm);
 	}
 	else
 	{
 	   	if(verbose_level >= 1)
 			print("Loading PCA model: %s ...\n", model_filename);
 
 		DiskXFile *file = NULL;
 		file = new DiskXFile(model_filename, "r");
 		pca_machine->loadXFile(file);
 		delete file;
 	}
 	   
 	//
 	pca_machine->setIOption("verbose_level", verbose_level);
 	pca_machine->setROption("variance", variance);
 	pca_machine->init();
 
 	if(variance <= 0.0 || n_output <= 0) n_output = pca_machine->n_outputs;
 	if(n_output > 0) pca_machine->n_outputs = n_output; 
 
 	//
 	if(verbose_level >= 1)
 		print("Projection bindata file into PCA space (%d -> %d) ...\n", n_input, n_output);
 
 	//
 	float *realinput = NULL;
         Sequence *seq;
 
         realinput = new float [n_input];
         seq = new Sequence(&realinput, 1, n_input);	
 
 	//
 	int dimIn;
 	int n_patterns;
 	
 	DiskXFile *pf = NULL;
 
 	//
 	for(int i = 0 ; i < filelist.n_files ; i++)
 	{
 		char *temp = strBaseName(filelist.file_names[i]);
 		char *file_name = strRemoveSuffix(temp);
 
 		if(verbose)
 			print("Processing file %s\n", file_name);
 
 		pf = new DiskXFile(filelist.file_names[i], "r");
 
 		pf->read(&n_patterns, sizeof(int), 1);
 		pf->read(&dimIn, sizeof(int), 1);
 
 		if(verbose_level >= 2)
 		{
 			print("n_inputs : %d\n", dimIn);
 			print("n_patterns : %d\n", n_patterns);  
 		}
 
 		if(n_input > dimIn)
 		{
 			error("Number of inputs specified (%d) bigger than into the file (%d)", n_input, dimIn);
 	   
 			delete pf;
 
 			return 0;
 		}
 
 
 
 		//
 		DiskXFile *pfOutput = NULL;
 		char filename_out[250];
 
 		sprintf(filename_out, "%s/%s.bindata", dir_name, file_name);
 		
 		pfOutput = new DiskXFile(filename_out, "w");
 	
 		int P = n_patterns + n_patterns*n_perturbations;
 	
 		pfOutput->write(&P, sizeof(int), 1);
 		int m = 0;
 		if(!noProjection) m = n_output;
 		if(doDFFS) m++;
 		if(doDIFS) m++;
 
 		if(m == 0)
 			error("Nothing to do then.");
 	
 		pfOutput->write(&m, sizeof(int), 1);
 		
 		real *machine_output = new real [n_output];
 
 		for(int p = 0 ; p < n_patterns ; p++)
 		{
 	   		//
 	  		for(int i = 0 ; i < n_input ; i++)
 				pf->read(&realinput[i], sizeof(float), 1);   
 
 	   		if(verbose_level >= 2)
 				print(" Seq =     [%2.3f %2.3f %2.3f ...]\n", realinput[0], realinput[1], realinput[2]);
 
 			if(image_normalize)
 				preprocess->preProcessInputs(seq);
 
 			if(normalize)
 				mv_norm->preProcessInputs(seq);
 
 			//
 			pca_machine->forward(seq);
 			
 	   		if(verbose_level >= 2)
 				print(" Output =   [%2.3f %2.3f %2.3f ...]\n", pca_machine->outputs->frames[0][0], pca_machine->outputs->frames[0][1], pca_machine->outputs->frames[0][2]);
 
 	  		for(int i = 0 ; i < n_output ; i++)
 				machine_output[i] = pca_machine->outputs->frames[0][i];
 
 
 			real dffs = 0.0;
 			real difs = 0.0;
 			real dfs = DFS(n_input, realinput, pca_machine->Xm, n_output, machine_output, pca_machine->eigenvalues, &dffs, &difs);
 
 			if(!noProjection)
 				pfOutput->write(machine_output, sizeof(float), n_output);   
 
 			if(doDFFS) pfOutput->write(&dffs, sizeof(float), 1);   
 			if(doDIFS) pfOutput->write(&difs, sizeof(float), 1);   
 
 			for(int j = 0 ; j < n_perturbations ; j++)
 			{
 				for(int i = 0 ; i < n_output ; i++)
 					if((dim == -1) || (i == dim)) machine_output[i] = pca_machine->outputs->frames[0][i] + Random::boundedUniform(nmin, nmax);
 					else machine_output[i] = pca_machine->outputs->frames[0][i];
 		
 				dffs = 0.0;
 				difs = 0.0;
 				dfs = DFS(n_input, realinput, pca_machine->Xm, n_output, machine_output, pca_machine->eigenvalues, &dffs, &difs);
 
 				if(!noProjection)
 					pfOutput->write(machine_output, sizeof(float), n_output);   
 
 				if(doDFFS) pfOutput->write(&dffs, sizeof(float), 1);   
 				if(doDIFS) pfOutput->write(&difs, sizeof(float), 1);   
 			}
 		}
 
 		//
 		delete pf;
 		delete pfOutput;
 		delete [] machine_output;
 	}
 
 	
 	//
 	delete seq;
 	delete [] realinput;
 
 	if(normalize) delete mv_norm;
 	if(image_normalize) delete preprocess;
 	delete pca_machine;
 	
   	return 0;
 }