Add procrustes tests, add ability to create ParticleShapeStatistics from Project.

akenmorris · akenmorris · commit 74d8c28bb66b · 2022-04-07T10:42:42.000-06:00
diff --git a/Libs/Particles/CMakeLists.txt b/Libs/Particles/CMakeLists.txt
@@ -33,6 +33,7 @@ target_link_libraries(Particles PUBLIC
   Mesh
   tinyxml
   Eigen3::Eigen
+  Project
   )
 
 # set
diff --git a/Libs/Particles/ParticleShapeStatistics.cpp b/Libs/Particles/ParticleShapeStatistics.cpp
@@ -1,6 +1,8 @@
 
 #include "ParticleShapeStatistics.h"
 #include "ShapeEvaluation.h"
+#include <Libs/Project/Project.h>
+
 
 #include <vnl/algo/vnl_symmetric_eigensystem.h>
 #include <tinyxml.h>
@@ -342,6 +344,28 @@ int ParticleShapeStatistics::ReadPointFiles(const std::string &s)
   return 0;
 }
 
+//---------------------------------------------------------------------------
+ParticleShapeStatistics::ParticleShapeStatistics(std::shared_ptr<Project> project) {
+
+  std::vector<Eigen::VectorXd> points;
+  std::vector<int> groups;
+  for (auto& s : project->get_subjects()) {
+    auto world_files = s->get_world_particle_filenames();
+    Eigen::VectorXd particles;
+    for (auto& file : world_files) {
+      Eigen::VectorXd domain_particles;
+      ParticleSystem::ReadParticleFile(file, domain_particles);
+      Eigen::VectorXd combined(particles.size() + domain_particles.size());
+      combined << particles, domain_particles;
+      particles = combined;
+    }
+    points.push_back(particles);
+    groups.push_back(1);
+  }
+  ImportPoints(points, groups);
+}
+
+//---------------------------------------------------------------------------
 int ParticleShapeStatistics::DoPCA(std::vector<std::vector<Point>> global_pts, int domainsPerShape)
 {
   this->m_domainsPerShape = domainsPerShape;
diff --git a/Libs/Particles/ParticleShapeStatistics.h b/Libs/Particles/ParticleShapeStatistics.h
@@ -16,6 +16,8 @@
 
 namespace shapeworks {
 
+class Project;
+
 /**
   * \class ParticleShapeStatistics
   * This class computes various statistics for a set of correspondence positions
@@ -27,8 +29,9 @@ class ParticleShapeStatistics {
 
   constexpr static int VDimension = 3;
 
-  ParticleShapeStatistics() {}
-  ~ParticleShapeStatistics() {}
+  ParticleShapeStatistics() {};
+  ParticleShapeStatistics(std::shared_ptr<Project> project);
+  ~ParticleShapeStatistics() {};
 
   int DoPCA(std::vector<std::vector<Point>> global_pts, int domainsPerShape = 1);
 
diff --git a/Libs/Particles/ParticleSystem.cpp b/Libs/Particles/ParticleSystem.cpp
@@ -120,4 +120,35 @@ bool ParticleSystem::EvaluationCompare(const ParticleSystem& other) const
   return good;
 }
 
+bool ParticleSystem::ReadParticleFile(std::string filename, Eigen::VectorXd &points)
+{
+  std::ifstream in(filename);
+  if (!in.good()) {
+    return false;
+  }
+  vtkSmartPointer<vtkPoints> vtk_points = vtkSmartPointer<vtkPoints>::New();
+  int num_points = 0;
+  while (in.good()) {
+    double x, y, z;
+    in >> x >> y >> z;
+    if (!in.good()) {
+      break;
+    }
+    vtk_points->InsertNextPoint(x, y, z);
+    num_points++;
+  }
+  in.close();
+  points.setZero();
+  points.resize(num_points * 3);
+
+  int idx = 0;
+  for (int i = 0; i < num_points; i++) {
+    double* pos = vtk_points->GetPoint(i);
+    points[idx++] = pos[0];
+    points[idx++] = pos[1];
+    points[idx++] = pos[2];
+  }
+  return true;
+}
+
 }
diff --git a/Libs/Particles/ParticleSystem.h b/Libs/Particles/ParticleSystem.h
@@ -42,6 +42,8 @@ class ParticleSystem
 
   bool EvaluationCompare(const ParticleSystem& other) const;
 
+  static bool ReadParticleFile(std::string filename, Eigen::VectorXd& points);
+
 private:
   friend struct SharedCommandData;
 
diff --git a/Studio/src/Data/Shape.cpp b/Studio/src/Data/Shape.cpp
@@ -477,33 +477,7 @@ void Shape::generate_meshes(std::vector<std::string> filenames, MeshGroup& mesh_
 
 //---------------------------------------------------------------------------
 bool Shape::import_point_file(QString filename, Eigen::VectorXd& points) {
-  std::ifstream in(filename.toStdString().c_str());
-  if (!in.good()) {
-    return false;
-  }
-  vtkSmartPointer<vtkPoints> vtk_points = vtkSmartPointer<vtkPoints>::New();
-  int num_points = 0;
-  while (in.good()) {
-    double x, y, z;
-    in >> x >> y >> z;
-    if (!in.good()) {
-      break;
-    }
-    vtk_points->InsertNextPoint(x, y, z);
-    num_points++;
-  }
-  in.close();
-  points.setZero();
-  points.resize(num_points * 3);
-
-  int idx = 0;
-  for (int i = 0; i < num_points; i++) {
-    double* pos = vtk_points->GetPoint(i);
-    points[idx++] = pos[0];
-    points[idx++] = pos[1];
-    points[idx++] = pos[2];
-  }
-  return true;
+  return ParticleSystem::ReadParticleFile(filename.toStdString(), points);
 }
 
 //---------------------------------------------------------------------------
diff --git a/Testing/OptimizeTests/OptimizeTests.cpp b/Testing/OptimizeTests/OptimizeTests.cpp
@@ -632,3 +632,94 @@ TEST(OptimizeTests, mesh_ffc_test) {
   bool good = check_constraint_violations(app, 3.0e-1);
   ASSERT_TRUE(good);
 }
+
+//---------------------------------------------------------------------------
+TEST(OptimizeTests, procrustes_disabled_test) {
+  TestUtils::Instance().prep_temp(std::string(TEST_DATA_DIR) + "/optimize/procrustes", "procrustes_disabled_test");
+
+  ProjectHandle project = std::make_shared<Project>();
+  project->load("procrustes.xlsx");
+  OptimizeParameters params(project);
+  params.set_use_procrustes(false);
+
+  Optimize app;
+  params.set_up_optimize(&app);
+
+  // run optimize
+  bool success = app.Run();
+  ASSERT_TRUE(success);
+
+  // compute stats
+  ParticleShapeStatistics stats(project);
+  stats.ComputeModes();
+  stats.PrincipalComponentProjections();
+
+  // print out eigenvalues (for debugging)
+  auto values = stats.Eigenvalues();
+  for (int i = 0; i < values.size(); i++) {
+    std::cerr << "Eigenvalue " << i << " : " << values[i] << "\n";
+  }
+  ASSERT_GT(values[values.size() - 1], 700.0);
+}
+
+//---------------------------------------------------------------------------
+TEST(OptimizeTests, procrustes_no_scale_test) {
+  TestUtils::Instance().prep_temp(std::string(TEST_DATA_DIR) + "/optimize/procrustes", "procrustes_no_scale_test");
+
+  ProjectHandle project = std::make_shared<Project>();
+  project->load("procrustes.xlsx");
+  OptimizeParameters params(project);
+  params.set_use_procrustes(true);
+  params.set_use_procrustes_scaling(false);
+
+  Optimize app;
+  params.set_up_optimize(&app);
+
+  // run optimize
+  bool success = app.Run();
+  ASSERT_TRUE(success);
+
+  // compute stats
+  ParticleShapeStatistics stats(project);
+  stats.ComputeModes();
+  stats.PrincipalComponentProjections();
+
+  // print out eigenvalues (for debugging)
+  auto values = stats.Eigenvalues();
+  for (int i = 0; i < values.size(); i++) {
+    std::cerr << "Eigenvalue " << i << " : " << values[i] << "\n";
+  }
+  ASSERT_GT(values[values.size() - 1], 150.0);
+  ASSERT_LT(values[values.size() - 1], 200.0);
+}
+
+//---------------------------------------------------------------------------
+TEST(OptimizeTests, procrustes_enabled_test) {
+  TestUtils::Instance().prep_temp(std::string(TEST_DATA_DIR) + "/optimize/procrustes", "procrustes_enabled_test");
+
+  ProjectHandle project = std::make_shared<Project>();
+  project->load("procrustes.xlsx");
+  OptimizeParameters params(project);
+  params.set_use_procrustes(true);
+  params.set_use_procrustes_scaling(true);
+
+  Optimize app;
+  params.set_up_optimize(&app);
+
+  // run optimize
+  bool success = app.Run();
+  ASSERT_TRUE(success);
+
+  // compute stats
+  ParticleShapeStatistics stats(project);
+  stats.ComputeModes();
+  stats.PrincipalComponentProjections();
+
+  // print out eigenvalues (for debugging)
+  auto values = stats.Eigenvalues();
+  for (int i = 0; i < values.size(); i++) {
+    std::cerr << "Eigenvalue " << i << " : " << values[i] << "\n";
+  }
+  // should be tiny with all of procrustes enabled
+  ASSERT_LT(values[values.size() - 1], 1.0);
+}
diff --git a/Testing/Testing.cpp b/Testing/Testing.cpp
@@ -33,7 +33,7 @@ TestUtils::~TestUtils() {
   if (!temp_base_.empty() && !should_keep_dir()) {
     try {
       boost::filesystem::remove_all(temp_base_);
-    } catch (std::exception &e) {
+    } catch (std::exception& e) {
       std::cerr << "Error removing test output directory: " << temp_base_ << "\n";
       std::cerr << "Error: " << e.what();
     }
@@ -59,6 +59,14 @@ std::string TestUtils::get_output_dir(std::string test_name) {
   return test_temp.generic_string();
 }
 
+//-----------------------------------------------------------------------------
+void TestUtils::prep_temp(std::string test_folder, std::string test_name) {
+  auto temp_dir = get_output_dir(test_name);
+
+  recursive_copy(test_folder, temp_dir);
+  boost::filesystem::current_path(temp_dir);
+}
+
 //-----------------------------------------------------------------------------
 bool TestUtils::should_keep_dir() {
   bool keep = false;
@@ -72,4 +80,22 @@ bool TestUtils::should_keep_dir() {
   return keep;
 }
 
+//-----------------------------------------------------------------------------
+void TestUtils::recursive_copy(const boost::filesystem::path& src, const boost::filesystem::path& dst) {
+  // adapted from https://gist.github.com/ssteffl/30055ac128bb92801568899ebad73365
+
+  if (boost::filesystem::is_directory(src)) {
+    boost::filesystem::create_directories(dst);
+    for (boost::filesystem::directory_entry& item : boost::filesystem::directory_iterator(src)) {
+      recursive_copy(item.path(), dst / item.path().filename());
+    }
+  } else if (boost::filesystem::is_regular_file(src)) {
+    boost::filesystem::copy(src, dst);
+  } else {
+    throw std::runtime_error(dst.generic_string() + " not dir or file");
+  }
+}
+
+//-----------------------------------------------------------------------------
+
 }  // namespace shapeworks
diff --git a/Testing/Testing.h b/Testing/Testing.h
@@ -20,9 +20,14 @@ class TestUtils {
   //! Return output directory
   std::string get_output_dir(std::string test_name);
 
+  //! Create a copy of a test data folder in a temp folder and chdir there
+  void prep_temp(std::string test_folder, std::string test_name);
+
  private:
   bool should_keep_dir();
 
+  void recursive_copy(const boost::filesystem::path &src, const boost::filesystem::path &dst);
+
   boost::filesystem::path temp_base_;
 };
 }  // namespace shapeworks
diff --git a/Testing/data/optimize/procrustes/la1.ply b/Testing/data/optimize/procrustes/la1.ply
diff --git a/Testing/data/optimize/procrustes/la2.ply b/Testing/data/optimize/procrustes/la2.ply
diff --git a/Testing/data/optimize/procrustes/la3.ply b/Testing/data/optimize/procrustes/la3.ply
diff --git a/Testing/data/optimize/procrustes/procrustes.xlsx b/Testing/data/optimize/procrustes/procrustes.xlsx

Original file line number	Diff line number	Diff line change
`@@ -33,6 +33,7 @@ target_link_libraries(Particles PUBLIC`
`33`	`33`	`Mesh`
`34`	`34`	`tinyxml`
`35`	`35`	`Eigen3::Eigen`
	`36`	`+ Project`
`36`	`37`	`)`
`37`	`38`
`38`	`39`	`# set`