Merge pull request #2127 from SCIInstitute/2124-vtk_particles

2124 vtk particles

Author: akenmorris

Applications/shapeworks/Commands.cpp

@@ -102,7 +102,10 @@ void OptimizeCommand::buildParser() {
 
   parser.add_option("--name").action("store").type("string").set_default("").help("Path to project file.");
   parser.add_option("--progress").action("store_true").set_default(false).help("Show progress [default: false].");
-  parser.add_option("--xmlconsole").action("store_true").set_default(false).help("XML console output [default: false].");
+  parser.add_option("--xmlconsole")
+      .action("store_true")
+      .set_default(false)
+      .help("XML console output [default: false].");
 
   Command::buildParser();
 }
@@ -133,7 +136,6 @@ bool OptimizeCommand::execute(const optparse::Values& options, SharedCommandData
         return false;
       }
 
-
       const auto oldBasePath = boost::filesystem::current_path();
       auto base = StringUtils::getPath(projectFile);
       if (base != projectFile) {
@@ -179,7 +181,10 @@ void GroomCommand::buildParser() {
 
   parser.add_option("--name").action("store").type("string").set_default("").help("Path to project file.");
   parser.add_option("--progress").action("store_true").set_default(false).help("Show progress [default: false].");
-  parser.add_option("--xmlconsole").action("store_true").set_default(false).help("XML console output [default: false].");
+  parser.add_option("--xmlconsole")
+      .action("store_true")
+      .set_default(false)
+      .help("XML console output [default: false].");
 
   Command::buildParser();
 }
@@ -236,8 +241,10 @@ void AnalyzeCommand::buildParser() {
 
   parser.add_option("--name").action("store").type("string").set_default("").help("Path to project file.");
   parser.add_option("--output").action("store").type("string").set_default("").help("Path to output file.");
-  parser.add_option("--range").action("store").type("float").set_default(3.0f).help("Standard deviation range for PCA [default: 3.0].");
-  parser.add_option("--steps").action("store").type("int").set_default(21).help("Number of steps to use for PCA [default: 21].");
+  parser.add_option("--range").action("store").type("float").set_default(3.0f).help(
+      "Standard deviation range for PCA [default: 3.0].");
+  parser.add_option("--steps").action("store").type("int").set_default(21).help(
+      "Number of steps to use for PCA [default: 21].");
 
   Command::buildParser();
 }

Libs/Optimize/Optimize.cpp

@@ -174,10 +174,6 @@ class Optimize {
   void SetOptimizationIterationsCompleted(int optimization_iterations_completed);
   //! Set the number of iterations per split
   void SetIterationsPerSplit(int iterations_per_split);
-  //! Set the init criterion (TODO: details)
-  void SetInitializationCriterion(double init_criterion);
-  //! Set the optimization criterion (TODO: details)
-  void SetOptimizationCriterion(double opt_criterion);
   //! Set if shape statistics should be used in initialization
   void SetUseShapeStatisticsInInit(bool use_shape_statistics_in_init);
   //! Set the interval for running procrustes (0 to disable)
@@ -208,9 +204,9 @@ class Optimize {
   //! Set if mixed effects should be used (TODO: details)
   void SetUseMixedEffects(bool use_mixed_effects);
 
-  //! Set if optimization should be done using disentangled spatiotemporal SSM method 
+  //! Set if optimization should be done using disentangled spatiotemporal SSM method
   void SetUseDisentangledSpatiotemporalSSM(bool use_disentangled_ssm_4d);
-  //! Returns true if optimization is to be done using disentangled spatiotemporal SSM method 
+  //! Returns true if optimization is to be done using disentangled spatiotemporal SSM method
   bool GetUseDisentangledSpatiotemporalSSM();
 
   //! For good/bad analysis, set the normal angle to use (TODO: details)
@@ -291,7 +287,7 @@ class Optimize {
   //! n * number_of_triangles
   void SetGeodesicsCacheSizeMultiplier(size_t n);
 
-  shapeworks::OptimizationVisualizer& GetVisualizer();
+  OptimizationVisualizer& GetVisualizer();
   void SetShowVisualizer(bool show);
   bool GetShowVisualizer();
 
@@ -302,6 +298,8 @@ class Optimize {
 
   void UpdateProgress();
 
+  void set_particle_format(std::string format) { particle_format_ = format; }
+
  protected:
   //! Set the iteration callback. Derived classes should override to set their own callback
   virtual void SetIterationCallback();
@@ -399,8 +397,6 @@ class Optimize {
   int m_optimization_iterations = 2000;
   int m_optimization_iterations_completed = 0;
   int m_iterations_per_split = 1000;
-  double m_initialization_criterion = 1e-6;
-  double m_optimization_criterion = 1e-6;
   bool m_use_shape_statistics_in_init = false;
   unsigned int m_procrustes_interval = 3;
   bool m_procrustes_scaling = true;
@@ -459,6 +455,7 @@ class Optimize {
   bool show_visualizer_ = false;
   shapeworks::OptimizationVisualizer visualizer_;
 
+  std::string particle_format_ = "particles";
   std::shared_ptr<Project> project_;
 
   std::chrono::system_clock::time_point m_start_time;

Libs/Optimize/Optimize.h

@@ -339,16 +339,6 @@ bool OptimizeParameterFile::set_optimization_parameters(TiXmlHandle* docHandle,
     optimize->SetIterationsPerSplit(atoi(elem->GetText()));
   }
 
-  elem = docHandle->FirstChild("init_criterion").Element();
-  if (elem) {
-    optimize->SetInitializationCriterion(atof(elem->GetText()));
-  }
-
-  elem = docHandle->FirstChild("opt_criterion").Element();
-  if (elem) {
-    optimize->SetOptimizationCriterion(atof(elem->GetText()));
-  }
-
   elem = docHandle->FirstChild("use_shape_statistics_in_init").Element();
   if (elem) {
     optimize->SetUseShapeStatisticsInInit((bool)atoi(elem->GetText()));

Libs/Optimize/OptimizeParameterFile.cpp

@@ -48,6 +48,7 @@ const std::string field_attributes = "field_attributes";
 const std::string field_attribute_weights = "field_attribute_weights";
 const std::string use_geodesics_to_landmarks = "use_geodesics_to_landmarks";
 const std::string geodesics_to_landmarks_weight = "geodesics_to_landmarks_weight";
+const std::string particle_format = "particle_format";
 }  // namespace Keys
 
 //---------------------------------------------------------------------------
@@ -87,9 +88,8 @@ OptimizeParameters::OptimizeParameters(ProjectHandle project) {
                                          Keys::use_geodesics_to_landmarks,
                                          Keys::geodesics_to_landmarks_weight,
                                          Keys::keep_checkpoints,
-                                         Keys::use_disentangled_ssm
-
-  };
+                                         Keys::use_disentangled_ssm,
+                                         Keys::particle_format};
 
   // check if params_ has any unknown keys
   for (auto& param : params_.get_map()) {
@@ -348,6 +348,7 @@ bool OptimizeParameters::set_up_optimize(Optimize* optimize) {
   optimize->SetOutputDir(get_output_prefix());
   optimize->SetMeshFFCMode(get_mesh_ffc_mode());
   optimize->SetUseDisentangledSpatiotemporalSSM(get_use_disentangled_ssm());
+  optimize->set_particle_format(get_particle_format());
 
   // TODO Remove this once Studio has controls for shared boundary
   optimize->SetSharedBoundaryEnabled(true);
@@ -616,9 +617,10 @@ bool OptimizeParameters::set_up_optimize(Optimize* optimize) {
 
       auto name = StringUtils::getBaseFilenameWithoutExtension(filename);
 
+      auto extension = get_particle_format();
       auto prefix = get_output_prefix();
-      local_particle_filenames.push_back(prefix + name + "_local.particles");
-      world_particle_filenames.push_back(prefix + name + "_world.particles");
+      local_particle_filenames.push_back(prefix + name + "_local." + extension);
+      world_particle_filenames.push_back(prefix + name + "_world." + extension);
     }
     s->set_local_particle_filenames(local_particle_filenames);
     s->set_world_particle_filenames(world_particle_filenames);
@@ -710,3 +712,9 @@ double OptimizeParameters::get_geodesic_to_landmarks_weight() {
 void OptimizeParameters::set_geodesic_to_landmarks_weight(double value) {
   params_.set(Keys::geodesics_to_landmarks_weight, value);
 }
+
+//---------------------------------------------------------------------------
+std::string OptimizeParameters::get_particle_format() { return params_.get(Keys::particle_format, "particles"); }
+
+//---------------------------------------------------------------------------
+void OptimizeParameters::set_particle_format(std::string format) { params_.set(Keys::particle_format, format); }

Libs/Optimize/OptimizeParameters.cpp

@@ -124,6 +124,10 @@ class OptimizeParameters {
   double get_geodesic_to_landmarks_weight();
   void set_geodesic_to_landmarks_weight(double value);
 
+  std::string get_particle_format();
+  void set_particle_format(std::string format);
+
+
  private:
   std::string get_output_prefix();
 

Libs/Optimize/OptimizeParameters.h

@@ -2,12 +2,11 @@
 #include "Sampler.h"
 
 #include <Logging.h>
+#include <Particles/ParticleFile.h>
 
 #include "Libs/Optimize/Domain/ContourDomain.h"
-#include "Libs/Optimize/Domain/ImageDomain.h"
 #include "Libs/Optimize/Utils/ObjectReader.h"
-#include "itkImageRegionIterator.h"
-#include "itkParticlePositionReader.h"
+
 
 namespace shapeworks {
 
@@ -73,7 +72,6 @@ void Sampler::AllocateDomainsAndNeighborhoods() {
   // *after* registering the attributes to the particle system since some of
   // them respond to AddDomain.
   // Here, the Constraints actually get added to the constraints class
-  int ctr = 0;
   for (unsigned int i = 0; i < this->m_DomainList.size(); i++) {
     auto domain = m_DomainList[i];
 
@@ -91,8 +89,6 @@ void Sampler::AllocateDomainsAndNeighborhoods() {
         }
       }
 
-      auto imageDomain = static_cast<ImplicitSurfaceDomain<ImageType::PixelType>*>(domain.get());
-
       // Adding free-form constraints to constraint object
       // std::cout << "m_FFCs.size() " << m_FFCs.size() << std::endl;
       if (m_FFCs.size() > i) {
@@ -133,10 +129,8 @@ void Sampler::ReadPointsFiles() {
   // If points file names have been specified, then read the initial points.
   for (unsigned int i = 0; i < m_PointsFiles.size(); i++) {
     if (m_PointsFiles[i] != "") {
-      ParticlePositionReader::Pointer reader = ParticlePositionReader::New();
-      reader->SetFileName(m_PointsFiles[i].c_str());
-      reader->Update();
-      this->GetParticleSystem()->AddPositionList(reader->GetOutput(), i);
+      auto points = particles::read_particles_as_vector(m_PointsFiles[i]);
+      this->GetParticleSystem()->AddPositionList(points, i);
     }
   }
 
@@ -148,7 +142,6 @@ void Sampler::ReadPointsFiles() {
 void Sampler::InitializeOptimizationFunctions() {
   // Set the minimum neighborhood radius and maximum sigma based on the
   // domain of the 1st input image.
-  unsigned int maxdim = 0;
   double maxradius = -1.0;
   double minimumNeighborhoodRadius = this->m_Spacing;
 

Libs/Optimize/Sampler.cpp

@@ -1,9 +1,10 @@
 #pragma once
 
+#include <Particles/ParticleFile.h>
+
 #include <string>
 
-#include "Libs/Optimize/Domain/ImageDomainWithGradients.h"
-#include "itkParticlePositionWriter.h"
+//#include "Libs/Optimize/Domain/ImageDomainWithGradients.h"
 #include "vnl/algo/vnl_symmetric_eigensystem.h"
 
 namespace shapeworks {
@@ -74,11 +75,7 @@ void ParticleGaussianModeWriter<VDimension>::Update() const {
       meanlist.push_back(p);
     }
 
-    std::string tmpstr = fn + ".mean";
-    ParticlePositionWriter::Pointer writer = ParticlePositionWriter::New();
-    writer->SetFileName(tmpstr.c_str());
-    writer->SetInput(meanlist);
-    writer->Update();
+    particles::write_particles_from_vector(fn + ".mean", meanlist);
 
     int modenum = 0;
     for (int mode = num_samples - 1; mode > num_samples - (m_NumberOfModes + 1); mode--, modenum++) {
@@ -99,10 +96,9 @@ void ParticleGaussianModeWriter<VDimension>::Update() const {
 
         char fp[255];
         ::sprintf(fp, ".%d.%d.mode", s, modenum);
-        tmpstr = fn + fp;
-        writer->SetFileName(tmpstr.c_str());
-        writer->SetInput(modelist);
-        writer->Update();
+        std::string tmpstr = fn + fp;
+
+        particles::write_particles_from_vector(tmpstr, modelist);
       }
     }
   }

Libs/Optimize/Utils/ParticleGaussianModeWriter.txx

@@ -1,23 +1,23 @@
 # Particles library
 
 set(Particles_sources
-        ParticleSystemEvaluation.cpp
- ParticleShapeStatistics.cpp
- ShapeEvaluation.cpp
- ReconstructSurface.cpp
- ParticleNormalEvaluation.cpp
- )
+  ParticleSystemEvaluation.cpp
+  ParticleShapeStatistics.cpp
+  ShapeEvaluation.cpp
+  ReconstructSurface.cpp
+  ParticleNormalEvaluation.cpp
+  ParticleFile.cpp
+  )
 
 set(Particles_headers
-        ParticleSystemEvaluation.h
- ParticleShapeStatistics.h
- itkParticlePositionReader.h
- itkParticlePositionWriter.h
- EvaluationUtil.h
- ShapeEvaluation.h
- ReconstructSurface.h
- ParticleNormalEvaluation.h
- )
+  ParticleSystemEvaluation.h
+  ParticleShapeStatistics.h
+  EvaluationUtil.h
+  ShapeEvaluation.h
+  ReconstructSurface.h
+  ParticleNormalEvaluation.h
+  ParticleFile.h
+  )
 
 add_library(Particles STATIC
   ${Particles_sources}