gistfile1.h

// MNIST: bag of keypoints
#include "nv_core.h"
#include "nv_ml.h"
#include "nv_io.h"
#include "nv_num.h"
#include "cxcore.h"
#include "cv.h"
#include "highgui.h"
#include "nv_face_feature.h"

static void conv_vec2img(nv_matrix_t *img, nv_matrix_t *vec, int vecm)
{
	int row, col, n;

	nv_matrix_zero(img);
	for (n = 0; n < vec->n; ++n) {
		row = (n / 28) + 18;
		col = (n % 28) + 18;
		NV_MAT3D_V(img, row, col, 0) = NV_MAT_V(vec, vecm, n);
	}
}

static int surf_desc(nv_matrix_t *data, int data_m, nv_matrix_t *desc_vec)
{
	nv_matrix_t *img = nv_matrix3d_alloc(1, 64, 64);
	int i, n, desc_m = 0;
	CvMemStorage* storage = cvCreateMemStorage(0);
	CvSeq *keypoints = 0, *descriptors = 0;
	IplImage *cv_img;
	IplImage *cv_img_big = cvCreateImage(cvSize(128, 128), IPL_DEPTH_8U, 1);
	CvSURFParams params = cvSURFParams(300, 1);

	nv_matrix_zero(img);

	conv_vec2img(img, data, data_m);
	cv_img = nv_conv_nv2ipl(img);
	cvResize(cv_img, cv_img_big, 1);
	cvThreshold(cv_img_big, cv_img_big, 128.0, 255.0, CV_THRESH_BINARY);

	cvExtractSURF(cv_img_big, NULL, &keypoints, &descriptors, storage, params, 0);
/*
	for (i = 0; i < descriptors->total; ++i) {
		const CvSURFPoint* p = (CvSURFPoint*)cvGetSeqElem(keypoints, i);
		cvCircle(cv_img_big, cvPoint(p->pt.x, p->pt.y), cvRound(p->size*1.2/9.*2), cvScalar(0x66, 0x66, 0x66, 0x66), 1, 8, 0);
	}
	cvNamedWindow("DEBUG",1);
	cvShowImage("DEBUG", cv_img_big);
	cvWaitKey(0);
*/

	for (i = 0; i < descriptors->total; ++i) {
		const float* desc = (const float*)cvGetSeqElem(descriptors, i);

		for (n = 0; n < 128; ++n) {
			NV_MAT_V(desc_vec, desc_m, n) = desc[n];
		}
		++desc_m;
	}

	cvReleaseImage(&cv_img);
	cvReleaseImage(&cv_img_big);
	cvReleaseMemStorage(&storage);
	nv_matrix_free(&img);

	return desc_m;
}


static void extract_surf()
{
	int i, j, m,x, y;
	//nv_matrix_t *train_data = nv_load_matrix("test_28x28.vec");
	nv_matrix_t *train_data = nv_load_matrix("train_28x28.vec");
	nv_matrix_t *data = nv_matrix_alloc(128, train_data->m * 40);
	nv_matrix_t *desc_vec = nv_matrix_alloc(128, 100);
	int data_m = 0;

	nv_matrix_zero(data);

	for (m = 0; m < train_data->m; ++m) {
		int desc_m;
		
		nv_matrix_zero(desc_vec);
		desc_m = surf_desc(train_data, m, desc_vec);
		printf("M: %d, N:%d\r", m, desc_m);
		for (i = 0; i < desc_m; ++i) {
			nv_vector_copy(data, data_m, desc_vec, i);
			++data_m;
		}
	}
	nv_matrix_m(data, data_m);
	printf("\n M:%d\n", data_m);

	nv_save_matrix_text("train_surf.vec", data);
}

#define KCLASS 2048
static void vecq_clustering()
{
	int c, m;
	nv_matrix_t *data = nv_load_matrix_text("train_surf.vec");
	nv_matrix_t *centroid = nv_matrix_alloc(data->n, KCLASS);
	nv_matrix_t *labels = nv_matrix_alloc(1, data->m);
	nv_matrix_t *maxmin = nv_matrix_alloc(2, KCLASS);
	nv_matrix_t *count = nv_matrix_alloc(1, KCLASS);

	nv_matrix_zero(centroid);
	nv_matrix_zero(count);
	nv_matrix_zero(labels);

	nv_kmeans(centroid, count, labels, data, KCLASS, 50);

	// 確認用
	for (c = 0; c < KCLASS; ++c) {
		NV_MAT_V(maxmin, c, 0) = 0.0;     // MAX
		NV_MAT_V(maxmin, c, 1) = FLT_MAX; // MIN
	}
	for (m = 0; m < data->m; ++m) {
		int label = (int)NV_MAT_V(labels, m, 0);
		float dist = nv_euclidean(centroid, label, data, m);
		if (NV_MAT_V(maxmin, label, 0) < dist) {
			NV_MAT_V(maxmin, label, 0) = dist;
		}
		if (NV_MAT_V(maxmin, label, 1) > dist) {
			NV_MAT_V(maxmin, label, 1) = dist;
		}
	}

	nv_save_matrix("surf_centroid.vec", centroid);
	nv_save_matrix_text("surf_count.vec", count);
	nv_save_matrix_text("surf_maxmin.vec", maxmin);
}

static void vw_vector(nv_matrix_t *vec, int vec_m, nv_matrix_t *centroid, nv_matrix_t *data, int m)
{
	nv_matrix_t *desc_vec = nv_matrix_alloc(128, 100);
	int i;
	int desc_m;
	
	nv_matrix_zero(desc_vec);

	desc_m = surf_desc(data, m, desc_vec);

	for (i = 0; i < desc_m; ++i) {
		int label = nv_vector_nn(centroid, desc_vec, i); //最近傍探索
		NV_MAT_V(vec, vec_m, label) += 1.0;
	}

	nv_matrix_free(&desc_vec);
}

static void make_vw_data()
{
	int j, m, n;
	//nv_matrix_t *train_data = nv_load_matrix("test_28x28.vec");
	nv_matrix_t *train_data = nv_load_matrix("train_28x28.vec");
	nv_matrix_t *data = nv_matrix_alloc(KCLASS, train_data->m);
	nv_matrix_t *centroid = nv_load_matrix("surf_centroid.vec");
	int data_m = 0;

	for (m = 0; m < train_data->m; ++m) {
		for (n = 0; n < train_data->n; ++n) {
			NV_MAT_V(data, m, n) = 0.1f;
		}
	}
#ifdef _OPENMP
#pragma omp parallel for
#endif
	for (m = 0; m < train_data->m; ++m) {
	   vw_vector(data, m, centroid, train_data, m);
	}
	nv_save_matrix_text("train_vecq.vec", data);
	//nv_save_matrix("test_vecq.vec", data);
}


static void train()
{
	nv_matrix_t *train_data = nv_load_matrix_text("train_vecq.vec");
	nv_matrix_t *train_label = nv_load_matrix("train_label.vec");
	//nv_mlp_t *mlp = nv_mlp_alloc(train_data->n, 64, 10);
	nv_mlp_t *mlp = nv_load_mlp("mnist_vecq.mlp");
	int i;

	printf("n: %d, m: %d\n", train_data->n, train_data->m);
	//nv_mlp_init(mlp);

	for (i = 500; i < 2000; i+= 5) {
		nv_mlp_train_ex(mlp, train_data, train_label, 0.001f, 0.00001f, i, i + 5, 0);
		nv_save_mlp("mnist_vecq.mlp", mlp);
	}
}

static void test()
{
	nv_matrix_t *test_data = nv_load_matrix("test_vecq.vec");
	nv_matrix_t *test_label = nv_load_matrix("test_label.vec");
	long t;
	nv_mlp_t *mlp = nv_load_mlp("mnist_vecq.mlp");
	char fname[100];
	int m, ok, ol, i;
	int result[10] = {0};

	ok = 0;
	ol = 0;
	t = nv_clock();
	for (m = 0; m < test_data->m; ++m) {
		int l = nv_mlp_predict_label(mlp, test_data, m);
		if (l == (int)NV_MAT_V(test_label, m, 0)) {
			++ok;
		}
		result[l]++;
		if ((float)ol != NV_MAT_V(test_label, m, 0)) {
			ol = (int)NV_MAT_V(test_label, m, 0);
			for (i = 0; i < 10; ++i) {
				if (i != 0) {
					printf(",");
				}
				printf("%04d", result[i]);
			}
			printf("\n");

			memset(result, 0, sizeof(result));
		}
	}
	for (i = 0; i < 10; ++i) {
		if (i != 0) {
			printf(",");
		}
		printf("%04d", result[i]);
	}
	printf("\n\nTEST ERROR RATE(%%): %f, %dms\n", 100.0 * (1.0-(float)ok / (test_data->m)), nv_clock() - t);
	printf("\n");
}


int mnist_surf(void)
{
#ifdef _OPENMP
	omp_set_num_threads(4);
#endif

	//extract_surf(); // SURF Descriptor抽出
	//vecq_clustering(); // ベクトル量子化のためのクラスタリング
	//make_vw_data(); // MNISTのデータを変換
	train(); // 分類器の学習
	//test(); // テスト
	return 0;
}