bindata2pca.cc

  1 const char *help = "\
  2 progname: bindata2pca.cc\n\
  3 code2html: This program reads a bindata file and a PCA model, projects patterns into it.\n\
  4 version: Torch3 vision2.0, 2003-2005\n\
  5 (c) Sebastien Marcel (marcel@idiap.ch)\n";
  6 
  7 #include "ImageGray.h"
  8 #include "PCAMachine.h"
  9 #include "MyMeanVarNorm.h"
 10 #include "HistoEqualSmoothImagePreProcessing.h"
 11 #include "Random.h"
 12 #include "CmdLine.h"
 13 
 14 using namespace Torch;
 15 
 16 real DFS(int n_input, real *realinput, real *mean, int n_output, real *projection, real *eigenvalues, real *dffs_, real *difs_)
 17 {
 18 	real alpha = 0.5;
 19 	real alpha2 = alpha * alpha;
 20 	
 21 	real dffs = 0.0;
 22 	for(int i = 0 ; i < n_input ; i++)
 23 	{
 24 		real z = realinput[i] - mean[i];
 25 			
 26 		dffs += z * z;
 27 	}
 28 
 29 	real difs = 0.0;
 30 
 31 	for(int i = 0 ; i < n_output ; i++)
 32 	{
 33 		real z = projection[i] * projection[i];
 34 
 35 		dffs -= z;
 36 		difs += z / (eigenvalues[i] * alpha2);
 37 	}
 38 
 39 	dffs /= (alpha * eigenvalues[n_output]);
 40 
 41 	*dffs_ = dffs;
 42 	*difs_ = difs;
 43 
 44 	real dfs = difs + dffs;
 45 	
 46 	return dfs;
 47 }
 48 
 49 int main(int argc, char *argv[])
 50 {
 51 	char *bindata_filename;
 52 	char *model_filename;
 53 	int n_input;
 54 	char *filename_out;
 55 	real variance;
 56 	int n_output;
 57 	bool verbose;
 58 	int verbose_level;
 59 	bool normalize;
 60 	bool image_normalize;
 61 	int width, height;
 62 	bool noProjection;
 63 	bool doDFFS;
 64 	bool doDIFS;
 65 	int the_seed;
 66 	int n_perturbations;
 67 	int dim;
 68 	real nmin, nmax;
 69 	CmdLine cmd;
 70 	cmd.setBOption("write log", false);
 71 	
 72 	cmd.info(help);
 73 	cmd.addText("\nArguments:");
 74 	cmd.addSCmdArg("bindata filename", &bindata_filename, "bindata filename");
 75 	cmd.addSCmdArg("model filename", &model_filename, "PCA model filename");
 76 	cmd.addICmdArg("n_input", &n_input, "number of inputs");
 77 	cmd.addText("\nOptions:");
 78 	cmd.addSCmdOption("-o", &filename_out, "pca.bindata", "bindata output file");
 79 	cmd.addRCmdOption("-variance", &variance, 0.95, "variance (-1 => 100\%)");
 80 	cmd.addICmdOption("-noutput", &n_output, -1, "number of outputs");
 81 	cmd.addBCmdOption("-verbose", &verbose, false, "verbose");
 82 	cmd.addICmdOption("-verbose_level", &verbose_level, 1, "level of verbose");
 83 	cmd.addBCmdOption("-normalize", &normalize, false, "normalize");
 84 	cmd.addBCmdOption("-imagenorm", &image_normalize, false, "considers the input pattern as an image and performs a photometric normalization");
 85 	cmd.addICmdOption("-width", &width, -1, "width");
 86 	cmd.addICmdOption("-height", &height, -1, "height");
 87 	cmd.addBCmdOption("-noProjection", &noProjection, false, "do not output projection");
 88 	cmd.addBCmdOption("-dffs", &doDFFS, false, "output DFFS");
 89 	cmd.addBCmdOption("-difs", &doDIFS, false, "output DIFS");
 90 	cmd.addICmdOption("-p", &n_perturbations, 0, "number of perturbations to generate");
 91 	cmd.addICmdOption("-dim", &dim, -1, "dimension to perturbate (all if -1)");
 92 	cmd.addICmdOption("-seed", &the_seed, 9503, "seed");
 93 	cmd.addRCmdOption("-nmin", &nmin, -0.1, "min random noise");
 94 	cmd.addRCmdOption("-nmax", &nmax, +0.1, "max random noise");
 95 	cmd.read(argc, argv);
 96 
 97 	if(verbose == false) verbose_level = 0;
 98 	
 99 	Random::manualSeed(the_seed);
100 
101 	if(image_normalize)
102 	{
103 		print("Perform photometric normalization ...\n");
104 
105 		print("The input pattern is an %dx%d image.\n", width, height);
106 	}
107 
108 	if(n_perturbations < 0) n_perturbations = 0;
109 
110 	//
111 	PCAMachine *pca_machine = NULL;
112 	pca_machine = new PCAMachine(n_input);
113 
114 	//
115 	PreProcessing *preprocess = NULL;
116 	if(image_normalize)
117 	{
118 		print("Pre-processing image %dx%d ...\n", width, height);
119 		preprocess = new HistoEqualSmoothImagePreProcessing(width, height);
120 	}
121 
122 	//
123 	MyMeanVarNorm *mv_norm = NULL;
124 
125 	if(normalize)
126 	{
127 	   	if(verbose_level >= 1)
128 			print("Loading PCA model and its normalisation: %s ...\n", model_filename);
129 
130 		mv_norm = new MyMeanVarNorm(n_input, 1);
131 
132 	   	pca_machine->load(model_filename, mv_norm);
133 	}
134 	else
135 	{
136 	   	if(verbose_level >= 1)
137 			print("Loading PCA model: %s ...\n", model_filename);
138 
139 		DiskXFile *file = NULL;
140 		file = new DiskXFile(model_filename, "r");
141 		pca_machine->loadXFile(file);
142 		delete file;
143 	}
144 	   
145 	//
146 	pca_machine->setIOption("verbose_level", verbose_level);
147 	pca_machine->setROption("variance", variance);
148 	pca_machine->init();
149 
150 	if(variance <= 0.0 || n_output <= 0) n_output = pca_machine->n_outputs;
151 	if(n_output > 0) pca_machine->n_outputs = n_output; 
152 
153 	//
154 	int dimIn;
155 	int n_patterns;
156 
157 	DiskXFile *pf = NULL;
158 	pf = new DiskXFile(bindata_filename, "r");
159 
160 	pf->read(&n_patterns, sizeof(int), 1);
161 	pf->read(&dimIn, sizeof(int), 1);
162 
163 	if(verbose_level >= 1)
164 	{
165 		print("n_inputs : %d\n", dimIn);
166 		print("n_patterns : %d\n", n_patterns);  
167 	}
168 
169 	if(n_input > dimIn)
170 	{
171 		error("Number of inputs specified (%d) bigger than into the file (%d)", n_input, dimIn);
172 	   
173 		delete pf;
174 
175 		return 0;
176 	}
177 
178 
179 	//
180 	float *realinput = NULL;
181         Sequence *seq;
182 
183         realinput = new float [n_input];
184         seq = new Sequence(&realinput, 1, n_input);	
185 
186 	//
187 	DiskXFile *pfOutput = NULL;
188 	pfOutput = new DiskXFile(filename_out, "w");
189 	
190 	//
191 	if(verbose_level >= 1)
192 		print("Projection bindata file into PCA space (%d -> %d) ...\n", n_input, n_output);
193 
194 	int P = n_patterns + n_patterns*n_perturbations;
195 	
196 	pfOutput->write(&P, sizeof(int), 1);
197 	int m = 0;
198 	if(!noProjection) m = n_output;
199 	if(doDFFS) m++;
200 	if(doDIFS) m++;
201 
202 	if(m == 0)
203 		error("Nothing to do then.");
204 	
205 	pfOutput->write(&m, sizeof(int), 1);
206 		
207 	real *machine_output = new real [n_output];
208 
209 	for(int p = 0 ; p < n_patterns ; p++)
210 	{
211 	   	//
212 	  	for(int i = 0 ; i < n_input ; i++)
213 			pf->read(&realinput[i], sizeof(float), 1);   
214 
215 	   	if(verbose_level >= 2)
216 			print(" Seq =     [%2.3f %2.3f %2.3f ...]\n", realinput[0], realinput[1], realinput[2]);
217 
218 		if(image_normalize)
219 			preprocess->preProcessInputs(seq);
220 
221 		if(normalize)
222 			mv_norm->preProcessInputs(seq);
223 
224 		//
225 		pca_machine->forward(seq);
226 			
227 	   	if(verbose_level >= 2)
228 			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]);
229 
230 	  	for(int i = 0 ; i < n_output ; i++)
231 			machine_output[i] = pca_machine->outputs->frames[0][i];
232 
233 
234 		real dffs = 0.0;
235 		real difs = 0.0;
236 		real dfs = DFS(n_input, realinput, pca_machine->Xm, n_output, machine_output, pca_machine->eigenvalues, &dffs, &difs);
237 
238 		if(!noProjection)
239 			pfOutput->write(machine_output, sizeof(float), n_output);   
240 
241 		if(doDFFS) pfOutput->write(&dffs, sizeof(float), 1);   
242 		if(doDIFS) pfOutput->write(&difs, sizeof(float), 1);   
243 
244 		for(int j = 0 ; j < n_perturbations ; j++)
245 		{
246 			for(int i = 0 ; i < n_output ; i++)
247 				if((dim == -1) || (i == dim)) machine_output[i] = pca_machine->outputs->frames[0][i] + Random::boundedUniform(nmin, nmax);
248 				else machine_output[i] = pca_machine->outputs->frames[0][i];
249 		
250 			dffs = 0.0;
251 			difs = 0.0;
252 			dfs = DFS(n_input, realinput, pca_machine->Xm, n_output, machine_output, pca_machine->eigenvalues, &dffs, &difs);
253 
254 			if(!noProjection)
255 				pfOutput->write(machine_output, sizeof(float), n_output);   
256 
257 			if(doDFFS) pfOutput->write(&dffs, sizeof(float), 1);   
258 			if(doDIFS) pfOutput->write(&difs, sizeof(float), 1);   
259 		}
260 	}
261 
262 	//
263 	delete pf;
264 	delete pfOutput;
265 	
266 	//
267 	delete seq;
268 	delete [] realinput;
269 	if(normalize) delete mv_norm;
270 	if(image_normalize) delete preprocess;
271 	delete pca_machine;
272 	
273   	return 0;
274 }
275