Improve mesh grooming robustness with integrity checking and repair capability

Libs/Groom/Groom.cpp

@@ -187,7 +187,7 @@ bool Groom::image_pipeline(std::shared_ptr<Subject> subject, size_t domain) {
 
   if (params.get_convert_to_mesh()) {
     Mesh mesh = image.toMesh(0.0);
-    run_mesh_pipeline(mesh, params);
+    run_mesh_pipeline(mesh, params, original);
     groomed_name = get_output_filename(original, DomainType::Mesh);
     // save the groomed mesh
     MeshUtils::threadSafeWriteMesh(groomed_name, mesh);
@@ -323,7 +323,7 @@ bool Groom::mesh_pipeline(std::shared_ptr<Subject> subject, size_t domain) {
   auto transform = vtkSmartPointer<vtkTransform>::New();
 
   if (!params.get_skip_grooming()) {
-    run_mesh_pipeline(mesh, params);
+    run_mesh_pipeline(mesh, params, original);
 
     // reflection
     if (params.get_reflect()) {
@@ -367,28 +367,17 @@ bool Groom::mesh_pipeline(std::shared_ptr<Subject> subject, size_t domain) {
 }
 
 //---------------------------------------------------------------------------
-bool Groom::run_mesh_pipeline(Mesh& mesh, GroomParameters params) {
-  // Extract largest component (should be first)
-  if (params.get_mesh_largest_component()) {
-    mesh.extractLargestComponent();
-    increment_progress();
-  }
+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()));
 
   if (params.get_fill_mesh_holes_tool()) {
     mesh.fillHoles();
+    check_and_fix_mesh(mesh, "fill_holes", filename);
     increment_progress();
   }
 
   if (params.get_remesh()) {
-    auto poly_data = mesh.getVTKMesh();
-    if (poly_data->GetNumberOfCells() == 0 || poly_data->GetCell(0)->GetNumberOfPoints() == 2) {
-      SW_DEBUG("Number of cells: {}", poly_data->GetNumberOfCells());
-      SW_DEBUG("Number of points in first cell: {}", poly_data->GetCell(0)->GetNumberOfPoints());
-      // write out to /tmp
-      std::string tmpname = "/tmp/bad_mesh.vtk";
-      MeshUtils::threadSafeWriteMesh(tmpname, mesh);
-      throw std::runtime_error("malformed mesh, mesh should be triangular");
-    }
     int total_vertices = mesh.getVTKMesh()->GetNumberOfPoints();
     int num_vertices = params.get_remesh_num_vertices();
     if (params.get_remesh_percent_mode()) {
@@ -397,6 +386,7 @@ bool Groom::run_mesh_pipeline(Mesh& mesh, GroomParameters params) {
     num_vertices = std::max<int>(num_vertices, 25);
     double gradation = clamp<double>(params.get_remesh_gradation(), 0.0, 2.0);
     mesh.remesh(num_vertices, gradation);
+    check_and_fix_mesh(mesh, "remesh", filename);
     increment_progress();
   }
 
@@ -406,6 +396,7 @@ bool Groom::run_mesh_pipeline(Mesh& mesh, GroomParameters params) {
     } else if (params.get_mesh_smoothing_method() == GroomParameters::GROOM_SMOOTH_VTK_WINDOWED_SINC_C) {
       mesh.smoothSinc(params.get_mesh_vtk_windowed_sinc_iterations(), params.get_mesh_vtk_windowed_sinc_passband());
     }
+    check_and_fix_mesh(mesh, "smooth", filename);
     increment_progress();
   }
   return true;
@@ -511,7 +502,6 @@ int Groom::get_total_ops() {
                     (project_->get_original_domain_types()[i] == DomainType::Image && params.get_convert_to_mesh());
 
     if (run_mesh) {
-      num_tools += params.get_mesh_largest_component() ? 1 : 0;
       num_tools += params.get_fill_mesh_holes_tool() ? 1 : 0;
       num_tools += params.get_mesh_smooth() ? 1 : 0;
       num_tools += params.get_remesh() ? 1 : 0;
@@ -1064,6 +1054,45 @@ void Groom::assign_transforms(std::vector<std::vector<double>> transforms, int d
   }
 }
 
+//---------------------------------------------------------------------------
+Mesh Groom::check_and_fix_mesh(Mesh& mesh, const std::string& step, const std::string& filename) {
+  auto issues = mesh.checkIntegrity();
+  if (!issues.empty()) {
+    SW_WARN("Mesh issues detected after step '{}', file '{}', issues: {}", step, filename, issues);
+    // try cleaning the mesh
+    SW_LOG("Cleaning mesh...");
+    mesh.clean();
+    issues = mesh.checkIntegrity();
+    if (!issues.empty()) {
+      SW_ERROR("Mesh issues remain after cleaning after step '{}', file '{}', issues: {}", step, filename, issues);
+      throw std::runtime_error(fmt::format("Error in mesh, step: {}, file: {}, issues: {}", step, filename, issues));
+    }
+  }
+  auto poly_data = mesh.getVTKMesh();
+
+  if (poly_data->GetNumberOfCells() == 0) {
+    throw std::runtime_error(fmt::format("Error in mesh, step: {}, file: {}, mesh has zero cells", step, filename));
+  }
+  if (poly_data->GetCell(0)->GetNumberOfPoints() == 2) {
+    // try cleaning the mesh
+    Mesh m(poly_data);
+    m.clean();
+    poly_data = m.getVTKMesh();
+
+    if (poly_data->GetNumberOfCells() == 0 || poly_data->GetCell(0)->GetNumberOfPoints() == 2) {
+      // still failed
+      SW_DEBUG("Number of cells: {}", poly_data->GetNumberOfCells());
+      SW_DEBUG("Number of points in first cell: {}", poly_data->GetCell(0)->GetNumberOfPoints());
+      // write out to /tmp
+      // std::string tmpname = "/tmp/bad_mesh.vtk";
+      // MeshUtils::threadSafeWriteMesh(tmpname, mesh);
+      throw std::runtime_error(
+          fmt::format("Error in mesh, step: {}, file: {}, mesh has invalid cells", step, filename));
+    }
+  }
+  return Mesh(poly_data);
+}
+
 //---------------------------------------------------------------------------
 vtkSmartPointer<vtkMatrix4x4> Groom::compute_landmark_transform(vtkSmartPointer<vtkPoints> source,
                                                                 vtkSmartPointer<vtkPoints> target) {

Libs/Groom/Groom.h

@@ -55,7 +55,7 @@ class Groom {
   //! Run the mesh based pipeline on a single subject
   bool mesh_pipeline(std::shared_ptr<Subject> subject, size_t domain);
 
-  bool run_mesh_pipeline(Mesh& mesh, GroomParameters params);
+  bool run_mesh_pipeline(Mesh& mesh, GroomParameters params, const std::string& filename);
 
   //! Run the contour based pipeline on a single subject
   bool contour_pipeline(std::shared_ptr<Subject> subject, size_t domain);
@@ -73,6 +73,8 @@ class Groom {
 
   void assign_transforms(std::vector<std::vector<double>> transforms, int domain, bool global = false);
 
+  Mesh check_and_fix_mesh(Mesh& mesh, const std::string& step, const std::string& filename);
+
   static std::vector<std::vector<double>> get_icp_transforms(const std::vector<Mesh> meshes, Mesh reference);
   static std::vector<std::vector<double>> get_landmark_transforms(
       const std::vector<vtkSmartPointer<vtkPoints>> landmarks, size_t reference);

Libs/Groom/GroomParameters.cpp

@@ -23,7 +23,6 @@ const std::string ISO_SPACING = "iso_spacing";
 const std::string SPACING = "spacing";
 const std::string CONVERT_MESH = "convert_to_mesh";
 const std::string FILL_MESH_HOLES = "fill_mesh_holes";
-const std::string MESH_LARGEST_COMPONENT = "mesh_largest_component";
 const std::string FILL_HOLES = "fill_holes";
 const std::string ISOLATE = "isolate";
 const std::string PAD = "pad";
@@ -78,7 +77,6 @@ const std::vector<double> spacing{0, 0, 0};
 const bool convert_mesh = false;
 const bool fill_holes = true;
 const bool fill_holes_mesh = false;
-const bool mesh_largest_component = true;
 const bool isolate = true;
 const bool pad = true;
 const int pad_value = 10;
@@ -158,7 +156,6 @@ GroomParameters::GroomParameters(ProjectHandle project, std::string domain_name)
       Keys::SPACING,
       Keys::CONVERT_MESH,
       Keys::FILL_MESH_HOLES,
-      Keys::MESH_LARGEST_COMPONENT,
       Keys::FILL_HOLES,
       Keys::ISOLATE,
       Keys::PAD,
@@ -196,12 +193,20 @@ GroomParameters::GroomParameters(ProjectHandle project, std::string domain_name)
       Keys::SHARED_BOUNDARIES,
   };
 
+  // Deprecated parameters that should be silently removed
+  std::vector<std::string> deprecated_params = {
+      "mesh_largest_component",  // Now handled by repair_mesh at start of pipeline
+  };
+
   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);
+      // Only warn if not a deprecated parameter
+      if (std::find(deprecated_params.begin(), deprecated_params.end(), param.first) == deprecated_params.end()) {
+        SW_WARN("Unknown Grooming parameter: " + param.first);
+      }
       to_remove.push_back(param.first);
     }
   }
@@ -236,14 +241,6 @@ bool GroomParameters::get_fill_mesh_holes_tool() {
 //---------------------------------------------------------------------------
 void GroomParameters::set_fill_mesh_holes_tool(bool value) { params_.set(Keys::FILL_MESH_HOLES, value); }
 
-//---------------------------------------------------------------------------
-bool GroomParameters::get_mesh_largest_component() {
-  return params_.get(Keys::MESH_LARGEST_COMPONENT, Defaults::mesh_largest_component);
-}
-
-//---------------------------------------------------------------------------
-void GroomParameters::set_mesh_largest_component(bool value) { params_.set(Keys::MESH_LARGEST_COMPONENT, value); }
-
 //---------------------------------------------------------------------------
 bool GroomParameters::get_auto_pad_tool() { return params_.get(Keys::PAD, Defaults::pad); }
 

Libs/Groom/GroomParameters.h

@@ -61,9 +61,6 @@ class GroomParameters {
   bool get_fill_mesh_holes_tool();
   void set_fill_mesh_holes_tool(bool value);
 
-  bool get_mesh_largest_component();
-  void set_mesh_largest_component(bool value);
-
   bool get_auto_pad_tool();
   void set_auto_pad_tool(bool value);
 

Libs/Mesh/Mesh.cpp

@@ -341,6 +341,9 @@ Mesh& Mesh::remesh(int numVertices, double adaptivity) {
   // Restore std::cout
   //  std::cout.clear();
 
+  clean();
+  checkIntegrity();
+
   this->poly_data_ = remesh->GetOutput();
 
   // must regenerate normals after smoothing
@@ -415,14 +418,12 @@ Mesh& Mesh::fillHoles(double hole_size) {
   filter->Update();
   this->poly_data_ = filter->GetOutput();
 
-  auto origNormal = poly_data_->GetPointData()->GetNormals();
+  // vtkFillHolesFilter can create zero-area triangles, remove them
+  this->poly_data_ = MeshUtils::remove_zero_area_triangles(this->poly_data_);
 
-  // Make the triangle window order consistent
+  // Make the triangle window order consistent and recompute normals
   computeNormals();
 
-  // Restore the original normals
-  poly_data_->GetPointData()->SetNormals(origNormal);
-
   this->invalidateLocators();
   return *this;
 }
@@ -728,6 +729,11 @@ Mesh& Mesh::clipClosedSurface(const Plane plane) {
 }
 
 Mesh& Mesh::computeNormals() {
+  // Validate mesh first
+  if (!poly_data_ || poly_data_->GetNumberOfPoints() == 0 || poly_data_->GetNumberOfCells() == 0) {
+    return *this;
+  }
+
   // Remove existing normals first
   poly_data_->GetPointData()->SetNormals(nullptr);
   poly_data_->GetCellData()->SetNormals(nullptr);
@@ -1758,6 +1764,133 @@ void Mesh::interpolate_scalars_to_mesh(std::string name, Eigen::VectorXd positio
   setField(name, scalars, Mesh::FieldType::Point);
 }
 
+std::string Mesh::checkIntegrity() const {
+  std::ostringstream issues;
+
+  if (!poly_data_) {
+    return "poly_data_ is null";
+  }
+
+  const vtkIdType numPoints = poly_data_->GetNumberOfPoints();
+  const vtkIdType numCells = poly_data_->GetNumberOfCells();
+
+  if (numPoints == 0) {
+    issues << "Mesh has no points\n";
+  }
+  if (numCells == 0) {
+    issues << "Mesh has no cells\n";
+  }
+  if (numPoints == 0 || numCells == 0) {
+    return issues.str();
+  }
+
+  // Check for NaN/Inf coordinates
+  vtkPoints* points = poly_data_->GetPoints();
+  vtkIdType nanCount = 0;
+  vtkIdType infCount = 0;
+  for (vtkIdType i = 0; i < numPoints; i++) {
+    double p[3];
+    points->GetPoint(i, p);
+    for (int j = 0; j < 3; j++) {
+      if (std::isnan(p[j])) nanCount++;
+      if (std::isinf(p[j])) infCount++;
+    }
+  }
+  if (nanCount > 0) {
+    issues << "Found " << nanCount << " NaN coordinates\n";
+  }
+  if (infCount > 0) {
+    issues << "Found " << infCount << " Inf coordinates\n";
+  }
+
+  // Check cell validity
+  vtkCellArray* polys = poly_data_->GetPolys();
+  vtkCellArray* strips = poly_data_->GetStrips();
+  vtkCellArray* lines = poly_data_->GetLines();
+  vtkCellArray* verts = poly_data_->GetVerts();
+
+  vtkIdType degenerateCells = 0;
+  vtkIdType invalidIndexCells = 0;
+  vtkIdType duplicatePointCells = 0;
+  vtkIdType zerAreaCells = 0;
+
+  for (vtkIdType cellId = 0; cellId < numCells; cellId++) {
+    vtkCell* cell = poly_data_->GetCell(cellId);
+    if (!cell) {
+      issues << "Null cell at index " << cellId << "\n";
+      continue;
+    }
+
+    vtkIdType npts = cell->GetNumberOfPoints();
+
+    // Check for degenerate cells (less than 3 points for a polygon)
+    if (cell->GetCellDimension() == 2 && npts < 3) {
+      degenerateCells++;
+      continue;
+    }
+
+    // Check for invalid point indices and duplicate points within cell
+    std::set<vtkIdType> pointsInCell;
+    bool hasInvalidIndex = false;
+    bool hasDuplicate = false;
+
+    for (vtkIdType i = 0; i < npts; i++) {
+      vtkIdType ptId = cell->GetPointId(i);
+      if (ptId < 0 || ptId >= numPoints) {
+        hasInvalidIndex = true;
+      }
+      if (pointsInCell.count(ptId) > 0) {
+        hasDuplicate = true;
+      }
+      pointsInCell.insert(ptId);
+    }
+
+    if (hasInvalidIndex) invalidIndexCells++;
+    if (hasDuplicate) duplicatePointCells++;
+
+    // Check for zero-area triangles
+    if (npts == 3 && !hasInvalidIndex) {
+      double p0[3], p1[3], p2[3];
+      points->GetPoint(cell->GetPointId(0), p0);
+      points->GetPoint(cell->GetPointId(1), p1);
+      points->GetPoint(cell->GetPointId(2), p2);
+      double area = vtkTriangle::TriangleArea(p0, p1, p2);
+      if (area < 1e-20) {
+        zerAreaCells++;
+      }
+    }
+  }
+
+  if (degenerateCells > 0) {
+    issues << "Found " << degenerateCells << " degenerate cells (< 3 points)\n";
+  }
+  if (invalidIndexCells > 0) {
+    issues << "Found " << invalidIndexCells << " cells with invalid point indices\n";
+  }
+  if (duplicatePointCells > 0) {
+    issues << "Found " << duplicatePointCells << " cells with duplicate point references\n";
+  }
+  if (zerAreaCells > 0) {
+    issues << "Found " << zerAreaCells << " zero-area triangles\n";
+  }
+
+  // Check for non-manifold edges (optional, can be slow)
+  auto featureEdges = vtkSmartPointer<vtkFeatureEdges>::New();
+  featureEdges->SetInputData(poly_data_);
+  featureEdges->BoundaryEdgesOff();
+  featureEdges->FeatureEdgesOff();
+  featureEdges->ManifoldEdgesOff();
+  featureEdges->NonManifoldEdgesOn();
+  featureEdges->Update();
+
+  vtkIdType nonManifoldEdges = featureEdges->GetOutput()->GetNumberOfCells();
+  if (nonManifoldEdges > 0) {
+    issues << "Found " << nonManifoldEdges << " non-manifold edges\n";
+  }
+
+  return issues.str();
+}
+
 double Mesh::getFFCValue(Eigen::Vector3d query) const {
   Point3 point(query.data());
   return interpolateFieldAtPoint("value", point);

Libs/Mesh/Mesh.h

@@ -302,6 +302,8 @@ class Mesh {
   //! given name
   void interpolate_scalars_to_mesh(std::string name, Eigen::VectorXd positions, Eigen::VectorXd scalar_values);
 
+  std::string checkIntegrity() const;
+
  private:
   friend struct SharedCommandData;
   Mesh()