Merge pull request #2434 from SCIInstitute/python_pca

Python PCA API

Author: akenmorris

Libs/Particles/ParticleShapeStatistics.cpp

@@ -476,7 +476,10 @@ int ParticleShapeStatistics::read_point_files(const std::string& s) {
 }
 
 //---------------------------------------------------------------------------
-ParticleShapeStatistics::ParticleShapeStatistics(std::shared_ptr<Project> project) {
+ParticleShapeStatistics::ParticleShapeStatistics(std::shared_ptr<Project> project) { load_from_project(project); }
+
+//---------------------------------------------------------------------------
+void ParticleShapeStatistics::load_from_project(std::shared_ptr<Project> project) {
   std::vector<Eigen::VectorXd> points;
   std::vector<int> groups;
   for (auto& s : project->get_subjects()) {
@@ -497,7 +500,7 @@ ParticleShapeStatistics::ParticleShapeStatistics(std::shared_ptr<Project> projec
 
 //---------------------------------------------------------------------------
 int ParticleShapeStatistics::do_pca(std::vector<std::vector<Point>> global_pts, int domainsPerShape) {
-  this->domains_per_shape_ = domainsPerShape;
+  domains_per_shape_ = domainsPerShape;
 
   // Assumes all the same size.
   num_samples_ = global_pts.size() / domains_per_shape_;
@@ -508,11 +511,6 @@ int ParticleShapeStatistics::do_pca(std::vector<std::vector<Point>> global_pts,
   mean_.resize(num_dimensions_);
   mean_.fill(0);
 
-  std::cout << "VDimension = " << values_per_particle_ << "-------------\n";
-  std::cout << "m_numSamples = " << num_samples_ << "-------------\n";
-  std::cout << "m_domainsPerShape = " << domains_per_shape_ << "-------------\n";
-  std::cout << "global_pts.size() = " << global_pts.size() << "-------------\n";
-
   // Compile the "meta shapes"
   for (unsigned int i = 0; i < num_samples_; i++) {
     for (unsigned int k = 0; k < domains_per_shape_; k++) {
@@ -571,6 +569,13 @@ int ParticleShapeStatistics::do_pca(ParticleSystemEvaluation ParticleSystemEvalu
   return do_pca(particlePoints, domainsPerShape);
 }
 
+//---------------------------------------------------------------------------
+int ParticleShapeStatistics::do_pca(std::shared_ptr<Project> project) {
+  load_from_project(project);
+  compute_modes();
+  return 0;
+}
+
 //---------------------------------------------------------------------------
 int ParticleShapeStatistics::compute_modes() {
   SW_DEBUG("computing PCA modes");
@@ -638,6 +643,13 @@ int ParticleShapeStatistics::principal_component_projections() {
   return 0;
 }
 
+//---------------------------------------------------------------------------
+Eigen::VectorXd ParticleShapeStatistics::project_new_sample(const Eigen::VectorXd& new_sample) {
+  // Subtract mean and project onto principal components
+  Eigen::VectorXd centered_sample = new_sample - mean_;
+  return eigenvectors_.transpose() * centered_sample;
+}
+
 //---------------------------------------------------------------------------
 int ParticleShapeStatistics::write_csv_file(const std::string& s) {
   // Write csv file

Libs/Particles/ParticleShapeStatistics.h

@@ -19,13 +19,15 @@ class Project;
 class ParticleShapeStatistics {
  public:
   ParticleShapeStatistics(){};
-  ParticleShapeStatistics(std::shared_ptr<Project> project);
+  explicit ParticleShapeStatistics(std::shared_ptr<Project> project);
   ~ParticleShapeStatistics(){};
 
   int do_pca(std::vector<std::vector<Point>> global_pts, int domainsPerShape = 1);
 
   int do_pca(ParticleSystemEvaluation particleSystem, int domainsPerShape = 1);
 
+  int do_pca(std::shared_ptr<Project> project);
+
   //! Loads a set of point files and pre-computes some statistics.
   int import_points(std::vector<Eigen::VectorXd> points, std::vector<int> group_ids);
 
@@ -57,6 +59,9 @@ class ParticleShapeStatistics {
   //!  principal componenent axes for each of the samples.  ComputeModes must be called first.
   int principal_component_projections();
 
+  //! Projects a new sample into the PCA space defined by the original samples.
+  Eigen::VectorXd project_new_sample(const Eigen::VectorXd& new_sample);
+
   //! Returns the sample size
   int get_num_samples() const { return num_samples_; }
   int get_group1_num_samples() const { return num_samples_group1_; }
@@ -135,6 +140,8 @@ class ParticleShapeStatistics {
   //! Set the meshes for each sample (used for some evaluation metrics)
   void set_meshes(const std::vector<Mesh>& meshes) { meshes_ = meshes; }
 
+  void load_from_project(std::shared_ptr<Project> project);
+
  private:
   unsigned int num_samples_group1_;
   unsigned int num_samples_group2_;

Libs/Python/ShapeworksPython.cpp

@@ -559,8 +559,7 @@ PYBIND11_MODULE(shapeworks_py, m) {
           [](Image& image, std::vector<double>& pt) -> decltype(auto) {
             return image.evaluate(Point({pt[0], pt[1], pt[2]}));
           },
-          "evaluate the image at any given point in space", "pt"_a)
-      ;
+          "evaluate the image at any given point in space", "pt"_a);
 
   // PhysicalRegion
   py::class_<PhysicalRegion>(m, "PhysicalRegion")
@@ -1281,9 +1280,14 @@ PYBIND11_MODULE(shapeworks_py, m) {
 
       .def(py::init<>())
 
+      // int do_pca(ParticleSystemEvaluation particleSystem, int domainsPerShape = 1);
       .def("PCA", py::overload_cast<ParticleSystemEvaluation, int>(&ParticleShapeStatistics::do_pca),
            "calculates the eigen values and eigen vectors of the data", "particleSystem"_a, "domainsPerShape"_a = 1)
 
+      // int do_pca(std::shared_ptr<Project> project);
+      .def("PCA", py::overload_cast<std::shared_ptr<Project>>(&ParticleShapeStatistics::do_pca),
+           "calculates the eigen values and eigen vectors of the data from a project", "project"_a)
+
       .def("principalComponentProjections", &ParticleShapeStatistics::principal_component_projections,
            "projects the original data on the calculated principal components")
 
@@ -1302,7 +1306,12 @@ PYBIND11_MODULE(shapeworks_py, m) {
        components are constructed")
 
       .def("percentVarByMode", &ParticleShapeStatistics::get_percent_variance_by_mode,
-           "return the variance accounted for by the principal components");
+           "return the variance accounted for by the principal components")
+
+      .def("projectNewSample", &ParticleShapeStatistics::project_new_sample, "project a new sample into the PCA space",
+           "newSample"_a)
+
+      .def("getMean", &ParticleShapeStatistics::get_mean, "returns the mean shape particles");
 
   define_python_analyze(m);
   define_python_groom(m);