From 6b41a203d0187179d302e49608bd763a5bc22cc2 Mon Sep 17 00:00:00 2001 From: Alan Morris Date: Thu, 18 Dec 2025 13:57:32 -0700 Subject: [PATCH 1/6] Add mesh integrity checking throughout grooming pipeline --- Libs/Groom/Groom.cpp | 60 +++++++++++++---- Libs/Groom/Groom.h | 4 +- Libs/Mesh/Mesh.cpp | 149 ++++++++++++++++++++++++++++++++++++++++++ Libs/Mesh/Mesh.h | 2 + Libs/Mesh/MeshUtils.h | 2 + 5 files changed, 204 insertions(+), 13 deletions(-) diff --git a/Libs/Groom/Groom.cpp b/Libs/Groom/Groom.cpp index 87a1f097fb..fea1d5c6e1 100644 --- a/Libs/Groom/Groom.cpp +++ b/Libs/Groom/Groom.cpp @@ -187,7 +187,7 @@ bool Groom::image_pipeline(std::shared_ptr 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, size_t domain) { auto transform = vtkSmartPointer::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,23 @@ bool Groom::mesh_pipeline(std::shared_ptr subject, size_t domain) { } //--------------------------------------------------------------------------- -bool Groom::run_mesh_pipeline(Mesh& mesh, GroomParameters params) { +bool Groom::run_mesh_pipeline(Mesh& mesh, GroomParameters params, const std::string& filename) { + check_and_fix_mesh(mesh, "initial", filename); + // Extract largest component (should be first) if (params.get_mesh_largest_component()) { mesh.extractLargestComponent(); + check_and_fix_mesh(mesh, "largest_component", filename); increment_progress(); } 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 +392,7 @@ bool Groom::run_mesh_pipeline(Mesh& mesh, GroomParameters params) { num_vertices = std::max(num_vertices, 25); double gradation = clamp(params.get_remesh_gradation(), 0.0, 2.0); mesh.remesh(num_vertices, gradation); + check_and_fix_mesh(mesh, "remesh", filename); increment_progress(); } @@ -406,6 +402,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; @@ -1064,6 +1061,45 @@ void Groom::assign_transforms(std::vector> 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 at 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 at 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 Groom::compute_landmark_transform(vtkSmartPointer source, vtkSmartPointer target) { diff --git a/Libs/Groom/Groom.h b/Libs/Groom/Groom.h index cec066375a..ff669a2c3a 100644 --- a/Libs/Groom/Groom.h +++ b/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, 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, size_t domain); @@ -73,6 +73,8 @@ class Groom { void assign_transforms(std::vector> transforms, int domain, bool global = false); + Mesh check_and_fix_mesh(Mesh& mesh, const std::string& step, const std::string& filename); + static std::vector> get_icp_transforms(const std::vector meshes, Mesh reference); static std::vector> get_landmark_transforms( const std::vector> landmarks, size_t reference); diff --git a/Libs/Mesh/Mesh.cpp b/Libs/Mesh/Mesh.cpp index b1d05de234..d105f3eac3 100644 --- a/Libs/Mesh/Mesh.cpp +++ b/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 @@ -728,10 +731,28 @@ Mesh& Mesh::clipClosedSurface(const Plane plane) { } Mesh& Mesh::computeNormals() { + + std::string issues = checkIntegrity(); + if (!issues.empty()) { + std::cerr << "Mesh integrity check failed:\n" << issues << std::endl; + // Either return early or try to fix + return *this; + } + // Remove existing normals first poly_data_->GetPointData()->SetNormals(nullptr); poly_data_->GetCellData()->SetNormals(nullptr); + // Validate mesh first + if (!poly_data_ || poly_data_->GetNumberOfPoints() == 0) { + std::cerr << "Warning: Empty mesh in computeNormals()" << std::endl; + return *this; + } + if (poly_data_->GetNumberOfCells() == 0) { + std::cerr << "Warning: Mesh has no cells in computeNormals()" << std::endl; + return *this; + } + auto normal = vtkSmartPointer::New(); normal->SetInputData(this->poly_data_); @@ -1758,6 +1779,134 @@ 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 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::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); diff --git a/Libs/Mesh/Mesh.h b/Libs/Mesh/Mesh.h index a9008ca71f..5040bbb988 100644 --- a/Libs/Mesh/Mesh.h +++ b/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() diff --git a/Libs/Mesh/MeshUtils.h b/Libs/Mesh/MeshUtils.h index ab991ae8d6..270ad99cf6 100644 --- a/Libs/Mesh/MeshUtils.h +++ b/Libs/Mesh/MeshUtils.h @@ -76,5 +76,7 @@ class MeshUtils { double pt3[3], double pt1[3], double pt2[3]); + + }; } // namespace shapeworks From 7dca439cf74bed38fd0af0fe75b5420118cf5f7b Mon Sep 17 00:00:00 2001 From: Alan Morris Date: Thu, 18 Dec 2025 14:09:31 -0700 Subject: [PATCH 2/6] Moved MeshWarper basic mesh routines to MeshUtils --- Libs/Mesh/MeshUtils.cpp | 146 ++++++++++++++++++++++++++++++ Libs/Mesh/MeshUtils.h | 10 +++ Libs/Mesh/MeshWarper.cpp | 148 +------------------------------ Libs/Mesh/MeshWarper.h | 12 --- Libs/Python/ShapeworksPython.cpp | 2 +- 5 files changed, 161 insertions(+), 157 deletions(-) diff --git a/Libs/Mesh/MeshUtils.cpp b/Libs/Mesh/MeshUtils.cpp index f7e35ed87e..2468230ecb 100644 --- a/Libs/Mesh/MeshUtils.cpp +++ b/Libs/Mesh/MeshUtils.cpp @@ -4,7 +4,11 @@ #include #include #include +#include #include +#include +#include +#include #include #include #include @@ -942,4 +946,146 @@ int MeshUtils::evaluate_triangle_position(const double x[3], double closestPoint return 0; } } + +//--------------------------------------------------------------------------- +vtkSmartPointer MeshUtils::clean_mesh(vtkSmartPointer mesh) { + vtkSmartPointer clean = vtkSmartPointer::New(); + clean->ConvertPolysToLinesOff(); + clean->ConvertLinesToPointsOff(); + clean->ConvertStripsToPolysOff(); + clean->PointMergingOn(); + clean->SetInputData(mesh); + clean->Update(); + return clean->GetOutput(); +} + +//--------------------------------------------------------------------------- +vtkSmartPointer MeshUtils::remove_zero_area_triangles(vtkSmartPointer mesh) { + // Create a new polydata to store the filtered triangles + auto filtered_mesh = vtkSmartPointer::New(); + + // Get the points from the input mesh + vtkSmartPointer points = mesh->GetPoints(); + filtered_mesh->SetPoints(points); + + // Create a new cell array to store triangles + auto new_triangles = vtkSmartPointer::New(); + + // Iterate through all cells + auto cell_point_ids = vtkSmartPointer::New(); + + // First pass: find max area to establish scale + double max_area = 0.0; + for (vtkIdType cellId = 0; cellId < mesh->GetNumberOfCells(); cellId++) { + mesh->GetCellPoints(cellId, cell_point_ids); + if (cell_point_ids->GetNumberOfIds() == 3) { + double p0[3], p1[3], p2[3]; + points->GetPoint(cell_point_ids->GetId(0), p0); + points->GetPoint(cell_point_ids->GetId(1), p1); + points->GetPoint(cell_point_ids->GetId(2), p2); + double area = vtkTriangle::TriangleArea(p0, p1, p2); + if (area > max_area) { + max_area = area; + } + } + } + + // Use relative epsilon based on max area, with absolute floor + const double RELATIVE_EPSILON = 1e-6; + const double ABSOLUTE_FLOOR = 1e-15; + const double epsilon = std::max(max_area * RELATIVE_EPSILON, ABSOLUTE_FLOOR); + + // Second pass: filter triangles + for (vtkIdType cellId = 0; cellId < mesh->GetNumberOfCells(); cellId++) { + mesh->GetCellPoints(cellId, cell_point_ids); + + // Only process triangles + if (cell_point_ids->GetNumberOfIds() == 3) { + // Get the coordinates of the three vertices + double p0[3], p1[3], p2[3]; + points->GetPoint(cell_point_ids->GetId(0), p0); + points->GetPoint(cell_point_ids->GetId(1), p1); + points->GetPoint(cell_point_ids->GetId(2), p2); + + // Calculate the area of the triangle using VTK's built-in function + double area = vtkTriangle::TriangleArea(p0, p1, p2); + + if (area > epsilon) { + // Add this triangle to the new cell array + new_triangles->InsertNextCell(cell_point_ids); + } + } + } + + // Set the triangles in the filtered mesh + filtered_mesh->SetPolys(new_triangles); + + // Copy the point data + filtered_mesh->GetPointData()->ShallowCopy(mesh->GetPointData()); + + // Use vtkCleanPolyData to remove unused points + auto cleaner = vtkSmartPointer::New(); + cleaner->SetInputData(filtered_mesh); + cleaner->PointMergingOff(); // Don't merge points, just remove unused ones + cleaner->Update(); + + return cleaner->GetOutput(); +} + +//--------------------------------------------------------------------------- +vtkSmartPointer MeshUtils::recreate_mesh(vtkSmartPointer mesh) { + vtkSmartPointer poly_data = vtkSmartPointer::New(); + + vtkNew points; + vtkNew polys; + + // copy points + for (vtkIdType i = 0; i < mesh->GetNumberOfPoints(); i++) { + points->InsertNextPoint(mesh->GetPoint(i)); + } + + // copy triangles + for (vtkIdType i = 0; i < mesh->GetNumberOfCells(); i++) { + vtkCell* cell = mesh->GetCell(i); + + if (cell->GetCellType() != VTK_EMPTY_CELL) { // VTK_EMPTY_CELL means it was deleted + // create an array of vtkIdType + vtkIdType* pts = new vtkIdType[cell->GetNumberOfPoints()]; + for (vtkIdType j = 0; j < cell->GetNumberOfPoints(); j++) { + pts[j] = cell->GetPointId(j); + } + polys->InsertNextCell(cell->GetNumberOfPoints(), pts); + delete[] pts; + } + } + + poly_data->SetPoints(points); + poly_data->SetPolys(polys); + return poly_data; +} + +//--------------------------------------------------------------------------- +vtkSmartPointer MeshUtils::repair_mesh(vtkSmartPointer mesh) { + auto triangle_filter = vtkSmartPointer::New(); + triangle_filter->SetInputData(mesh); + triangle_filter->PassLinesOff(); + triangle_filter->Update(); + + auto connectivity = vtkSmartPointer::New(); + connectivity->SetInputConnection(triangle_filter->GetOutputPort()); + connectivity->SetExtractionModeToLargestRegion(); + connectivity->Update(); + + auto cleaned = MeshUtils::clean_mesh(connectivity->GetOutput()); + + auto fixed = Mesh(cleaned).fixNonManifold(); + + // remove deleted triangles and clean up + auto recreated = MeshUtils::recreate_mesh(fixed.getVTKMesh()); + + auto final = MeshUtils::remove_zero_area_triangles(recreated); + + return final; +} + } // namespace shapeworks diff --git a/Libs/Mesh/MeshUtils.h b/Libs/Mesh/MeshUtils.h index 270ad99cf6..7cf1b01eb0 100644 --- a/Libs/Mesh/MeshUtils.h +++ b/Libs/Mesh/MeshUtils.h @@ -77,6 +77,16 @@ class MeshUtils { double pt1[3], double pt2[3]); + /// Clean mesh (merge points and remove degenerate cells) + static vtkSmartPointer clean_mesh(vtkSmartPointer mesh); + /// Remove zero area triangles from a mesh + static vtkSmartPointer remove_zero_area_triangles(vtkSmartPointer mesh); + + /// Recreate mesh, dropping deleted cells + static vtkSmartPointer recreate_mesh(vtkSmartPointer mesh); + + /// Repair mesh: triangulate, extract largest component, clean, fix non-manifold, remove zero-area triangles + static vtkSmartPointer repair_mesh(vtkSmartPointer mesh); }; } // namespace shapeworks diff --git a/Libs/Mesh/MeshWarper.cpp b/Libs/Mesh/MeshWarper.cpp index 6ce79d645c..0a2d18ba8d 100644 --- a/Libs/Mesh/MeshWarper.cpp +++ b/Libs/Mesh/MeshWarper.cpp @@ -1,5 +1,6 @@ #include #include +#include #include #include #include @@ -257,7 +258,7 @@ void MeshWarper::add_particle_vertices(Eigen::MatrixXd& vertices) { } // remove deleted triangles and clean up - reference_mesh_ = MeshWarper::recreate_mesh(reference_mesh_); + reference_mesh_ = MeshUtils::recreate_mesh(reference_mesh_); } } @@ -369,147 +370,6 @@ void MeshWarper::split_cell_on_edge(int cell_id, int new_vertex, int v0, int v1, reference_mesh_->DeleteCell(cell_id); } -//--------------------------------------------------------------------------- -vtkSmartPointer MeshWarper::prep_mesh(vtkSmartPointer mesh) { - auto triangle_filter = vtkSmartPointer::New(); - triangle_filter->SetInputData(mesh); - triangle_filter->PassLinesOff(); - triangle_filter->Update(); - - auto connectivity = vtkSmartPointer::New(); - connectivity->SetInputConnection(triangle_filter->GetOutputPort()); - connectivity->SetExtractionModeToLargestRegion(); - connectivity->Update(); - - auto cleaned = MeshWarper::clean_mesh(connectivity->GetOutput()); - - auto fixed = Mesh(cleaned).fixNonManifold(); - - // remove deleted triangles and clean up - auto recreated = MeshWarper::recreate_mesh(fixed.getVTKMesh()); - - auto final = MeshWarper::remove_zero_area_triangles(recreated); - - return final; -} - -//--------------------------------------------------------------------------- -vtkSmartPointer MeshWarper::clean_mesh(vtkSmartPointer mesh) { - vtkSmartPointer clean = vtkSmartPointer::New(); - clean->ConvertPolysToLinesOff(); - clean->ConvertLinesToPointsOff(); - clean->ConvertStripsToPolysOff(); - clean->PointMergingOn(); - clean->SetInputData(mesh); - clean->Update(); - return clean->GetOutput(); -} - -//--------------------------------------------------------------------------- -vtkSmartPointer MeshWarper::remove_zero_area_triangles(vtkSmartPointer mesh) { - // Create a new polydata to store the filtered triangles - auto filtered_mesh = vtkSmartPointer::New(); - - // Get the points from the input mesh - vtkSmartPointer points = mesh->GetPoints(); - filtered_mesh->SetPoints(points); - - // Create a new cell array to store triangles - auto new_triangles = vtkSmartPointer::New(); - - // Iterate through all cells - auto cell_point_ids = vtkSmartPointer::New(); - - // First pass: find max area to establish scale - double max_area = 0.0; - for (vtkIdType cellId = 0; cellId < mesh->GetNumberOfCells(); cellId++) { - mesh->GetCellPoints(cellId, cell_point_ids); - if (cell_point_ids->GetNumberOfIds() == 3) { - double p0[3], p1[3], p2[3]; - points->GetPoint(cell_point_ids->GetId(0), p0); - points->GetPoint(cell_point_ids->GetId(1), p1); - points->GetPoint(cell_point_ids->GetId(2), p2); - double area = vtkTriangle::TriangleArea(p0, p1, p2); - if (area > max_area) { - max_area = area; - } - } - } - - // Use relative epsilon based on max area, with absolute floor - const double RELATIVE_EPSILON = 1e-6; - const double ABSOLUTE_FLOOR = 1e-15; - const double epsilon = std::max(max_area * RELATIVE_EPSILON, ABSOLUTE_FLOOR); - - // Second pass: filter triangles - for (vtkIdType cellId = 0; cellId < mesh->GetNumberOfCells(); cellId++) { - mesh->GetCellPoints(cellId, cell_point_ids); - - // Only process triangles - if (cell_point_ids->GetNumberOfIds() == 3) { - // Get the coordinates of the three vertices - double p0[3], p1[3], p2[3]; - points->GetPoint(cell_point_ids->GetId(0), p0); - points->GetPoint(cell_point_ids->GetId(1), p1); - points->GetPoint(cell_point_ids->GetId(2), p2); - - // Calculate the area of the triangle using VTK's built-in function - double area = vtkTriangle::TriangleArea(p0, p1, p2); - - if (area > epsilon) { - // Add this triangle to the new cell array - new_triangles->InsertNextCell(cell_point_ids); - } - } - } - - // Set the triangles in the filtered mesh - filtered_mesh->SetPolys(new_triangles); - - // Copy the point data - filtered_mesh->GetPointData()->ShallowCopy(mesh->GetPointData()); - - // Use vtkCleanPolyData to remove unused points - auto cleaner = vtkSmartPointer::New(); - cleaner->SetInputData(filtered_mesh); - cleaner->PointMergingOff(); // Don't merge points, just remove unused ones - cleaner->Update(); - - return cleaner->GetOutput(); -} - -//--------------------------------------------------------------------------- -vtkSmartPointer MeshWarper::recreate_mesh(vtkSmartPointer mesh) { - vtkSmartPointer poly_data = vtkSmartPointer::New(); - - vtkNew points; - vtkNew polys; - - // copy points - for (vtkIdType i = 0; i < mesh->GetNumberOfPoints(); i++) { - points->InsertNextPoint(mesh->GetPoint(i)); - } - - // copy triangles - for (vtkIdType i = 0; i < mesh->GetNumberOfCells(); i++) { - vtkCell* cell = mesh->GetCell(i); - - if (cell->GetCellType() != VTK_EMPTY_CELL) { // VTK_EMPTY_CELL means it was deleted - // create an array of vtkIdType - vtkIdType* pts = new vtkIdType[cell->GetNumberOfPoints()]; - for (vtkIdType j = 0; j < cell->GetNumberOfPoints(); j++) { - pts[j] = cell->GetPointId(j); - } - polys->InsertNextCell(cell->GetNumberOfPoints(), pts); - delete[] pts; - } - } - - poly_data->SetPoints(points); - poly_data->SetPolys(polys); - return poly_data; -} - //--------------------------------------------------------------------------- bool MeshWarper::generate_warp() { if (is_contour_) { @@ -523,12 +383,12 @@ bool MeshWarper::generate_warp() { for (int attempt = 0; attempt < MAX_ATTEMPTS; ++attempt) { // Prepare mesh (remesh on retry) if (attempt == 0) { - reference_mesh_ = MeshWarper::prep_mesh(incoming_reference_mesh_); + reference_mesh_ = MeshUtils::repair_mesh(incoming_reference_mesh_); } else { SW_DEBUG("Initial igl::biharmonic failed, attempting again with remeshed reference mesh"); Mesh mesh(reference_mesh_); mesh.remeshPercent(0.75); - reference_mesh_ = MeshWarper::prep_mesh(mesh.getVTKMesh()); + reference_mesh_ = MeshUtils::repair_mesh(mesh.getVTKMesh()); } // Prep points diff --git a/Libs/Mesh/MeshWarper.h b/Libs/Mesh/MeshWarper.h index 500e6deb71..610cb2e95d 100644 --- a/Libs/Mesh/MeshWarper.h +++ b/Libs/Mesh/MeshWarper.h @@ -66,9 +66,6 @@ class MeshWarper { //! Return the reference particles const Eigen::MatrixXd& get_reference_particles() const { return this->reference_particles_; } - //! Prep incoming mesh - static vtkSmartPointer prep_mesh(vtkSmartPointer mesh); - protected: //! For overriding to handle progress updates virtual void update_progress(float p) {} @@ -93,15 +90,6 @@ class MeshWarper { //! Construct the map from landmarks to vertex ids bool find_landmarks_vertices_on_ref_mesh(); - //! Clean mesh (remove deleted) - static vtkSmartPointer clean_mesh(vtkSmartPointer mesh); - - //! Remove zero area triangles from a mesh - static vtkSmartPointer remove_zero_area_triangles(vtkSmartPointer mesh); - - //! Recreate mesh, dropping deleted cells - static vtkSmartPointer recreate_mesh(vtkSmartPointer mesh); - //! Generate the warp matrix bool generate_warp_matrix(Eigen::MatrixXd target_vertices, Eigen::MatrixXi target_faces, const Eigen::MatrixXd& references_vertices, Eigen::MatrixXd& warp); diff --git a/Libs/Python/ShapeworksPython.cpp b/Libs/Python/ShapeworksPython.cpp index ec1bc9312c..ba606bcbf1 100644 --- a/Libs/Python/ShapeworksPython.cpp +++ b/Libs/Python/ShapeworksPython.cpp @@ -1119,7 +1119,7 @@ PYBIND11_MODULE(shapeworks_py, m) { .def_static( "prepareMesh", - [](const Mesh& mesh) -> decltype(auto) { return Mesh(MeshWarper::prep_mesh(mesh.getVTKMesh())); }, + [](const Mesh& mesh) -> decltype(auto) { return Mesh(MeshUtils::repair_mesh(mesh.getVTKMesh())); }, "Return the prepared mesh used for warping (before vertices were inserted).", "mesh"_a) .def( From 93148ada1b1c67bcaefba9f8d720d7cbf7593ac0 Mon Sep 17 00:00:00 2001 From: Alan Morris Date: Thu, 18 Dec 2025 14:16:49 -0700 Subject: [PATCH 3/6] Remove blocking integrity check from computeNormals() --- Libs/Mesh/Mesh.cpp | 17 ++--------------- 1 file changed, 2 insertions(+), 15 deletions(-) diff --git a/Libs/Mesh/Mesh.cpp b/Libs/Mesh/Mesh.cpp index d105f3eac3..5f56cd8ce5 100644 --- a/Libs/Mesh/Mesh.cpp +++ b/Libs/Mesh/Mesh.cpp @@ -731,11 +731,8 @@ Mesh& Mesh::clipClosedSurface(const Plane plane) { } Mesh& Mesh::computeNormals() { - - std::string issues = checkIntegrity(); - if (!issues.empty()) { - std::cerr << "Mesh integrity check failed:\n" << issues << std::endl; - // Either return early or try to fix + // Validate mesh first + if (!poly_data_ || poly_data_->GetNumberOfPoints() == 0 || poly_data_->GetNumberOfCells() == 0) { return *this; } @@ -743,16 +740,6 @@ Mesh& Mesh::computeNormals() { poly_data_->GetPointData()->SetNormals(nullptr); poly_data_->GetCellData()->SetNormals(nullptr); - // Validate mesh first - if (!poly_data_ || poly_data_->GetNumberOfPoints() == 0) { - std::cerr << "Warning: Empty mesh in computeNormals()" << std::endl; - return *this; - } - if (poly_data_->GetNumberOfCells() == 0) { - std::cerr << "Warning: Mesh has no cells in computeNormals()" << std::endl; - return *this; - } - auto normal = vtkSmartPointer::New(); normal->SetInputData(this->poly_data_); From 35f5f7c61e3c07219f687ea54ec6565ac00aca48 Mon Sep 17 00:00:00 2001 From: Alan Morris Date: Thu, 18 Dec 2025 19:55:58 -0700 Subject: [PATCH 4/6] Replace extractLargestComponent grooming step with repair_mesh, remove UI option, silently ignore deprecated parameter --- Libs/Groom/Groom.cpp | 15 ++++-------- Libs/Groom/GroomParameters.cpp | 21 +++++++--------- Libs/Groom/GroomParameters.h | 3 --- Libs/Mesh/Mesh.cpp | 10 ++------ Studio/Groom/GroomTool.cpp | 2 -- Studio/Groom/GroomTool.ui | 44 ---------------------------------- 6 files changed, 15 insertions(+), 80 deletions(-) diff --git a/Libs/Groom/Groom.cpp b/Libs/Groom/Groom.cpp index fea1d5c6e1..5246e48e97 100644 --- a/Libs/Groom/Groom.cpp +++ b/Libs/Groom/Groom.cpp @@ -368,14 +368,8 @@ bool Groom::mesh_pipeline(std::shared_ptr subject, size_t domain) { //--------------------------------------------------------------------------- bool Groom::run_mesh_pipeline(Mesh& mesh, GroomParameters params, const std::string& filename) { - check_and_fix_mesh(mesh, "initial", filename); - - // Extract largest component (should be first) - if (params.get_mesh_largest_component()) { - mesh.extractLargestComponent(); - check_and_fix_mesh(mesh, "largest_component", filename); - increment_progress(); - } + // 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(); @@ -508,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; @@ -1065,13 +1058,13 @@ void Groom::assign_transforms(std::vector> 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 at step '{}', file '{}', issues: {}", step, filename, issues); + 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 at step '{}', file '{}', issues: {}", step, filename, issues); + 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)); } } diff --git a/Libs/Groom/GroomParameters.cpp b/Libs/Groom/GroomParameters.cpp index 79035407b7..aaa2c57347 100644 --- a/Libs/Groom/GroomParameters.cpp +++ b/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 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 deprecated_params = { + "mesh_largest_component", // Now handled by repair_mesh at start of pipeline + }; + std::vector 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); } diff --git a/Libs/Groom/GroomParameters.h b/Libs/Groom/GroomParameters.h index 94c63fe955..346d0f88f9 100644 --- a/Libs/Groom/GroomParameters.h +++ b/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); diff --git a/Libs/Mesh/Mesh.cpp b/Libs/Mesh/Mesh.cpp index 5f56cd8ce5..9ef169168d 100644 --- a/Libs/Mesh/Mesh.cpp +++ b/Libs/Mesh/Mesh.cpp @@ -418,14 +418,9 @@ Mesh& Mesh::fillHoles(double hole_size) { filter->Update(); this->poly_data_ = filter->GetOutput(); - auto origNormal = poly_data_->GetPointData()->GetNormals(); - - // 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; } @@ -1766,8 +1761,7 @@ void Mesh::interpolate_scalars_to_mesh(std::string name, Eigen::VectorXd positio setField(name, scalars, Mesh::FieldType::Point); } -std::string Mesh::checkIntegrity() const -{ +std::string Mesh::checkIntegrity() const { std::ostringstream issues; if (!poly_data_) { diff --git a/Studio/Groom/GroomTool.cpp b/Studio/Groom/GroomTool.cpp index 4b30d82c4e..68ea36c6ea 100644 --- a/Studio/Groom/GroomTool.cpp +++ b/Studio/Groom/GroomTool.cpp @@ -489,7 +489,6 @@ void GroomTool::set_ui_from_params(GroomParameters params) { ui_->blur_checkbox->setChecked(params.get_blur_tool()); ui_->isolate_checkbox->setChecked(params.get_isolate_tool()); ui_->fill_holes_checkbox->setChecked(params.get_fill_holes_tool()); - ui_->mesh_largest_component->setChecked(params.get_mesh_largest_component()); ui_->mesh_fill_holes->setChecked(params.get_fill_mesh_holes_tool()); ui_->antialias_iterations->setValue(params.get_antialias_iterations()); ui_->blur_sigma->setValue(params.get_blur_amount()); @@ -587,7 +586,6 @@ void GroomTool::store_params() { params.set_blur_amount(ui_->blur_sigma->value()); params.set_isolate_tool(ui_->isolate_checkbox->isChecked()); params.set_fill_holes_tool(ui_->fill_holes_checkbox->isChecked()); - params.set_mesh_largest_component(ui_->mesh_largest_component->isChecked()); params.set_fill_mesh_holes_tool(ui_->mesh_fill_holes->isChecked()); params.set_antialias_iterations(ui_->antialias_iterations->value()); params.set_groom_output_prefix(preferences_.get_groom_file_template().toStdString()); diff --git a/Studio/Groom/GroomTool.ui b/Studio/Groom/GroomTool.ui index d28c7a7767..d3e79ecb53 100644 --- a/Studio/Groom/GroomTool.ui +++ b/Studio/Groom/GroomTool.ui @@ -1397,50 +1397,6 @@ QWidget#domain_panel { 0 - - - - - 0 - 30 - - - - - 0 - - - 6 - - - 0 - - - 6 - - - 0 - - - - - Extract Largest Component - - - true - - - - - - - - - - Qt::Horizontal - - - From 363ab515f211d6492a40c8923abc86b0945db8c7 Mon Sep 17 00:00:00 2001 From: Alan Morris Date: Thu, 18 Dec 2025 20:02:05 -0700 Subject: [PATCH 5/6] Fix fillHoles to remove zero-area triangles created by vtkFillHolesFilter --- Libs/Mesh/Mesh.cpp | 3 +++ Testing/data/fillholes.vtk | 4 ++-- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/Libs/Mesh/Mesh.cpp b/Libs/Mesh/Mesh.cpp index 9ef169168d..5eeb3cc3ac 100644 --- a/Libs/Mesh/Mesh.cpp +++ b/Libs/Mesh/Mesh.cpp @@ -418,6 +418,9 @@ Mesh& Mesh::fillHoles(double hole_size) { filter->Update(); this->poly_data_ = filter->GetOutput(); + // 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 and recompute normals computeNormals(); diff --git a/Testing/data/fillholes.vtk b/Testing/data/fillholes.vtk index 224041f0b8..8c885aa4b9 100644 --- a/Testing/data/fillholes.vtk +++ b/Testing/data/fillholes.vtk @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:e912e0b0af3d2117546b7126ada975d6fbf9d12b57d35f9b1e48716e1e328396 -size 488838 +oid sha256:7d978c50ae5c577fb42ce32b54baf1561c36e72c8be4a75731cbcf2027d4d4f1 +size 1036658 From a36663c0e0d8ce6c1774c69b8ca90891730c6fe8 Mon Sep 17 00:00:00 2001 From: Alan Morris Date: Thu, 18 Dec 2025 20:31:40 -0700 Subject: [PATCH 6/6] Fix tab order on mesh grooming page for sinc smoothing options --- Studio/Groom/GroomTool.ui | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Studio/Groom/GroomTool.ui b/Studio/Groom/GroomTool.ui index d3e79ecb53..fe73fe92b2 100644 --- a/Studio/Groom/GroomTool.ui +++ b/Studio/Groom/GroomTool.ui @@ -2156,6 +2156,8 @@ QWidget#domain_panel { mesh_smooth_method laplacian_iterations laplacian_relaxation + sinc_iterations + sinc_passband remesh_checkbox remesh_percent_checkbox remesh_percent_slider @@ -2169,8 +2171,6 @@ QWidget#domain_panel { reflect_axis alignment_image_checkbox alignment_box - sinc_passband - sinc_iterations