C++ PCA类代码示例

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();
}

void test_pca9685(PCA & pca)
{
	//Rise
	float val = 0.0f;
	for(; val < 1.0; val += 0.001f)
	{		
		for(unsigned i = 0; i < 15; ++i)
			pca.set_channel_pulse_width(i, val);
	}
	
	//Fall
	for(; val > 0.0; val -= 0.001f)
		for(unsigned i = 0; i < 15; ++i)
			pca.set_channel_pulse_width(i, val);
}

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;
}

void compressFeature( string filename, std::vector< std::vector<float> > &models, const int dim, bool ascii ){
  PCA pca;
  pca.read( filename.c_str(), ascii );
  VectorXf variance = pca.getVariance();
  MatrixXf tmpMat = pca.getAxis();
  MatrixXf tmpMat2 = tmpMat.block(0,0,tmpMat.rows(),dim);
  const int num = (int)models.size();
  for( int i=0; i<num; i++ ){
    Map<VectorXf> vec( &(models[i][0]), models[i].size() );
    //vec = tmpMat2.transpose() * vec;
    VectorXf tmpvec = tmpMat2.transpose() * vec;
    models[i].resize( dim );
    if( WHITENING ){
      for( int t=0; t<dim; t++ )
	models[i][t] = tmpvec[t] / sqrt( variance( t ) );
    }
    else{
      for( int t=0; t<dim; t++ )
	models[i][t] = tmpvec[t];
    }
  }
}

void PCANormalEstimator::estimateNormals(vtkPolyData* data, double neighRadius)
{
	vtkPointLocator* locator = vtkPointLocator::New();
	  locator->SetDataSet(data);
	  locator->BuildLocator();

	int i, numOfPoints = data->GetNumberOfPoints();
	vtkDoubleArray* normals = vtkDoubleArray::New();
	  normals->SetNumberOfComponents(3);
	  normals->SetNumberOfTuples(numOfPoints);

	PCA pca;
	double p[3], com[3], eigenvals[3], eigenvecs[3][3];
	vtkPoints* points = data->GetPoints();
	vtkIdList* neighs = vtkIdList::New();

	for ( i = 0 ; i < numOfPoints ; ++i )
	{
		points->GetPoint(i, p);
		neighs->Reset();
		locator->FindPointsWithinRadius(neighRadius, p, neighs);
		if ( neighs->GetNumberOfIds() < 3 )
		{
			normals->SetTuple3(i, 0.0, 0.0, 0.0);
			continue;
		}
		// Perform PCA
		pca.doPCA(points, neighs, eigenvecs, eigenvals, com);
		// Save the normal
		normals->SetTuple3(i, eigenvecs[0][2], eigenvecs[1][2], eigenvecs[2][2]);
	}

	data->GetPointData()->SetNormals(normals);

	// Clean up
	locator->Delete();
	normals->Delete();
	neighs->Delete();
}

//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);
}

void MTT::doPCA(RoadGraph* roads, PCA& pca) {
	// 頂点データを使って行列を生成する
	cv::Mat vmat = cv::Mat(GraphUtil::getNumVertices(roads), 2, CV_64FC1);

	int count = 0;
	RoadVertexIter vi, vend;
	for (boost::tie(vi, vend) = boost::vertices(roads->graph); vi != vend; ++vi) {
		if (!roads->graph[*vi]->valid) continue;

		vmat.at<double>(count, 0) = roads->graph[*vi]->getPt().x();
		vmat.at<double>(count, 1) = roads->graph[*vi]->getPt().y();

		count++;
	}

	// PCAを実施
	pca.pca(vmat, false);
}

	// 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;

	}

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];
}

//********************************
//* main
int main(int argc, char* argv[]) {
  if( (argc != 13) && (argc != 15) ){
    std::cerr << "usage: " << argv[0] << " [path] <rank_num> <exist_voxel_num_threshold> [model_pca_filename] <dim_model> <size1> <size2> <size3> <detect_th> <distance_th> <model_num> /input:=/camera/rgb/points" << std::endl;
    exit( EXIT_FAILURE );
  }
  char tmpname[ 1000 ];
  ros::init (argc, argv, "detect_object_vosch_multi", ros::init_options::AnonymousName);

  // read the length of voxel side
  sprintf( tmpname, "%s/param/parameters.txt", argv[1] );
  voxel_size = Param::readVoxelSize( tmpname );

  detect_th = atof( argv[9] );
  distance_th = atof( argv[10] );
  model_num = atoi( argv[11] );
  rank_num = atoi( argv[2] );

  // set marker color
  const int marker_model_num = 6;
  if( model_num > marker_model_num ){
    std::cerr << "Please set more marker colors for detection of more than " << marker_model_num << " objects." << std::endl;
    exit( EXIT_FAILURE );
  }
  marker_color_r = new float[ marker_model_num ];
  marker_color_g = new float[ marker_model_num ];
  marker_color_b = new float[ marker_model_num ];
  marker_color_r[ 0 ] = 1.0; marker_color_g[ 0 ] = 0.0; marker_color_b[ 0 ] = 0.0;  // red
  marker_color_r[ 1 ] = 0.0; marker_color_g[ 1 ] = 1.0; marker_color_b[ 1 ] = 0.0;  // green
  marker_color_r[ 2 ] = 0.0; marker_color_g[ 2 ] = 0.0; marker_color_b[ 2 ] = 1.0;  // blue
  marker_color_r[ 3 ] = 1.0; marker_color_g[ 3 ] = 1.0; marker_color_b[ 3 ] = 0.0;  // yellow
  marker_color_r[ 4 ] = 0.0; marker_color_g[ 4 ] = 1.0; marker_color_b[ 4 ] = 1.0;  // cyan
  marker_color_r[ 5 ] = 1.0; marker_color_g[ 5 ] = 0.0; marker_color_b[ 5 ] = 1.0;  // magenta
  // marker_color_r[ 0 ] = 0.0; marker_color_g[ 0 ] = 1.0; marker_color_b[ 0 ] = 0.0; // green
  // marker_color_r[ 1 ] = 0.0; marker_color_g[ 1 ] = 0.0; marker_color_b[ 1 ] = 1.0; // blue
  // marker_color_r[ 2 ] = 0.0; marker_color_g[ 2 ] = 1.0; marker_color_b[ 2 ] = 1.0; // cyan
  // marker_color_r[ 3 ] = 1.0; marker_color_g[ 3 ] = 0.0; marker_color_b[ 3 ] = 0.0; // pink

  // read the number of voxels in each subdivision's side of scene
  box_size = Param::readBoxSizeScene( tmpname );

  // read the dimension of compressed feature vectors
  dim = Param::readDim( tmpname );
  const int dim_model = atoi(argv[5]);
  if( dim <= dim_model ){
    std::cerr << "ERR: dim_model should be less than dim(in dim.txt)" << std::endl;
    exit( EXIT_FAILURE );
  }

  // read the threshold for RGB binalize
  sprintf( tmpname, "%s/param/color_threshold.txt", argv[1] );
  Param::readColorThreshold( color_threshold_r, color_threshold_g, color_threshold_b, tmpname );

  // determine the size of sliding box
  region_size = box_size * voxel_size;
  float tmp_val = atof(argv[6]) / region_size;
  int size1 = (int)tmp_val;
  if( ( ( tmp_val - size1 ) >= 0.5 ) || ( size1 == 0 ) ) size1++;
  tmp_val = atof(argv[7]) / region_size;
  int size2 = (int)tmp_val;
  if( ( ( tmp_val - size2 ) >= 0.5 ) || ( size2 == 0 ) ) size2++;
  tmp_val = atof(argv[8]) / region_size;
  int size3 = (int)tmp_val;
  if( ( ( tmp_val - size3 ) >= 0.5 ) || ( size3 == 0 ) ) size3++;
  sliding_box_size_x = size1 * region_size;
  sliding_box_size_y = size2 * region_size;
  sliding_box_size_z = size3 * region_size;

  // set variables
  search_obj.setModelNum( model_num );
#ifdef CCHLAC_TEST
  sprintf( tmpname, "%s/param/max_c.txt", argv[1] );
#else
  sprintf( tmpname, "%s/param/max_r.txt", argv[1] );
#endif
  search_obj.setNormalizeVal( tmpname );
  search_obj.setRange( size1, size2, size3 );
  search_obj.setRank( rank_num * model_num ); // for removeOverlap()
  search_obj.setThreshold( atoi(argv[3]) );

  // read projection axes of the target objects' subspace
  FILE *fp = fopen( argv[4], "r" );
  char **model_file_names = new char* [ model_num ];
  char line[ 1000 ];
  for( int i=0; i<model_num; i++ ){
    model_file_names[ i ] = new char [ 1000 ];
    if( fgets( line, sizeof(line), fp ) == NULL ) std::cerr<<"fgets err"<<std::endl;
    line[ strlen( line ) - 1 ] = '\0';
    //fscanf( fp, "%s\n", model_file_names + i );
    //model_file_names[ i ] = line;
    sprintf( model_file_names[ i ], "%s", line );
    //std::cout << model_file_names[ i ] << std::endl;
  }
  fclose(fp);
  search_obj.readAxis( model_file_names, dim, dim_model, ASCII_MODE_P, MULTIPLE_SIMILARITY );

  // read projection axis for feature compression
  PCA pca;
  sprintf( tmpname, "%s/models/compress_axis", argv[1] );
//.........这里部分代码省略.........

int main(int argc, char** argv)
{
	
	string filter = "Gaussian";
	string descriptor = "HAOG";
	string database = "CUFSF";
	uint count = 0;
	
	vector<string> extraPhotos, photos, sketches;
	
	loadImages(argv[1], photos);
	loadImages(argv[2], sketches);
	//loadImages(argv[7], extraPhotos);
	
	uint nPhotos = photos.size(),
	nSketches = sketches.size(),
	nExtra = extraPhotos.size();
	
	uint nTraining = 2*nPhotos/3;
	
	cout << "Read " << nSketches << " sketches." << endl;
	cout << "Read " << nPhotos + nExtra << " photos." << endl;
	
	vector<Mat*> sketchesDescriptors(nSketches), photosDescriptors(nPhotos), extraDescriptors(nExtra);
	
	Mat img, temp;
	
	int size=32, delta=16;
	
	#pragma omp parallel for private(img, temp)
	for(uint i=0; i<nSketches; i++){
		img = imread(sketches[i],0);
		sketchesDescriptors[i] = new Mat();
		
		#pragma omp critical
		temp = extractDescriptors(img, size, delta, filter, descriptor);
		
		*(sketchesDescriptors[i]) = temp.clone();
	}
	
	#pragma omp parallel for private(img, temp)
	for(uint i=0; i<nPhotos; i++){
		img = imread(photos[i],0);
		photosDescriptors[i] = new Mat();
		
		#pragma omp critical
		temp = extractDescriptors(img, size, delta, filter, descriptor);
		
		*(photosDescriptors[i]) = temp.clone();
	}
	
	#pragma omp parallel for private(img, temp)
	for(uint i=0; i<nExtra; i++){
		img = imread(extraPhotos[i],0);
		extraDescriptors[i] = new Mat();
		
		#pragma omp critical
		temp = extractDescriptors(img, size, delta, filter, descriptor);
		
		*(extraDescriptors[i]) = temp.clone();
	}
	
	auto seed = unsigned(count);
	
	srand(seed);
	random_shuffle(sketchesDescriptors.begin(), sketchesDescriptors.end());
	srand(seed);
	random_shuffle(photosDescriptors.begin(), photosDescriptors.end());
	
	//training
	vector<Mat*> trainingSketchesDescriptors, trainingPhotosDescriptors;
	
	trainingSketchesDescriptors.insert(trainingSketchesDescriptors.end(), sketchesDescriptors.begin(), sketchesDescriptors.begin()+nTraining);
	trainingPhotosDescriptors.insert(trainingPhotosDescriptors.end(), photosDescriptors.begin(), photosDescriptors.begin()+nTraining);
	
	//testing
	vector<Mat*> testingSketchesDescriptors, testingPhotosDescriptors;
	
	testingSketchesDescriptors.insert(testingSketchesDescriptors.end(), sketchesDescriptors.begin()+nTraining, sketchesDescriptors.end());
	testingPhotosDescriptors.insert(testingPhotosDescriptors.end(), photosDescriptors.begin()+nTraining, photosDescriptors.end());
	testingPhotosDescriptors.insert(testingPhotosDescriptors.end(), extraDescriptors.begin(), extraDescriptors.end());
	
	PCA pca;
	LDA lda;
	vector<int> labels;
	
	uint nTestingSketches = testingSketchesDescriptors.size(),
	nTestingPhotos = testingPhotosDescriptors.size();
	
	for(uint i=0; i<nTraining; i++){
		labels.push_back(i);
	}
	labels.insert(labels.end(),labels.begin(),labels.end());
	
	//bags
	vector<Mat*> testingSketchesDescriptorsBag(nTestingSketches), testingPhotosDescriptorsBag(nTestingPhotos);
	
	for(int b=0; b<200; b++){
		
		vector<int> bag_indexes = gen_bag(154, 0.1);
//.........这里部分代码省略.........

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;
    }
}