Fix shared boundary donut: skip extract-largest-component for shared boundary domains, add empty mesh crash guards

Author: akenmorris

Libs/Groom/Groom.cpp

@@ -385,8 +385,10 @@ bool Groom::mesh_pipeline(std::shared_ptr<Subject> subject, size_t domain) {
 
 //---------------------------------------------------------------------------
 bool Groom::run_mesh_pipeline(Mesh& mesh, GroomParameters params, const std::string& filename) {
-  // Repair mesh: triangulate, extract largest component, clean, fix non-manifold, remove zero-area triangles
-  mesh = Mesh(MeshUtils::repair_mesh(mesh.getVTKMesh()));
+  // Repair mesh: triangulate, clean, fix non-manifold, remove zero-area triangles
+  // Skip extract-largest-component for shared boundary domains to avoid removing fragments at the cut surface
+  bool extract_largest = !params.is_shared_boundary_domain();
+  mesh = Mesh(MeshUtils::repair_mesh(mesh.getVTKMesh(), extract_largest));
 
   if (params.get_fill_mesh_holes_tool()) {
     mesh.fillHoles();

Libs/Groom/GroomParameters.cpp

@@ -544,6 +544,19 @@ bool GroomParameters::get_shared_boundaries_enabled() {
 //---------------------------------------------------------------------------
 void GroomParameters::set_shared_boundaries_enabled(bool enabled) { params_.set(Keys::SHARED_BOUNDARY, enabled); }
 
+//---------------------------------------------------------------------------
+bool GroomParameters::is_shared_boundary_domain() {
+  if (!get_shared_boundaries_enabled()) {
+    return false;
+  }
+  for (const auto& boundary : get_shared_boundaries()) {
+    if (boundary.first_domain == domain_name_ || boundary.second_domain == domain_name_) {
+      return true;
+    }
+  }
+  return false;
+}
+
 //---------------------------------------------------------------------------
 std::vector<GroomParameters::SharedBoundary> GroomParameters::get_shared_boundaries() {
   std::vector<SharedBoundary> boundaries;

Libs/Groom/GroomParameters.h

@@ -151,6 +151,9 @@ class GroomParameters {
   bool get_shared_boundaries_enabled();
   void set_shared_boundaries_enabled(bool enabled);
 
+  //! Check if the current domain participates in any shared boundary
+  bool is_shared_boundary_domain();
+
   std::vector<SharedBoundary> get_shared_boundaries();
   void set_shared_boundaries(const std::vector<SharedBoundary>& boundaries);
   void add_shared_boundary(const std::string& first_domain, const std::string& second_domain, double tolerance);

Libs/Mesh/MeshUtils.cpp

@@ -389,6 +389,12 @@ static std::tuple<Eigen::MatrixXd, Eigen::MatrixXi, Eigen::MatrixXd, Eigen::Matr
   Eigen::MatrixXi out_F;
   Eigen::MatrixXd rem_V;
   Eigen::MatrixXi rem_F;
+
+  // If either mesh is empty, there can be no shared surface
+  if (is_empty(src_V, src_F) || is_empty(other_V, other_F)) {
+    return std::make_tuple(out_V, out_F, src_V, src_F);
+  }
+
   igl::AABB<Eigen::MatrixXd, 3> tree;
   tree.init(other_V, other_F);
 
@@ -482,6 +488,10 @@ std::array<Mesh, 3> MeshUtils::shared_boundary_extractor(const Mesh& mesh_l, con
   V_r = mesh_r.points();
   F_r = mesh_r.faces();
 
+  if (is_empty(V_l, F_l) || is_empty(V_r, F_r)) {
+    throw std::runtime_error("Input mesh is empty. Cannot extract shared boundary from empty meshes");
+  }
+
   Eigen::MatrixXd shared_V_l, shared_V_r, rem_V_l, rem_V_r;
   Eigen::MatrixXi shared_F_l, shared_F_r, rem_F_l, rem_F_r;
   std::tie(shared_V_l, shared_F_l, rem_V_l, rem_F_l) = find_shared_surface(V_l, F_l, V_r, F_r, tol);
@@ -1073,18 +1083,23 @@ vtkSmartPointer<vtkPolyData> MeshUtils::recreate_mesh(vtkSmartPointer<vtkPolyDat
 }
 
 //---------------------------------------------------------------------------
-vtkSmartPointer<vtkPolyData> MeshUtils::repair_mesh(vtkSmartPointer<vtkPolyData> mesh) {
+vtkSmartPointer<vtkPolyData> MeshUtils::repair_mesh(vtkSmartPointer<vtkPolyData> mesh, bool extract_largest) {
   auto triangle_filter = vtkSmartPointer<vtkTriangleFilter>::New();
   triangle_filter->SetInputData(mesh);
   triangle_filter->PassLinesOff();
   triangle_filter->Update();
 
-  auto connectivity = vtkSmartPointer<vtkPolyDataConnectivityFilter>::New();
-  connectivity->SetInputConnection(triangle_filter->GetOutputPort());
-  connectivity->SetExtractionModeToLargestRegion();
-  connectivity->Update();
+  vtkSmartPointer<vtkPolyData> triangulated = triangle_filter->GetOutput();
+
+  if (extract_largest) {
+    auto connectivity = vtkSmartPointer<vtkPolyDataConnectivityFilter>::New();
+    connectivity->SetInputData(triangulated);
+    connectivity->SetExtractionModeToLargestRegion();
+    connectivity->Update();
+    triangulated = connectivity->GetOutput();
+  }
 
-  auto cleaned = MeshUtils::clean_mesh(connectivity->GetOutput());
+  auto cleaned = MeshUtils::clean_mesh(triangulated);
 
   auto fixed = Mesh(cleaned).fixNonManifold();
 

Libs/Mesh/MeshUtils.h

@@ -86,7 +86,7 @@ class MeshUtils {
     /// Recreate mesh, dropping deleted cells
     static vtkSmartPointer<vtkPolyData> recreate_mesh(vtkSmartPointer<vtkPolyData> mesh);
 
-    /// Repair mesh: triangulate, extract largest component, clean, fix non-manifold, remove zero-area triangles
-    static vtkSmartPointer<vtkPolyData> repair_mesh(vtkSmartPointer<vtkPolyData> mesh);
+    /// Repair mesh: triangulate, optionally extract largest component, clean, fix non-manifold, remove zero-area triangles
+    static vtkSmartPointer<vtkPolyData> repair_mesh(vtkSmartPointer<vtkPolyData> mesh, bool extract_largest = true);
 };
 } // namespace shapeworks

Libs/Optimize/OptimizeParameters.cpp

@@ -328,6 +328,9 @@ std::vector<std::vector<itk::Point<double>>> OptimizeParameters::get_initial_poi
     for (auto s : subjects) {
       if (s->is_fixed()) {
         count++;
+        if (d >= s->get_world_particle_filenames().size()) {
+          throw std::runtime_error("Subject " + s->get_display_name() + " does not have enough world particle files");
+        }
         // 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);