Merge pull request #2169 from SCIInstitute/coarse_geodesics

Coarse geodesics

Author: akenmorris

Libs/Analyze/Analyze.cpp

@@ -69,7 +69,7 @@ static void write_good_bad_angles(json& json_object, ProjectHandle project, Anal
   std::vector<json> good_bad_angles;
 
   for (int d = 0; d < project->get_number_of_domains_per_subject(); d++) {
-    std::vector<std::shared_ptr<VtkMeshWrapper>> meshes;
+    std::vector<std::shared_ptr<MeshWrapper>> meshes;
     for (auto& shape : analyze.get_shapes()) {
       meshes.push_back(shape->get_groomed_mesh_wrappers()[d]);
     }

Libs/Analyze/Shape.cpp

@@ -16,7 +16,7 @@
 #include <vtkKdTreePointLocator.h>
 #include <vtkPointData.h>
 
-#include "Libs/Optimize/Domain/VtkMeshWrapper.h"
+#include "Libs/Optimize/Domain/MeshWrapper.h"
 
 using ReaderType = itk::ImageFileReader<ImageType>;
 
@@ -819,14 +819,14 @@ bool Shape::has_planes() {
 }
 
 //---------------------------------------------------------------------------
-std::vector<std::shared_ptr<VtkMeshWrapper>> Shape::get_groomed_mesh_wrappers() {
+std::vector<std::shared_ptr<MeshWrapper>> Shape::get_groomed_mesh_wrappers() {
   if (!groomed_mesh_wrappers_.empty()) {
     return groomed_mesh_wrappers_;
   }
 
   auto group = get_groomed_meshes(true /* wait */);
   for (auto& mesh : group.meshes()) {
-    auto wrapper = std::make_shared<VtkMeshWrapper>(mesh->get_poly_data());
+    auto wrapper = std::make_shared<MeshWrapper>(mesh->get_poly_data());
     groomed_mesh_wrappers_.push_back(wrapper);
   }
   return groomed_mesh_wrappers_;

Libs/Analyze/Shape.h

@@ -22,7 +22,7 @@ namespace shapeworks {
 class Shape;
 using ShapeHandle = std::shared_ptr<Shape>;
 using ShapeList = std::vector<ShapeHandle>;
-class VtkMeshWrapper;
+class MeshWrapper;
 
 //! Representation of a single shape/patient/subject.
 class Shape {
@@ -175,7 +175,7 @@ class Shape {
 
   bool has_planes();
 
-  std::vector<std::shared_ptr<VtkMeshWrapper>> get_groomed_mesh_wrappers();
+  std::vector<std::shared_ptr<MeshWrapper>> get_groomed_mesh_wrappers();
 
  private:
   void generate_meshes(std::vector<std::string> filenames, MeshGroup& mesh_list, bool save_transform,
@@ -191,7 +191,7 @@ class Shape {
   MeshGroup original_meshes_;
   MeshGroup groomed_meshes_;
   MeshGroup reconstructed_meshes_;
-  std::vector<std::shared_ptr<VtkMeshWrapper>> groomed_mesh_wrappers_;
+  std::vector<std::shared_ptr<MeshWrapper>> groomed_mesh_wrappers_;
 
   std::string override_feature_;
 

Libs/Groom/GroomParameters.cpp

@@ -141,12 +141,18 @@ GroomParameters::GroomParameters(ProjectHandle project, std::string domain_name)
                                          Keys::CENTER,
                                          Keys::ICP};
 
+  std::vector<std::string> to_remove;
+
   // check if params_ has any unknown keys
   for (auto& param : params_.get_map()) {
     if (std::find(all_params.begin(), all_params.end(), param.first) == all_params.end()) {
       SW_WARN("Unknown Grooming parameter: " + param.first);
+      to_remove.push_back(param.first);
     }
   }
+  for (auto& param : to_remove) {
+    params_.remove_entry(param);
+  }
 }
 
 //---------------------------------------------------------------------------

Libs/Mesh/Mesh.cpp

@@ -60,7 +60,7 @@
 
 #include "FEFixMesh.h"
 #include "Image.h"
-#include "Libs/Optimize/Domain/VtkMeshWrapper.h"
+#include "Libs/Optimize/Domain/MeshWrapper.h"
 #include "Logging.h"
 #include "MeshComputeThickness.h"
 #include "MeshUtils.h"
@@ -768,7 +768,7 @@ double Mesh::geodesicDistance(int source, int target) const {
     throw std::invalid_argument("requested point ids outside range of points available in mesh");
   }
 
-  VtkMeshWrapper wrap(this->poly_data_, true);
+  MeshWrapper wrap(this->poly_data_, true);
   return wrap.ComputeDistance(getPoint(source), -1, getPoint(target), -1);
 }
 
@@ -778,7 +778,7 @@ Field Mesh::geodesicDistance(const Point3 landmark) const {
   distance->SetNumberOfTuples(numPoints());
   distance->SetName("GeodesicDistanceToLandmark");
 
-  VtkMeshWrapper wrap(this->poly_data_, true);
+  MeshWrapper wrap(this->poly_data_, true);
 
   for (int i = 0; i < numPoints(); i++) {
     distance->SetValue(i, wrap.ComputeDistance(landmark, -1, getPoint(i), -1));

Libs/Mesh/Mesh.h

@@ -270,7 +270,7 @@ class Mesh {
 
   //! Formats mesh into an IGL format
   MeshPoints getIGLMesh(Eigen::MatrixXd& V, Eigen::MatrixXi& F)
-      const;  // Copied directly from VtkMeshWrapper. this->poly_data_ becomes this->mesh. // WARNING: Copied directly
+      const;  // Copied directly from MeshWrapper. this->poly_data_ becomes this->mesh. // WARNING: Copied directly
               // from Meshwrapper. TODO: When refactoring, take this into account.
 
   //! Clips the mesh according to a field value

Libs/Optimize/Domain/MeshDomain.cpp

@@ -1,16 +1,8 @@
-/*=========================================================================
-  Copyright (c) 2009 Scientific Computing and Imaging Institute.
-  See ShapeWorksLicense.txt for details.
-
-     This software is distributed WITHOUT ANY WARRANTY; without even
-     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
-     PURPOSE.  See the above copyright notices for more information.
-=========================================================================*/
-
 #include "MeshDomain.h"
 
 namespace shapeworks {
 
+//-------------------------------------------------------------------
 bool MeshDomain::ApplyConstraints(PointType &p, int idx, bool dbg) const {
   if (!mesh_wrapper_) {
     return true;
@@ -19,17 +11,20 @@ bool MeshDomain::ApplyConstraints(PointType &p, int idx, bool dbg) const {
   return true;
 }
 
+//-------------------------------------------------------------------
 bool MeshDomain::ApplyVectorConstraints(vnl_vector_fixed<double, DIMENSION> &gradE, const PointType &pos) const {
   // TODO need to refactor into update particle method
   // TODO reduce magnitude of vector so it doesn't violate constraints
   return true;
 }
 
+//-------------------------------------------------------------------
 vnl_vector_fixed<double, DIMENSION> MeshDomain::ProjectVectorToSurfaceTangent(
     vnl_vector_fixed<double, DIMENSION> &gradE, const PointType &pos, int idx) const {
   return mesh_wrapper_->ProjectVectorToSurfaceTangent(pos, idx, gradE);
 }
 
+//-------------------------------------------------------------------
 MeshDomain::PointType MeshDomain::UpdateParticlePosition(const PointType &point, int idx,
                                                          vnl_vector_fixed<double, DIMENSION> &update) const {
   vnl_vector_fixed<double, DIMENSION> negativeUpdate;
@@ -40,6 +35,7 @@ MeshDomain::PointType MeshDomain::UpdateParticlePosition(const PointType &point,
   return newPoint;
 }
 
+//-------------------------------------------------------------------
 double MeshDomain::GetMaxDiameter() const {
   // todo should this not be the length of the bounding box diagonal?
   PointType boundingBoxSize = mesh_wrapper_->GetMeshUpperBound() - mesh_wrapper_->GetMeshLowerBound();
@@ -50,6 +46,79 @@ double MeshDomain::GetMaxDiameter() const {
   return max;
 }
 
+//-------------------------------------------------------------------
+vnl_vector_fixed<float, 3> MeshDomain::SampleGradientAtPoint(const PointType &point, int idx) const {
+  return mesh_wrapper_->SampleNormalAtPoint(point, idx);
+}
+
+//-------------------------------------------------------------------
+vnl_vector_fixed<float, 3> MeshDomain::SampleNormalAtPoint(const PointType &point, int idx) const {
+  return mesh_wrapper_->SampleNormalAtPoint(point, idx);
+}
+
+//-------------------------------------------------------------------
+ParticleDomain::GradNType MeshDomain::SampleGradNAtPoint(const PointType &p, int idx) const {
+  return mesh_wrapper_->SampleGradNAtPoint(p, idx);
+}
+
+//-------------------------------------------------------------------
+double MeshDomain::Distance(const PointType &a, int idx_a, const PointType &b, int idx_b,
+                            vnl_vector_fixed<double, 3> *out_grad) const {
+  return geodesics_mesh_->ComputeDistance(a, idx_a, b, idx_b, out_grad);
+}
+
+//-------------------------------------------------------------------
+double MeshDomain::SquaredDistance(const PointType &a, int idx_a, const PointType &b, int idx_b) const {
+  double dist = geodesics_mesh_->ComputeDistance(a, idx_a, b, idx_b);
+  return dist * dist;
+}
+
+//-------------------------------------------------------------------
+bool MeshDomain::IsWithinDistance(const PointType &a, int idx_a, const PointType &b, int idx_b, double test_dist,
+                                  double &dist) const {
+  return geodesics_mesh_->IsWithinDistance(a, idx_a, b, idx_b, test_dist, dist);
+}
+
+//-------------------------------------------------------------------
+void MeshDomain::SetMesh(std::shared_ptr<MeshWrapper> mesh, double geodesic_remesh_percent) {
+  m_FixedDomain = false;
+  mesh_wrapper_ = mesh;
+  sw_mesh_ = std::make_shared<Mesh>(mesh_wrapper_->GetPolydata());
+
+  if (geodesic_remesh_percent >= 100.0) { // no remeshing
+    geodesics_mesh_ = mesh_wrapper_;
+  } else {
+    auto poly_data = mesh_wrapper_->GetPolydata();
+    Mesh mesh_copy(poly_data);
+    mesh_copy.remeshPercent(geodesic_remesh_percent, 1.0);
+    geodesics_mesh_ = std::make_shared<MeshWrapper>(mesh_copy.getVTKMesh());
+  }
+}
+
+//-------------------------------------------------------------------
 void MeshDomain::InvalidateParticlePosition(int idx) const { this->mesh_wrapper_->InvalidateParticle(idx); }
 
+//-------------------------------------------------------------------
+ParticleDomain::PointType MeshDomain::GetZeroCrossingPoint() const {
+  // TODO Hong
+  // Apply constraints somehow
+  if (mesh_wrapper_ == nullptr) {
+    // Fixed domain. Unsure if this is the correct thing to do, but it preserves existing behaviour.
+    PointType p;
+    p[0] = p[1] = p[2] = 0;
+    return p;
+  }
+  return mesh_wrapper_->GetPointOnMesh();
+}
+
+//-------------------------------------------------------------------
+ParticleDomain::PointType MeshDomain::GetValidLocationNear(PointType p) const {
+  PointType valid;
+  valid[0] = p[0];
+  valid[1] = p[1];
+  valid[2] = p[2];
+  ApplyConstraints(valid, -1);
+  return valid;
+}
+
 }  // namespace shapeworks

Libs/Optimize/Domain/MeshDomain.h

@@ -40,26 +40,9 @@ class MeshDomain : public ParticleDomain {
   const PointType &GetLowerBound() const override { return mesh_wrapper_->GetMeshLowerBound(); }
   const PointType &GetUpperBound() const override { return mesh_wrapper_->GetMeshUpperBound(); }
 
-  PointType GetZeroCrossingPoint() const override {
-    // TODO Hong
-    // Apply constraints somehow
-    if (mesh_wrapper_ == nullptr) {
-      // Fixed domain. Unsure if this is the correct thing to do, but it preserves existing behaviour.
-      PointType p;
-      p[0] = p[1] = p[2] = 0;
-      return p;
-    }
-    return mesh_wrapper_->GetPointOnMesh();
-  }
+  PointType GetZeroCrossingPoint() const override;
 
-  PointType GetValidLocationNear(PointType p) const override {
-    PointType valid;
-    valid[0] = p[0];
-    valid[1] = p[1];
-    valid[2] = p[2];
-    ApplyConstraints(valid, -1);
-    return valid;
-  }
+  PointType GetValidLocationNear(PointType p) const override;
 
   double GetSurfaceArea() const override {
     // TODO return actual surface area
@@ -68,32 +51,19 @@ class MeshDomain : public ParticleDomain {
 
   double GetMaxDiameter() const override;
 
-  inline vnl_vector_fixed<float, DIMENSION> SampleGradientAtPoint(const PointType &point, int idx) const override {
-    return mesh_wrapper_->SampleNormalAtPoint(point, idx);
-  }
+  vnl_vector_fixed<float, DIMENSION> SampleGradientAtPoint(const PointType &point, int idx) const override;
 
-  inline vnl_vector_fixed<float, DIMENSION> SampleNormalAtPoint(const PointType &point, int idx) const override {
-    return mesh_wrapper_->SampleNormalAtPoint(point, idx);
-  }
+  vnl_vector_fixed<float, DIMENSION> SampleNormalAtPoint(const PointType &point, int idx) const override;
 
-  inline GradNType SampleGradNAtPoint(const PointType &p, int idx) const override {
-    return mesh_wrapper_->SampleGradNAtPoint(p, idx);
-  }
+  GradNType SampleGradNAtPoint(const PointType &p, int idx) const override;
 
-  inline double Distance(const PointType &a, int idx_a, const PointType &b, int idx_b,
-                         vnl_vector_fixed<double, DIMENSION> *out_grad = nullptr) const override {
-    return mesh_wrapper_->ComputeDistance(a, idx_a, b, idx_b, out_grad);
-  }
+  double Distance(const PointType &a, int idx_a, const PointType &b, int idx_b,
+                  vnl_vector_fixed<double, DIMENSION> *out_grad = nullptr) const override;
 
-  inline double SquaredDistance(const PointType &a, int idx_a, const PointType &b, int idx_b) const override {
-    double dist = mesh_wrapper_->ComputeDistance(a, idx_a, b, idx_b);
-    return dist * dist;
-  }
+  double SquaredDistance(const PointType &a, int idx_a, const PointType &b, int idx_b) const override;
 
-  inline bool IsWithinDistance(const PointType &a, int idx_a, const PointType &b, int idx_b, double test_dist,
-                               double &dist) const override {
-    return mesh_wrapper_->IsWithinDistance(a, idx_a, b, idx_b, test_dist, dist);
-  }
+  bool IsWithinDistance(const PointType &a, int idx_a, const PointType &b, int idx_b, double test_dist,
+                        double &dist) const override;
 
   void DeleteImages() override {
     // TODO Change this to a generic delete function
@@ -103,18 +73,15 @@ class MeshDomain : public ParticleDomain {
     // TODO Change this to a generic delete function
   }
 
-  void SetMesh(std::shared_ptr<shapeworks::MeshWrapper> mesh_) {
-    m_FixedDomain = false;
-    mesh_wrapper_ = mesh_;
-    sw_mesh_ = std::make_shared<Mesh>(mesh_wrapper_->GetPolydata());
-  }
+  void SetMesh(std::shared_ptr<MeshWrapper> mesh_, double geodesic_remesh_percent);
 
   std::shared_ptr<Mesh> GetSWMesh() const { return sw_mesh_; }
 
   void UpdateZeroCrossingPoint() override {}
 
  private:
   std::shared_ptr<MeshWrapper> mesh_wrapper_;
+  std::shared_ptr<MeshWrapper> geodesics_mesh_;
   std::shared_ptr<Mesh> sw_mesh_;
   PointType zero_crossing_point_;
 };