C++ PCA::project方法代码示例

int main(int argc, char** argv)
{
	
	float a[]={	0.234400724097, 0.445210153051, 0.420883079381, 0.0584111370634, 0.930917795284, 0.463946380108, 0.827477442854, 0.195052690912,
		0.224843236267, 0.011674592046, 0.778465234345, 0.795119566607, 0.834330061452, 0.250878254601, 0.907848368295, 0.159768191396,
		0.359447753375, 0.694377176768, 0.323279688498, 0.590454463022, 0.32053508251, 0.25926247011, 0.473382632749, 0.680857359827,
		0.871843303433, 0.347550207092, 0.807721675262, 0.51342440135, 0.633862634367, 0.588847708996, 0.604920986251, 0.9485023141,
		0.511286105241, 0.780677021392, 0.346168472115, 0.408572254219, 0.977881372787, 0.994457177414, 0.553713182589, 0.181657338197,
		0.188679332574, 0.138351555791, 0.549762090688, 0.763422732648, 0.270469815182, 0.368094710756, 0.28652717945, 0.344130955251,
		0.808703681865, 0.48242375244, 0.0961390490465, 0.585178232015, 0.0947071702324, 0.00663925147531, 0.409282147388, 0.865532591897,
		0.233760414088, 0.399258033215, 0.547551739688, 0.078241816204, 0.672857401346, 0.083814529556, 0.68575517509, 0.213487218459	};
		
		vector<int> labels = {1,1,1,1,1,2,1,2};
		
		Mat data(Size(8,8), CV_32F, a);
		
		//cout << data << endl;
		
		PCA pca;
		LDA lda;
		
		pca(data, Mat(), CV_PCA_DATA_AS_ROW, 7);
		
		//cout << pca.project(data) << endl;
		
		lda.compute(pca.project(data), labels);
		
		cout << lda.project(pca.project(data)) << endl;
		//cout << lda.reconstruct(lda.project(pca.project(data))) << endl;
		
		return 0;
}

//Rate a location on how likely it is to be a bubble
double rateBubble(Mat& det_img_gray, Point bubble_location, PCA& my_PCA){
	Mat query_pixels, pca_components;
	getRectSubPix(det_img_gray, Point(14,18), bubble_location, query_pixels);
	query_pixels.reshape(0,1).convertTo(query_pixels, CV_32F);
	pca_components = my_PCA.project(query_pixels);
	//The rating is the SSD of query pixels and their back projection
	Mat out = my_PCA.backProject(pca_components)- query_pixels;
	return sum(out.mul(out)).val[0];
}

void columbiaTest(int testId = 0)
{
	CascadeClassifier classifier;
	classifier.load("haarcascades/haarcascade_frontalface_alt_tree.xml");

	ShapePredictor predictor;
	predictor.load("model/helen.txt");

	PCA pca;
	FileStorage fs("model/pca.xml", FileStorage::READ);
	fs["mean"] >> pca.mean;
	fs["eigenvals"] >> pca.eigenvalues;
	fs["eigenvecs"] >> pca.eigenvectors;
	fs.release();

	/*LDA lda;
	lda.load("model/lda.xml");*/

	SVM svm;
	svm.load("model/svm.xml");

	cout << "\nmodel loaded" << endl;

	// test prediction
	cout << "\nbegin test" << endl;
	int corr = 0, total = 0;

	Mat_<float> labels, multihog, ldmks;
	collectData(testId, classifier, predictor,
				labels, multihog, ldmks);

	for (int i = 0; i < multihog.rows; i++) {
		Mat_<float> pcaVec = pca.project(multihog.row(i));
		Mat_<float> datVec(1, pcaVec.cols + ldmks.cols);
		for (int j = 0; j < pcaVec.cols; j++)
			datVec(0, j) = pcaVec(0, j);
		for (int j = 0; j < ldmks.cols; j++)
			datVec(0, j + pcaVec.cols) = ldmks(i, j);
		//Mat_<float> ldaVec = lda.project(datVec);

		float pred = svm.predict(datVec);
		if ((int)pred == (int)labels(i, 0))
			corr++;

		total++;
	}
	cout << "testId = " << testId << endl;
	cout << "corr = " << corr << " , total = " << total << endl;
	cout << "percentage: " << (double)corr / total << endl;

	ofstream fout("data/testId" + to_string(testId) + ".txt");
	fout << "corr = " << corr << " , total = " << total << endl;
	fout << "percentage: " << (double)corr / total << endl;
	fout.close();
}

//We should probably encapsulate the PCA stuff...
void train_PCA_classifier(Mat& train_img_gray, PCA& my_PCA, Mat& comparison_vectors){
	//Goes though all the selected bubble locations and puts their pixels into rows of
	//a giant matrix called so we can perform PCA on them (we need at least 3 locations to do this)
	Mat PCA_set = Mat::zeros(bubble_locations.size(), 18*14, CV_32F);
	for(size_t i = 0; i < bubble_locations.size(); i+=1) {
		Mat PCA_set_row;
		getRectSubPix(train_img_gray, Point(14,18), bubble_locations[i], PCA_set_row);
		PCA_set_row.convertTo(PCA_set_row, CV_32F);
		PCA_set.row(i) += PCA_set_row.reshape(0,1);
	}

	my_PCA = PCA(PCA_set, Mat(), CV_PCA_DATA_AS_ROW, 5);
	comparison_vectors = my_PCA.project(PCA_set);
}

	// according to http://www.pcl-users.org/Finding-oriented-bounding-box-of-a-cloud-td4024616.html
	FitCube PCLTools::fitBox(PointCloud<PointXYZ>::Ptr cloud) {

		FitCube retCube;
		PCA<PointXYZ> pca;
		PointCloud<PointXYZ> proj;

		pca.setInputCloud(cloud);
		pca.project(*cloud, proj);

		PointXYZ proj_min;
		PointXYZ proj_max;
		getMinMax3D(proj, proj_min, proj_max);

		PointXYZ min;
		PointXYZ max;
		pca.reconstruct(proj_min, min);
		pca.reconstruct(proj_max, max);

		// Rotation of PCA
		Eigen::Matrix3f rot_mat = pca.getEigenVectors();

		// Translation of PCA
		Eigen::Vector3f cl_translation = pca.getMean().head(3);

		Eigen::Matrix3f affine_trans;

		// Reordering of principal components
		affine_trans.col(2) << (rot_mat.col(0).cross(rot_mat.col(1))).normalized();
		affine_trans.col(0) << rot_mat.col(0);
		affine_trans.col(1) << rot_mat.col(1);

		retCube.rotation = Eigen::Quaternionf(affine_trans);
		Eigen::Vector4f t = pca.getMean();

		retCube.translation = Eigen::Vector3f(t.x(), t.y(), t.z());

		retCube.width = fabs(proj_max.x - proj_min.x);
		retCube.height = fabs(proj_max.y - proj_min.y);
		retCube.depth = fabs(proj_max.z - proj_min.z);

		return retCube;

	}

int
predict(PCA &pca, Mat &train, Mat face, double threshold = 0.3, double distance = 44000){
    Mat w = pca.project(face);
    Mat predicted = pca.backProject(w);

    int label = -1;
    double min = distance * threshold;

    for(int i = 0; i < 15; i++){
        for(int j = 0; j < 8; j++){
            double d = norm(predicted, train.row((i*8) + j));
            if(d < min){
                min = d;
                label = i;
            }
        }
    }

    return label;
}

vector<float> ReconnaissanceHandler::recognisePCA(vector<float>& caracteristicVector, PCA pca, vector<vector<float>>& classes)
{
	Mat v(caracteristicVector.size(),1, CV_32F);
	int i = 0;
	for (float f : caracteristicVector) {
		v.at<float>(i, 0) = f;
		++i;
	}
	Mat reduceVector = pca.project(v);
	
	vector<float> caractReduced;
	for (i = 0; i < reduceVector.rows; ++i) {
		caractReduced.push_back(reduceVector.at<float>(i, 0));
	}

	vector<float> dist;
	for (int i = 0; i < classes.size(); ++i)
	{
		dist.push_back(distanceVector(caractReduced, classes.at(i)));
	}
	return dist;
}

//Compare the bubbles with all the bubbles used in the classifier.
bubble_val checkBubble(Mat& det_img_gray, Point bubble_location, PCA& my_PCA, Mat& comparison_vectors, Point search_window=Point(5,5)){
	Mat query_pixels;

	//This bit of code finds the location in the search_window most likely to be a bubble 
	//then it checks that rather than the exact specified location.
	Mat out = Mat::zeros(Size(search_window.y, search_window.x) , CV_32FC1);
	Point offset = Point(bubble_location.x - search_window.x/2, bubble_location.y - search_window.y/2);
	for(size_t i = 0; i < search_window.y; i+=1) {
		for(size_t j = 0; j < search_window.x; j+=1) {
			out.row(i).col(j) += rateBubble(det_img_gray, Point(j,i) + offset, my_PCA);
		}
	}
	Point min_location;
	minMaxLoc(out, NULL,NULL, &min_location);

	getRectSubPix(det_img_gray, Point(14,18), min_location + offset, query_pixels);

	query_pixels.reshape(0,1).convertTo(query_pixels, CV_32F);
	Point max_location;
	matchTemplate(comparison_vectors, my_PCA.project(query_pixels), out, CV_TM_CCOEFF_NORMED);
	minMaxLoc(out, NULL,NULL,NULL, &max_location);
	return bubble_values[max_location.y];
}

void
FMR_FNMR(PCA &pca, Mat &train, Mat &test, double step = 0.05, double distance = 44000.0){

    struct stat {
        double total{0}, match{0};
        double res() { return match/total; }
    };

    Mat w, r;

    cout << "Threshold;" << "FMR;" << "FNMR;" << "ROC TP;" << "ROC FP"  << endl;

    for(double t = 0.0; t <= 1.000001; t += step) {
        double min = distance * t;

        stat fmr, fnmr;

        for(int i = 0; i < 15; i++){
            for(int j = 0; j < 3; j++){

                Mat face = test.row((i*3) + j);

                pca.project(face, w);
                pca.backProject(w, r);

                for(int k = 0; k < 15; k++) {
                    if(k == i){
                        // FNMR
                        for(int l = 0; l < 8; l++) {
                            double d = norm(r, train.row((k*8) + l));
                            if(d > min){
                                fnmr.match++;
                            }
//                            else {
                                // ROC TP
//                            }
                            fnmr.total++;
                        }

                    } else {
                        // FMR
                        for(int l = 0; l < 8; l++){
                            double d = norm(r, train.row((k*8) + l));
                            if(d < min){
                                fmr.match++;

                                // ROC FP
                            }
                            fmr.total++;
                        }
                    }
                }
            }
        }

        cout << 1-t << ";"
             << fmr.res() << ";"
             << fnmr.res() << ";"
             << ((fnmr.total-fnmr.match)/fnmr.total) << ";"
             << fmr.res() << endl;
    }
}

int main(int argc, char* argv[])
{
	string trainlistfile;
	string testlistfile;
	string modelfile;
	string rawfeaturefile;
	string resultpath;
	string kmeansfilepath;
	int maxComponent = 32;
	int cluster_num = 512;
	int feature_dimention = 128;
	if (argc == 1)
	{
		help();
		return -1;
	}
	for (int i = 1; i < argc;++i)
	{
		
		if (string(argv[i])=="--trainlist")
		{
			trainlistfile = argv[++i];
		}
		else if (string(argv[i]) == "--testlist")
		{
			testlistfile = argv[++i];
		}
		else if (string(argv[i]) == "--clusternum")
		{
			cluster_num = std::stoi(argv[++i]);
		}
		else if (string(argv[i]) == "--featuredim")
		{
			feature_dimention = std::stoi(argv[++i]);
		}
		else if (string(argv[i]) == "--modelfile")
		{
			modelfile = argv[++i];
		}
		else if (string(argv[i]) == "--rawfeature")
		{
			rawfeaturefile = argv[++i];
		}
		else if (string(argv[i])=="--resultpath")
		{
			resultpath = argv[++i];

		}
		else if (string(argv[i]) == "--maxComponent")
		{
			maxComponent = std::stoi(argv[++i]);

	}
		else if (string(argv[i]) == "--kmeansfilepath")
		{
			kmeansfilepath = argv[++i];

		}
	}
#ifdef kmeans_pca
	vlad::configure(256,128);
	PCA pca;
	vlad::loadPCAmodel("pca_model.yml",pca);
	cout<<pca.eigenvalues<<endl;
#endif
#ifdef kmeans
	vlad::getKmeansModel(cluster_num,feature_dimention,rawfeaturefile,resultpath);
#endif
#ifdef VLAD
	Stopwatch watch;
	watch.Start();
	vlad::GetVladFeatureFromSift(kmeansfilepath,testlistfile);
	watch.Stop();
	cout<<"getkmeans use "<<watch.GetTime()<<endl;
	getchar();
#endif
#ifdef ReduceMatrix
	ifstream rawfeature(rawfeaturefile.c_str());
	ofstream outPut(resultpath.c_str());
	string single_raw_fea;
	PCA pca;
	vlad::loadPCAmodel(modelfile, pca);
	while (getline(rawfeature,single_raw_fea))
	{
		vector<string> ww;
		cv::Mat rawfeature_mat = cv::Mat(1,LENGTH_OF_SINGLE_DATA, CV_32F);
		split_words(single_raw_fea, " ", ww);
		for (int j = 0; j < LENGTH_OF_SINGLE_DATA; ++j){
			rawfeature_mat.at<float>(0, j) = atof(ww[j].c_str());
		}
		cv::Mat matric_reduced;
		matric_reduced = pca.project(rawfeature_mat);
		for (int i=0;i<matric_reduced.cols;++i)
		{
			outPut<<matric_reduced.at<float>(0,i);
		}
		outPut <<endl;
        
	}
	rawfeature.close();
//.........这里部分代码省略.........

void test(set<int> &testSet, int code)
{
	CascadeClassifier classifier;
	classifier.load("haarcascades/haarcascade_frontalface_alt_tree.xml");

	ShapePredictor predictor;
	predictor.load("model/helen.txt");

	PCA pca;
	FileStorage fs("model/girl_pca.xml", FileStorage::READ);
	fs["mean"] >> pca.mean;
	fs["eigenvals"] >> pca.eigenvalues;
	fs["eigenvecs"] >> pca.eigenvectors;

	SVM svm;
	svm.load("model/girl_svm.xml");

	cout << "\nmodel loaded" << endl;

	ifstream fin("img/labels.txt");
	ofstream fout("data/out_" + 
				  to_string(code) + ".txt");
	VideoWriter writer("data/out.avi", 0, 10, Size(1920, 1080), true);

	string line;
	int corr = 0, total = 0;
	while (getline(fin, line)) {
		stringstream ss(line);
		int frame, label;
		ss >> frame >> label;
		label -= 49;

		if (testSet.find(frame) == testSet.end())
			continue;

		Mat vis = imread("img/" + to_string(frame) + ".jpg",
						 CV_LOAD_IMAGE_UNCHANGED);
		Mat_<uchar> img;
		cvtColor(vis, img, COLOR_BGR2GRAY);
		BBox bbox = getTestBBox(img, classifier);
		if (EmptyBox(bbox)) continue;

		Mat_<double> shape = predictor(img, bbox);
		Geom G;	initGeom(shape, G);
		Pose P; calcPose(G, P);

		Mat_<uchar> lEye, rEye;
		regularize(img, bbox, P, shape, lEye, rEye);

		vector<float> lRlt;
		vector<float> rRlt;
		calcMultiHog(lEye, lRlt);
		calcMultiHog(rEye, rRlt);

		Mat_<float> pcaVec, ldmks;

		vector<float> _hog2nd_vec;
		for (int k = 0; k < lRlt.size(); k++)
			_hog2nd_vec.push_back(lRlt[k]);
		for (int k = 0; k < rRlt.size(); k++)
			_hog2nd_vec.push_back(rRlt[k]);
		Mat_<float> multihog = Mat_<float>(_hog2nd_vec).reshape(1, 1);
		pcaVec = pca.project(multihog);

		vector<float> _ldmks;
		for (int i = 28; i < 48; i++) {
			_ldmks.push_back((shape(i, 0) - bbox.cx) / bbox.w);
			_ldmks.push_back((shape(i, 1) - bbox.cy) / bbox.h);
		}
		float mouthx = (shape(51, 0) + shape(62, 0) + shape(66, 0) + shape(57, 0)) / 4;
		float mouthy = (shape(51, 1) + shape(62, 1) + shape(66, 1) + shape(57, 1)) / 4;
		_ldmks.push_back((mouthx - bbox.cx) / bbox.w);
		_ldmks.push_back((mouthy - bbox.cy) / bbox.h);
		float maxVal = *std::max_element(_ldmks.begin(), _ldmks.end());
		for (int i = 0; i < _ldmks.size(); i++) _ldmks[i] *= 1.0 / maxVal; // scale to [-1, 1]
		ldmks = Mat_<float>(_ldmks).reshape(1, 1);

		Mat_<float> sample(1, pcaVec.cols + ldmks.cols);
		for (int j = 0; j < pcaVec.cols; j++)
			sample(0, j) = pcaVec(0, j);
		for (int j = 0; j < ldmks.cols; j++)
			sample(0, j + pcaVec.cols) = ldmks(0, j);

		int pred = svm.predict(sample);
		if (pred == label) corr++;
		total++;

		fout << frame << ' ' << label << ' ' << pred << endl;

		string s1, s2;
		switch (label) {
		case 0: s1 = "annotation: Eye"; break;
		case 1: s1 = "annotation: Face"; break;
		case 2: s1 = "annotation: NOF"; break;
		}
		switch (pred) {
		case 0: s2 = "prediction: Eye"; break;
		case 1: s2 = "prediction: Face"; break;
		case 2: s2 = "prediction: NOF"; break;
		}
//.........这里部分代码省略.........