facenet.cpp

Gaussian Mixture Algorithm

#include "facenet.h"

#include <stdio.h>
#include <string.h>
#include <math.h>

using namespace std;
using namespace ebl;


//////////////////////////////////////////////////////
///
//////////////////////////////////////////////////////

Idx<intg> custom_table(int in, int out){
	int out2 = out-2;
	Idx<intg> tbl((3.5 * out2)+(2*in), 2);
    for(int i = 0; i < out2/2; i++){
    	for(int j = 0; j <= 2; j++){
    		tbl.set(i, (3*i)+j, 1);
    		tbl.set((i+j)% in, 3*i +j, 0);
    	}
    }
    for(int i = out2/2; i < out2; i++){
    	for(int j = 0; j <= 3; j++){
    		tbl.set(i, i*4 - out2/2 + j, 1);
    		tbl.set((i+j)% in,  i*4 - out2/2 + j, 0);
    	}
    }
    for(int j = 0; j < in; j++){
    	tbl.set(j, 3.5*out2 + j, 0);
    	tbl.set(out2,  3.5*out2 + j, 1);
    }
    for(int j = 0; j < in; j++){
    	tbl.set(j, 3.5*out2 + j + in, 0);
    	tbl.set(out-1,  3.5*out2 + j + in, 1);
    }
    return tbl;
}


Idx<double> mexican_hat(double s, int n){
	Idx<double> m(n, n);
	double vinv = 1/(s*s);
	double total = 0;
	int cx = n/2;
	int cy = n/2;
	for(int x = 0; x < n; x++){
		for(int y = 0; y < n; y++){
			int dx = x - cx;
			int dy = y - cy;
			m.set( -exp(-sqrt(vinv*(dx*dx + dy*dy))) , x, y);
			total += m.get(x, y);
		}
	}
	//! set center valus so it's zero sum
	m.set( m.get(cx, cy) - total, cx, cy);
	//! normalize so that energy is 1
	double energy = sqrt(idx_sumsqr(m));
	idx_dotc(m, 1/energy, m);
	return m;
}


//////////////////////////////////////////////////////////////////

facenet_datasource::facenet_datasource(Idx<double> *inputs, Idx<ubyte> *labels, double m, double d, bool pretreatment) :
	LabeledDataSource<double, ubyte>(inputs, labels)
	{
		height = data->dim(1);
		width = data->dim(2);
		mymean = m;
		mydeviation = d;
		if(pretreatment){
			mexican_filter(2, 5);
			mymean = idx3_mean(data);
			mydeviation = idx3_coef(data, mymean);
			normalize();
		}
	}


void facenet_datasource::mexican_filter(double s, int n){
	Idx<double> filter = mexican_hat(s, n);
	for(int layer = 0; layer < data->dim(0); layer++){
		Idx<double> big(height + floor(n/2), width + floor(n/2));
		idx_fill(big, 0.0);
		Idx<double> bla = data->select(0, layer);
		Idx<double> bla2 = big.narrow(0, height, floor(n/2)).narrow(1, width, floor(n/2));
		idx_copy( bla, bla2);
		idx_2dconvol(big, filter, bla);
	}
	printf("High-pass filtering performed \n");
}

void facenet_datasource::save(string dir, string name){
	string s1 = dir + "/" + name + "_inputs.mat";
	save_matrix(*data, s1.c_str());
	string s2 = dir + "/" + name + "_labels.mat";
	save_matrix(*labels, s2.c_str());
}

void facenet_datasource::normalize(){
	for(int layer = 0; layer < data->dim(0); layer++){
		Idx<double> bla = data->select(0, layer);
		Idx<double> bla2(height, width);
		idx_addc(bla, -mymean, bla2);
		idx_dotc(bla2, 1/mydeviation, bla);
	}
	printf("Normalization performed \n");
}


/////////////////////////////////////////////////////////////////

facenet::facenet(const char *paramfile){
	theparam = new parameter(60000);

	int nclasses = 2;
	//! the target values for mean-squared error training
	//! are +target for the output corresponding to the class
	//! being shown and -target for the other outputs.
	double target = 1;
	labels = new Idx<ubyte>(nclasses);
	targets = new Idx<double>(nclasses, nclasses);
	idx_fill( *targets, -target);
	for(int i = 0; i < nclasses; i++){
	     targets->set(target, i, i);
	     labels->set(i, i);
	}
	idx_dotc(*targets, 1.5, *targets);

	table0 = full_table(1, 6);
	table1 = custom_table(6, 16);
	table2 = custom_table(16, 80);

	//! initializes the net (cscscf) and charges the weights if needed
	init(theparam, 42, 42, 7, 7, &table0, 2, 2, 7, 7, &table1, 2, 2, 6, 6, &table2, 2);
	if((paramfile != NULL) && (theparam->load(paramfile))){
		printf("Weights loaded. \n");
	}
	else{
		printf("Weights loaded \n");
	}

	//! creates the kernels for smoothing and high-pass filtering
	double sker[3][3] = {{0.3, 0.5, 0.3}, {0.5, 1, 0.5}, {0.3, 0.5, 0.3}};
	smoothing_kernel = Idx<double>(3, 3);
	memcpy(smoothing_kernel.idx_ptr(), sker, sizeof (sker));

	highpass_kernel = mexican_hat(2, 5);
}

facenet::~facenet() {
	delete theparam;
	{ idx_bloop3(in, inputs, void*, out, outputs, void*, r, results, void*) {
		delete((state_idx*) in.get());
		delete((state_idx*) out.get());
		delete((Idx<double>*) r.get());
	}}
	delete labels;
	delete targets;

}


void facenet::build_databases( int sz, const string directory, float ratio){
	Idx<double> dataface(1,1,1), datanonface(1,1,1);
	Idx<ubyte> lblface(1), lblnonface(1);

	//! loads the various images idxs
	string s1 = directory + "/data_face.mat";
	load_matrix( dataface, s1.c_str());
	string s2 = directory + "/labels_face.mat";
	load_matrix( lblface, s2.c_str());
	string s3 = directory + "/data_nonface.mat";
	load_matrix( datanonface, s3.c_str());
	string s4 = directory + "/labels_nonface.mat";
	load_matrix( lblnonface, s4.c_str());

	//! initializes the various sets of images and according labels
	Idx<double> training_data( dataface.dim(0)+ datanonface.dim(0), sz, sz);
	Idx<double> validation_data( dataface.dim(0)+ datanonface.dim(0), sz, sz);
	Idx<ubyte> training_lbl( lblface.dim(0)+ lblnonface.dim(0));
	Idx<ubyte> validation_lbl( lblface.dim(0)+ lblnonface.dim(0));
	int n_train = 0, n_valid = 0;
	int imax = min(dataface.dim(0)-1, datanonface.dim(0) - 1);

	//! distributes randomly the images into the training and validation datasets, according to the ratio
	dseed(0);
	for(int i = 0; i <= imax; i++){
		if(drand() < ratio){
			Idx<double> bla = datanonface.select(0,i);
			Idx<double> bla2 = training_data.select(0, n_train);
			idx_copy( bla, bla2);
			training_lbl.set( lblnonface.get(i) , n_train);
			n_train++;
			bla = dataface.select(0,i);
			bla2 = training_data.select(0, n_train);
			idx_copy( bla, bla2);
			training_lbl.set( lblface.get(i) , n_train);
			n_train++;
		} else {
			Idx<double> bla = datanonface.select(0,i);
			Idx<double> bla2 = validation_data.select(0, n_valid);
			idx_copy( bla, bla2);
			validation_lbl.set( lblnonface.get(i) , n_valid);
			n_valid++;
			bla = dataface.select(0,i);
			bla2 = validation_data.select(0, n_valid);
			idx_copy( bla, bla2);
			validation_lbl.set( lblface.get(i) , n_valid);
			n_valid++;
		}
	}

	training_data = training_data.narrow(0, n_train, 0);
	training_lbl = training_lbl.narrow(0, n_train, 0);
	validation_data = validation_data.narrow(0, n_valid, 0);
	validation_lbl = validation_lbl.narrow(0, n_valid, 0);

	//! saves the idxs for future use
	facenet_datasource* trainingdb = new facenet_datasource( &training_data, &training_lbl, 0, 0, true);
	facenet_datasource* testingdb = new facenet_datasource( &validation_data, &validation_lbl, 0, 0, true);

	trainingdb->save(directory, "training");
	testingdb->save(directory, "testing");

	printf("Databases built \n");

}


void facenet::train(const string directory){
	//! create the network
	idx3_supervised_module* thenet = new idx3_supervised_module(this, new edist_cost(labels, 1, 1, targets), new max_classer(labels));
	//! create the trainer
	supervised_gradient* thetrainer = new supervised_gradient( thenet, theparam);
	//! a classifier-meter measures classification errors
	classifier_meter* trainmeter = new classifier_meter();
	classifier_meter* testmeter = new classifier_meter();

	//! initialize the network weights
	forget_param_linear forgetparam(1, 1/2);
	init_drand(0);
	thenet->machine->forget(forgetparam);

	//! load the training and testing databases
	Idx<double> training_input(1,1,1);
	Idx<ubyte> training_label(1);
	Idx<double> testing_input(1,1,1);
	Idx<ubyte> testing_label(1);

	string s1 = directory + "/training_inputs.mat";
	load_matrix(training_input, s1.c_str());
	string s2 = directory + "/training_labels.mat";
	load_matrix(training_label, s2.c_str());
	string s3 = directory + "/testing_inputs.mat";
	load_matrix(testing_input, s3.c_str());
	string s4 = directory + "/testing_labels.mat";
	load_matrix(testing_label, s4.c_str());

	facenet_datasource* trainingdb = new facenet_datasource(&training_input, &training_label);
	facenet_datasource* testingdb = new facenet_datasource( &testing_input, &testing_label);
	printf("Databases loaded \n");

	//! do training iterations (here only one)
	printf("training with %d training samples and %d test samples\n", trainingdb->size(), testingdb->size());
	gd_param* gdp = new gd_param(0.0001, 0, 0, 0, 0, 0, 0, 0, 0);

	//! perform only one pass to check that the training is working
	train_online(thetrainer, trainingdb, testingdb, trainmeter, testmeter, 1, gdp, 0, directory);

	//! basic training scheme
	//train_online(thetrainer, trainingdb, testingdb, trainmeter, testmeter, 3, gdp, 0, directory);
	//gdp->eta = 0.00003;
	//train_online(thetrainer, trainingdb, testingdb, trainmeter, testmeter, 2, gdp, 1, directory);
	//gdp->eta = 0.00001;
	//train_online(thetrainer, trainingdb, testingdb, trainmeter, testmeter, 2, gdp, 2, directory);


}


void facenet::train_online(supervised_gradient* trainer, facenet_datasource* trainingdb, facenet_datasource* testingdb, classifier_meter* trainmeter, classifier_meter* testmeter, int n, gd_param* eta, int pass, const string directory){
	int i = 0;
	for(int j = 0; j < n; j++){
		time_t t = clock();
	    //! estimate second derivative on 100 iterations, using mu=0.02 and set individual espilons
	    printf("computing diagonal hessian and learning rates\n");
	    trainer->compute_diaghessian(trainingdb, 100, 0.02);
	    trainer->train(trainingdb, trainmeter, eta, 0);
	    time_t t2 = clock() - t;
	    printf("time : %li \n", (int)t2/CLOCKS_PER_SEC);
	    printf("training: ");
	    trainer->test(trainingdb, trainmeter);
	    trainmeter->display();
	    printf(" testing: ");
	    trainer->test(testingdb, testmeter);
	    testmeter->display();
		//! save the generated weights
	    std::stringstream ss;
	    ss << pass << i;
	    string s = directory + "/weights" + ss.str() + ".mat";
	    save_matrix( theparam->x, s.c_str());
	    i++;
	}
}



Idx<double> facenet::detect(Idx<ubyte> *img, int h, int w, Idx<int> &sz, float zoom, double threshold, int objsize) {
	height = h;
	width = w;
	grabbed = Idx<ubyte>(height, width);

	sizes = sz;
	//! initialize input and output states and result matrices for each size
	inputs = Idx<void*>(sizes.nelements());
	outputs = Idx<void*>(sizes.nelements());
	results = Idx<void*>(sizes.nelements());
	{ idx_bloop4(size, sizes, int, in, inputs, void*, out, outputs, void*, r, results, void*) {
		in.set((void*) new state_idx(2, 38 + (size.get() * 4), 38 + (size.get() * 4)));
		out.set((void*) new state_idx(labels->nelements(), size.get(), size.get()));
		r.set((void*) new Idx<double>(size.get(), size.get(), 2)); // (class,score)
	}}

	//! prepares the multiple scales
	Idx<ubyte> display = this->multi_res_prep(img, 0.5);
	//! performs the fprop
	Idx<double> res = this->multi_res_fprop(threshold, objsize);
	//! draws the rectangles around the hot spots detected
	{ idx_bloop1(re, res, double) {
		image_draw_box(display, (ubyte)255,
				(unsigned int) (zoom * (re.get(3) - (0.5 * re.get(5)))),
				(unsigned int) (zoom * (re.get(2) - (0.5 * re.get(4)))),
				(unsigned int) (zoom * re.get(5)),
				(unsigned int) (zoom * re.get(4)));
	}}

	idx_copy(grabbed, *img);
	return res;
}




void facenet::mark_maxima(Idx<double> &in, Idx<double> &inc, Idx<double> &res, double threshold) {
	idx_clear(res);
	int tr[] = { 1, 2, 0 };
	Idx<double> s(inc.transpose(tr));
	Idx<double> z(in.transpose(tr));
	z = z.unfold(0, 3, 1);
	z = z.unfold(1, 3, 1);
	{ idx_bloop3(se, s, double, ze, z, double, re, res, double) {
		{ idx_bloop3(see, se, double, zee, ze, double, ree, re, double)  {
			// find winning class
			intg w = 1;
			double c = see.get(w);
			// look if above threshold and local maximum
			ree.set(-1.0, 0),
			ree.set(-100.0, 1);
			if ((c > threshold) &&
					(c > zee.get(w, 0, 0)) && (c > zee.get(w, 0, 1)) && (c > zee.get(w, 0, 2)) &&
					(c > zee.get(w, 1, 0)) && (c > zee.get(w, 1, 2)) &&
					(c > zee.get(w, 2, 0)) && (c > zee.get(w, 2, 1)) && (c > zee.get(w, 2, 2))) {
				ree.set(w, 0);
				ree.set(c, 1);
			}
		}}
	}}
}

//! produce a score between 0 and 1 for each class at each location
Idx<double> facenet::postprocess_output(double threshold, int objsize) {
	//! find candidates at each scale
	{ idx_bloop2(out, outputs, void*, resu, results, void*) {
		Idx<double> outx = ((state_idx*) out.get())->x;
		int c = outx.dim(0);
		int h = outx.dim(1);
		int w = outx.dim(2);
		Idx<double> in(c, 2 + w, 2 + h);
		Idx<double> inc(in);
		inc = inc.narrow(1, w, 1);
		inc = inc.narrow(2, h, 1);
		Idx<double> m(c, h, w);
		idx_fill(in, 0.0);
		idx_clip(outx, 0.0, inc);
		//! smooth
		{ idx_bloop2(ie, in, double, me, m, double) {
			idx_2dconvol(ie, smoothing_kernel, me);
		}}
		idx_copy(m, inc);
		//! find points that are local maxima spatial and class-wise
		//! write result in m. rescale result to [0 1]
		mark_maxima(in, inc, *((Idx<double>*) resu.get()), threshold);
	}}
	//! now prune by scale and make a list
	Idx<double> rlist(1, 6);
	rlist.resize(0, rlist.dim(1));
	Idx<double> in0x(((state_idx*) inputs.get(0))->x);
	intg s0j = in0x.dim(2);
	{ idx_bloop3(input, inputs, void*, output, outputs, void*, r, results, void*) {
		Idx<double> inx(((state_idx*) input.get())->x);