Merge pull request #2170 from SCIInstitute/exclude_mode

Add exclude capability

Author: akenmorris

Libs/Analyze/Particles.cpp

@@ -248,6 +248,9 @@ void Particles::save_particles_file(std::string filename, const Eigen::VectorXd&
   out.close();
 }
 
+//---------------------------------------------------------------------------
+int Particles::get_number_of_domains() { return local_particles_.size(); }
+
 //---------------------------------------------------------------------------
 void Particles::set_alignment_type(int alignment) { alignment_type_ = alignment; }
 

Libs/Analyze/Particles.h

@@ -54,6 +54,7 @@ class Particles {
 
   static void save_particles_file(std::string filename, const Eigen::VectorXd& points);
 
+  int get_number_of_domains();
 
  private:
   void transform_global_particles();

Libs/Analyze/Shape.cpp

@@ -71,15 +71,7 @@ void Shape::set_subject(std::shared_ptr<Subject> subject) {
   groomed_meshes_.set_number_of_meshes(subject->get_number_of_domains());
   reconstructed_meshes_.set_number_of_meshes(subject->get_number_of_domains());
 
-  if (!subject_->get_original_filenames().empty()) {
-    /// TODO: Show multiple lines of filenames for multiple domains?
-    std::string filename = subject_->get_original_filenames()[0];
-    corner_annotations_[0] = StringUtils::getFilename(filename);
-  }
-
-  if (subject->get_display_name() != "") {
-    corner_annotations_[0] = subject->get_display_name();
-  }
+  update_annotations();
 }
 
 //---------------------------------------------------------------------------
@@ -88,6 +80,12 @@ bool Shape::is_subject() { return subject_ != nullptr; }
 //---------------------------------------------------------------------------
 std::shared_ptr<Subject> Shape::get_subject() { return subject_; }
 
+//---------------------------------------------------------------------------
+bool Shape::is_fixed() { return is_subject() && subject_->is_fixed(); }
+
+//---------------------------------------------------------------------------
+bool Shape::is_excluded() { return is_subject() && subject_->is_excluded(); }
+
 //---------------------------------------------------------------------------
 void Shape::import_original_file(const std::string& filename) {
   subject_->set_original_filenames({filename});
@@ -320,10 +318,26 @@ int Shape::get_id() { return id_; }
 void Shape::set_id(int id) { id_ = id; }
 
 //---------------------------------------------------------------------------
-void Shape::update_name() {
+void Shape::update_annotations() {
+  if (!subject_->get_original_filenames().empty()) {
+    /// TODO: Show multiple lines of filenames for multiple domains?
+    std::string filename = subject_->get_original_filenames()[0];
+    corner_annotations_[0] = StringUtils::getFilename(filename);
+  }
+
   if (subject_->get_display_name() != "") {
     corner_annotations_[0] = subject_->get_display_name();
   }
+
+  std::vector<std::string> corner2_items;
+  if (subject_->is_fixed()) {
+    corner2_items.push_back("Fixed");
+  }
+  if (subject_->is_excluded()) {
+    corner2_items.push_back("Excluded");
+  }
+
+  corner_annotations_[2] = StringUtils::join(corner2_items, ", ");
 }
 
 //---------------------------------------------------------------------------

Libs/Analyze/Shape.h

@@ -56,6 +56,12 @@ class Shape {
   //! Return the pointer to the subject object
   std::shared_ptr<shapeworks::Subject> get_subject();
 
+  //! Helper to ask if this shape is fixed or not
+  bool is_fixed();
+
+  /// Helper to ask if this shape is excluded
+  bool is_excluded();
+
   //! Import the original raw mesh or image file
   void import_original_file(const std::string& filename);
 
@@ -119,7 +125,7 @@ class Shape {
   void set_id(int id);
 
   //! Update the name of this shape
-  void update_name();
+  void update_annotations();
 
   std::vector<std::string> get_original_filenames();
   std::vector<std::string> get_original_filenames_with_path();
@@ -201,7 +207,7 @@ class Shape {
   std::map<std::string, Eigen::VectorXd> point_features_;
   Particles particles_;
 
-  std::shared_ptr<shapeworks::Subject> subject_;
+  std::shared_ptr<Subject> subject_;
 
   vtkSmartPointer<vtkTransform> transform_ = vtkSmartPointer<vtkTransform>::New();
 

Libs/Common/Logging.cpp

@@ -10,11 +10,17 @@ namespace spd = spdlog;
 namespace shapeworks {
 
 //-----------------------------------------------------------------------------
-static std::string create_header(const int line, const char* filename) {
+static std::string create_header(const int line, const char* filename, const char* function = "") {
   const char* name = (strrchr(filename, '/') ? strrchr(filename, '/') + 1 : filename);
   const char* name2 = (strrchr(name, '\\') ? strrchr(name, '\\') + 1 : name);
-  std::string header = "[" + std::string(name2) + "|" + std::to_string(line) + "]";
-  return header;
+  const char* function_name = (strrchr(function, ':') ? strrchr(function, ':') + 1 : function);
+  if (!function) {
+    std::string header = "[" + std::string(name2) + "|" + std::to_string(line) + "]";
+    return header;
+  } else {
+    std::string header = "[" + std::string(name2) + "|" + std::string(function_name) + "|" + std::to_string(line) + "]";
+    return header;
+  }
 }
 
 //-----------------------------------------------------------------------------
@@ -98,8 +104,8 @@ void Logging::show_status(const std::string& message, const int line, const char
 }
 
 //-----------------------------------------------------------------------------
-void Logging::log_debug(const std::string& message, const int line, const char* file) const {
-  std::string str = create_header(line, file) + " " + message;
+void Logging::log_debug(const std::string& message, const int line, const char* file, const char* function) const {
+  std::string str = create_header(line, file, function) + " " + message;
   spd::debug(str);
   if (log_open_) {
     spd::get("file")->debug(str);
@@ -155,6 +161,8 @@ void Logging::set_debug_callback(const std::function<void(std::string)>& callbac
 void Logging::set_status_callback(const std::function<void(std::string)>& callback) { status_callback_ = callback; }
 
 //-----------------------------------------------------------------------------
-void Logging::set_progress_callback(const std::function<void(double, std::string)>& callback) { progress_callback_ = callback; }
+void Logging::set_progress_callback(const std::function<void(double, std::string)>& callback) {
+  progress_callback_ = callback;
+}
 
 }  // namespace shapeworks

Libs/Common/Logging.h

@@ -105,7 +105,7 @@ class Logging {
   void show_progress(double value, const std::string& message);
 
   //! Log a debug message, use SW_DEBUG macro
-  void log_debug(const std::string& message, const int line, const char* file) const;
+  void log_debug(const std::string& message, const int line, const char* file, const char *function) const;
 
   //! Log a warning message, use SW_WARN macro
   void log_warning(const std::string& message, const int line, const char* file) const;
@@ -168,7 +168,7 @@ class Logging {
 
 //! Log debug macro
 #define SW_DEBUG(message, ...) \
-  shapeworks::Logging::Instance().log_debug(fmt::format(message, ##__VA_ARGS__), __LINE__, __FILE__)
+  shapeworks::Logging::Instance().log_debug(fmt::format(message, ##__VA_ARGS__), __LINE__, __FILE__, __FUNCTION__)
 
 //! Variable trace macro (e.g. output variable name = <variable value>)
 #define SW_TRACE(x) SW_DEBUG(#x " = {}", x);

Libs/Groom/Groom.cpp

@@ -52,6 +52,12 @@ bool Groom::run() {
           continue;
         }
 
+        if (subjects[i]->is_excluded()) {
+          // clear groomed filenames
+          subjects[i]->set_groomed_filenames(std::vector<std::string>());
+          return;
+        }
+
         bool is_image = project_->get_original_domain_types()[domain] == DomainType::Image;
         bool is_mesh = project_->get_original_domain_types()[domain] == DomainType::Mesh;
         bool is_contour = project_->get_original_domain_types()[domain] == DomainType::Contour;
@@ -511,9 +517,11 @@ bool Groom::run_alignment() {
         vtkSmartPointer<vtkTransform> transform = ProjectUtils::convert_transform(list);
         mesh.applyTransform(transform);
 
+        // if fixed subjects are present, only add the fixed subjects
         if (subjects[i]->is_fixed() || !project_->get_fixed_subjects_present()) {
-          // if fixed subjects are present, only add the fixed subjects
-          reference_meshes.push_back(mesh);
+          if (!subjects[i]->is_excluded()) {
+            reference_meshes.push_back(mesh);
+          }
         }
         meshes.push_back(mesh);
       }
@@ -550,8 +558,11 @@ bool Groom::run_alignment() {
       auto list = subjects[i]->get_groomed_transforms()[0];
       vtkSmartPointer<vtkTransform> transform = ProjectUtils::convert_transform(list);
       mesh.applyTransform(transform);
+      // if fixed subjects are present, only add the fixed subjects
       if (subjects[i]->is_fixed() || !project_->get_fixed_subjects_present()) {
-        meshes.push_back(mesh);
+        if (!subjects[i]->is_excluded()) {
+          meshes.push_back(mesh);
+        }
       }
     }
 
@@ -590,6 +601,9 @@ void Groom::assign_transforms(std::vector<std::vector<double>> transforms, int d
   auto subjects = project_->get_subjects();
 
   for (size_t i = 0; i < subjects.size(); i++) {
+    if (subjects[i]->is_excluded()) {
+      continue;
+    }
     auto subject = subjects[i];
 
     int base_domain = domain;

Libs/Optimize/Domain/ImageDomain.h

@@ -137,7 +137,7 @@ class ImageDomain : public ParticleRegionDomain {
     } else {
       std::ostringstream message;
       message << "Domain " << m_DomainID << ": " << m_DomainName << " : Distance transform queried for a Point, " << p
-              << ", outside the given image domain. Consider increasing the narrow band";
+              << ", outside the given image domain. Consider increasing the padding in grooming or the narrow band optimization parameter";
       throw std::runtime_error(message.str());
     }
   }

Libs/Optimize/Matrix/LegacyShapeMatrix.h

@@ -32,14 +32,14 @@ class LegacyShapeMatrix : public vnl_matrix<double>, public Observer {
   /** Method for creation through the object factory. */
   itkNewMacro(Self)
 
-  /** Run-time type information (and related methods). */
-  itkTypeMacro(LegacyShapeMatrix, Observer)
-
-  /** Callbacks that may be defined by a subclass.  If a subclass defines one
-       of these callback methods, the corresponding flag in m_DefinedCallbacks
-       should be set to true so that the ParticleSystem will know to register
-       the appropriate event with this method. */
-  virtual void DomainAddEventCallback(Object*, const itk::EventObject& e) {
+      /** Run-time type information (and related methods). */
+      itkTypeMacro(LegacyShapeMatrix, Observer)
+
+      /** Callbacks that may be defined by a subclass.  If a subclass defines one
+           of these callback methods, the corresponding flag in m_DefinedCallbacks
+           should be set to true so that the ParticleSystem will know to register
+           the appropriate event with this method. */
+      virtual void DomainAddEventCallback(Object*, const itk::EventObject& e) {
     const ParticleDomainAddEvent& event = dynamic_cast<const ParticleDomainAddEvent&>(e);
     unsigned int d = event.GetDomainIndex();
 
@@ -67,16 +67,28 @@ class LegacyShapeMatrix : public vnl_matrix<double>, public Observer {
     const typename ParticleSystem::PointType pos = ps->GetTransformedPosition(idx, d);
 
     int numRows = 0;
-    for (int i = 0; i < m_DomainsPerShape; i++) numRows += VDimension * ps->GetNumberOfParticles(i);
+    for (int i = 0; i < m_DomainsPerShape; i++) {
+      numRows += VDimension * ps->GetNumberOfParticles(i);
+    }
 
-    if (numRows > this->rows()) this->ResizeMatrix(numRows, this->cols());
+    if (numRows > this->rows()) {
+      this->ResizeMatrix(numRows, this->cols());
+    }
 
     unsigned int k = 0;
     int dom = d % m_DomainsPerShape;
-    for (int i = 0; i < dom; i++) k += VDimension * ps->GetNumberOfParticles(i);
+    for (int i = 0; i < dom; i++) {
+      k += VDimension * ps->GetNumberOfParticles(i);
+    }
     k += idx * VDimension;
 
-    for (unsigned int i = 0; i < VDimension; i++) this->operator()(i + k, d / m_DomainsPerShape) = pos[i];
+    for (unsigned int i = 0; i < VDimension; i++) {
+      if (i + k >= this->rows()) {
+        throw std::runtime_error(
+            "PositionSetEventCallback: index out of bounds! Different number of particles per shape?");
+      }
+      this->operator()(i + k, d / m_DomainsPerShape) = pos[i];
+    }
   }
 
   virtual void PositionSetEventCallback(Object* o, const itk::EventObject& e) {
@@ -89,10 +101,18 @@ class LegacyShapeMatrix : public vnl_matrix<double>, public Observer {
 
     unsigned int k = 0;
     int dom = d % m_DomainsPerShape;
-    for (int i = 0; i < dom; i++) k += VDimension * ps->GetNumberOfParticles(i);
+    for (int i = 0; i < dom; i++) {
+      k += VDimension * ps->GetNumberOfParticles(i);
+    }
     k += idx * VDimension;
 
-    for (unsigned int i = 0; i < VDimension; i++) this->operator()(i + k, d / m_DomainsPerShape) = pos[i];
+    for (unsigned int i = 0; i < VDimension; i++) {
+      if (i + k >= this->rows()) {
+        throw std::runtime_error(
+            "PositionSetEventCallback: index out of bounds! Different number of particles per shape?");
+      }
+      this->operator()(i + k, d / m_DomainsPerShape) = pos[i];
+    }
   }
 
   virtual void PositionRemoveEventCallback(Object*, const itk::EventObject&) {

Libs/Optimize/OptimizeParameters.cpp

@@ -297,8 +297,8 @@ std::vector<std::vector<itk::Point<double>>> OptimizeParameters::get_initial_poi
     for (auto s : subjects) {
       if (s->is_fixed()) {
         count++;
-        // read the local points
-        auto filename = s->get_local_particle_filenames()[d];
+        // read the world points that are in the shared coordinate space
+        auto filename = s->get_world_particle_filenames()[d];
         auto particles = read_particles_as_vector(filename);
         if (domain_sum.size() == 0) {
           domain_sum = particles;
@@ -531,13 +531,19 @@ bool OptimizeParameters::set_up_optimize(Optimize* optimize) {
     }
 
     optimize->SetFixedDomains(fixed_domains);
-    optimize->SetInitialPoints(get_initial_points());
+    if (!get_use_landmarks()) {  // can't use both initial points and landmarks
+      SW_DEBUG("Setting Initial Points");
+      optimize->SetInitialPoints(get_initial_points());
+    }
   }
 
   for (auto s : subjects) {
     if (abort_load_) {
       return false;
     }
+    if (s->is_excluded()) {
+      continue;
+    }
     auto files = s->get_groomed_filenames();
     if (files.empty()) {
       throw std::invalid_argument("No groomed inputs for optimization");
@@ -552,6 +558,7 @@ bool OptimizeParameters::set_up_optimize(Optimize* optimize) {
       point_files.insert(std::end(point_files), std::begin(landmarks), std::end(landmarks));
     }
     if (!point_files.empty()) {
+      SW_DEBUG("Setting Initial Points as landmarks");
       optimize->SetPointFiles(point_files);
     }
   }
@@ -619,6 +626,11 @@ bool OptimizeParameters::set_up_optimize(Optimize* optimize) {
     if (abort_load_) {
       return false;
     }
+
+    if (s->is_excluded()) {
+      continue;
+    }
+
     auto files = s->get_groomed_filenames();
     if (files.empty()) {
       throw std::invalid_argument("No groomed inputs for optimization");
@@ -681,9 +693,9 @@ bool OptimizeParameters::set_up_optimize(Optimize* optimize) {
       } else {
         Image image(filename);
         if (s->is_fixed()) {
-          optimize->AddImage(nullptr);
+          optimize->AddImage(nullptr, filename);
         } else {
-          optimize->AddImage(image);
+          optimize->AddImage(image, filename);
         }
       }