From a4ba636effac0bb6f3237480f29eadea3076e3a5 Mon Sep 17 00:00:00 2001 From: Alan Morris Date: Fri, 22 Aug 2025 00:57:30 -0600 Subject: [PATCH 01/18] Add UI elements for multiple shared boundaries. --- Studio/Groom/GroomTool.ui | 53 ++++++++++++++++++++++++++------------- 1 file changed, 35 insertions(+), 18 deletions(-) diff --git a/Studio/Groom/GroomTool.ui b/Studio/Groom/GroomTool.ui index bd3ea50878..07cf119084 100644 --- a/Studio/Groom/GroomTool.ui +++ b/Studio/Groom/GroomTool.ui @@ -7,7 +7,7 @@ 0 0 500 - 1262 + 1361 @@ -351,10 +351,37 @@ QWidget#domain_panel { 6 - + + + + + + + + + + Tolerance + + + + + + + First Domain + + + + + + + Second Domain + + + + @@ -370,30 +397,20 @@ QWidget#domain_panel { - - - - Tolerance - - - - - + + - First Domain + Delete Selected - - + + - Second Domain + Add - - - From c8298864986fa83e87ad5ba35bd63da4a2612aea Mon Sep 17 00:00:00 2001 From: Alan Morris Date: Sat, 23 Aug 2025 11:03:15 -0600 Subject: [PATCH 02/18] Adding support for multiple shared boundaries --- Libs/Groom/Groom.cpp | 201 +++++++++++++------------- Libs/Groom/GroomParameters.cpp | 209 +++++++++++++++++++-------- Libs/Groom/GroomParameters.h | 28 ++-- Libs/Mesh/MeshUtils.cpp | 150 +++++++++---------- Studio/Groom/GroomTool.cpp | 124 ++++++++++++++-- Studio/Groom/GroomTool.h | 4 + Studio/Groom/GroomTool.ui | 14 +- Studio/Utils/StudioUtils.cpp | 12 +- Studio/Utils/StudioUtils.h | 3 +- Studio/Visualization/PlaneWidget.cpp | 2 +- Studio/Visualization/Viewer.cpp | 3 + 11 files changed, 483 insertions(+), 267 deletions(-) diff --git a/Libs/Groom/Groom.cpp b/Libs/Groom/Groom.cpp index b076a98307..2347643e87 100644 --- a/Libs/Groom/Groom.cpp +++ b/Libs/Groom/Groom.cpp @@ -698,7 +698,7 @@ bool Groom::run_alignment() { bool Groom::run_shared_boundaries() { // only need to check on one domain auto params = GroomParameters(project_, project_->get_domain_names()[0]); - if (!params.get_shared_boundary()) { + if (!params.get_shared_boundaries_enabled()) { return true; } @@ -708,123 +708,124 @@ bool Groom::run_shared_boundaries() { auto original_domain_types = project_->get_original_domain_types(); auto groomed_domain_types = project_->get_groomed_domain_types(); - std::string first_domain_name = params.get_shared_boundary_first_domain(); - std::string second_domain_name = params.get_shared_boundary_second_domain(); - std::string shared_surface_name = "shared_surface"; - std::string shared_boundary_name = "shared_boundary"; - - // if domain_names doesn't have the shared domains, add them - if (std::find(domain_names.begin(), domain_names.end(), shared_surface_name) == domain_names.end()) { - domain_names.push_back(shared_surface_name); - original_domain_types.push_back(DomainType::Mesh); - groomed_domain_types.push_back(DomainType::Mesh); - } - if (std::find(domain_names.begin(), domain_names.end(), shared_boundary_name) == domain_names.end()) { - domain_names.push_back(shared_boundary_name); - original_domain_types.push_back(DomainType::Contour); - groomed_domain_types.push_back(DomainType::Contour); - } - - int shared_domain_index = - std::find(domain_names.begin(), domain_names.end(), shared_surface_name) - domain_names.begin(); - int shared_boundary_index = - std::find(domain_names.begin(), domain_names.end(), shared_boundary_name) - domain_names.begin(); - - // find the index of the first domain and second domain - int first_domain = 0; - int second_domain = 1; - for (int i = 0; i < domain_names.size(); i++) { - if (domain_names[i] == first_domain_name) { - first_domain = i; + for (auto shared_boundary : params.get_shared_boundaries()) { + std::string first_domain_name = shared_boundary.first_domain; + std::string second_domain_name = shared_boundary.second_domain; + std::string shared_surface_name = "shared_surface_" + first_domain_name + "_" + second_domain_name; + std::string shared_boundary_name = "shared_boundary_" + first_domain_name + "_" + second_domain_name; + + // if domain_names doesn't have the shared domains, add them + if (std::find(domain_names.begin(), domain_names.end(), shared_surface_name) == domain_names.end()) { + domain_names.push_back(shared_surface_name); + original_domain_types.push_back(DomainType::Mesh); + groomed_domain_types.push_back(DomainType::Mesh); } - if (domain_names[i] == second_domain_name) { - second_domain = i; + if (std::find(domain_names.begin(), domain_names.end(), shared_boundary_name) == domain_names.end()) { + domain_names.push_back(shared_boundary_name); + original_domain_types.push_back(DomainType::Contour); + groomed_domain_types.push_back(DomainType::Contour); } - } - - groomed_domain_types[first_domain] = DomainType::Mesh; - groomed_domain_types[second_domain] = DomainType::Mesh; - project_->set_domain_names(domain_names); - project_->set_original_domain_types(original_domain_types); - project_->set_groomed_domain_types(groomed_domain_types); + int shared_domain_index = + std::find(domain_names.begin(), domain_names.end(), shared_surface_name) - domain_names.begin(); + int shared_boundary_index = + std::find(domain_names.begin(), domain_names.end(), shared_boundary_name) - domain_names.begin(); + + // find the index of the first domain and second domain + int first_domain = 0; + int second_domain = 1; + for (int i = 0; i < domain_names.size(); i++) { + if (domain_names[i] == first_domain_name) { + first_domain = i; + } + if (domain_names[i] == second_domain_name) { + second_domain = i; + } + } - auto subjects = project_->get_subjects(); + groomed_domain_types[first_domain] = DomainType::Mesh; + groomed_domain_types[second_domain] = DomainType::Mesh; - std::mutex progress_mutex; - int progress = 0; + project_->set_domain_names(domain_names); + project_->set_original_domain_types(original_domain_types); + project_->set_groomed_domain_types(groomed_domain_types); - tbb::parallel_for(tbb::blocked_range{0, subjects.size()}, [&](const tbb::blocked_range& r) { - for (size_t i = r.begin(); i < r.end(); ++i) { - if (abort_) { - return; - } + auto subjects = project_->get_subjects(); - auto subject = subjects[i]; + std::mutex progress_mutex; + int progress = 0; - Mesh first_mesh = get_mesh(i, first_domain, false); - Mesh second_mesh = get_mesh(i, second_domain, false); + tbb::parallel_for(tbb::blocked_range{0, subjects.size()}, [&](const tbb::blocked_range& r) { + for (size_t i = r.begin(); i < r.end(); ++i) { + if (abort_) { + return; + } - auto [extracted_l, extracted_r, extracted_s] = - MeshUtils::sharedBoundaryExtractor(first_mesh, second_mesh, params.get_shared_boundary_tolerance()); + auto subject = subjects[i]; - extracted_l.remeshPercent(.99, 1.0); - extracted_r.remeshPercent(.99, 1.0); + Mesh first_mesh = get_mesh(i, first_domain, false); + Mesh second_mesh = get_mesh(i, second_domain, false); - auto output_contour = MeshUtils::boundaryLoopExtractor(extracted_s); + auto [extracted_l, extracted_r, extracted_s] = + MeshUtils::sharedBoundaryExtractor(first_mesh, second_mesh, shared_boundary.tolerance); - // overwrite groomed filename for first and second with extracted_l and extracted_r - auto first_filename = subject->get_groomed_filenames()[first_domain]; - auto second_filename = subject->get_groomed_filenames()[second_domain]; - // change file extension to .vtk if it's not already - first_filename = StringUtils::removeExtension(first_filename) + "_extracted.vtk"; - second_filename = StringUtils::removeExtension(second_filename) + "_extracted.vtk"; - MeshUtils::threadSafeWriteMesh(first_filename, extracted_l); - MeshUtils::threadSafeWriteMesh(second_filename, extracted_r); + extracted_l.remeshPercent(.99, 1.0); + extracted_r.remeshPercent(.99, 1.0); - auto shared_surface_filename = get_output_filename("shared_surface", DomainType::Mesh); - auto shared_boundary_filename = get_output_filename("shared_boundary", DomainType::Contour); + auto output_contour = MeshUtils::boundaryLoopExtractor(extracted_s); - extracted_s.write(shared_surface_filename); - output_contour.write(shared_boundary_filename); + // overwrite groomed filename for first and second with extracted_l and extracted_r + auto first_filename = subject->get_groomed_filenames()[first_domain]; + auto second_filename = subject->get_groomed_filenames()[second_domain]; + // change file extension to .vtk if it's not already + first_filename = StringUtils::removeExtension(first_filename) + "_extracted.vtk"; + second_filename = StringUtils::removeExtension(second_filename) + "_extracted.vtk"; + MeshUtils::threadSafeWriteMesh(first_filename, extracted_l); + MeshUtils::threadSafeWriteMesh(second_filename, extracted_r); - // store filenames - auto groomed_filenames = subject->get_groomed_filenames(); - groomed_filenames[first_domain] = first_filename; - groomed_filenames[second_domain] = second_filename; - if (shared_domain_index >= groomed_filenames.size()) { - groomed_filenames.resize(shared_domain_index + 1); - } - if (shared_boundary_index >= groomed_filenames.size()) { - groomed_filenames.resize(shared_boundary_index + 1); - } - groomed_filenames[shared_domain_index] = shared_surface_filename; - groomed_filenames[shared_boundary_index] = shared_boundary_filename; - subject->set_groomed_filenames(groomed_filenames); - subject->set_number_of_domains(domain_names.size()); + auto shared_surface_filename = get_output_filename("shared_surface", DomainType::Mesh); + auto shared_boundary_filename = get_output_filename("shared_boundary", DomainType::Contour); - // also store the two new ones as original - auto original_filenames = subject->get_original_filenames(); - if (original_filenames.size() <= shared_domain_index) { - original_filenames.resize(shared_domain_index + 1); - } - if (original_filenames.size() <= shared_boundary_index) { - original_filenames.resize(shared_boundary_index + 1); - } - original_filenames[shared_domain_index] = shared_surface_filename; - original_filenames[shared_boundary_index] = shared_boundary_filename; - subject->set_original_filenames(original_filenames); + extracted_s.write(shared_surface_filename); + output_contour.write(shared_boundary_filename); - { - // lock - std::scoped_lock lock(progress_mutex); - progress++; - progress_ = static_cast(progress) / static_cast(subjects.size()) * 100.0; - SW_PROGRESS(progress_, fmt::format("Processing shared boundaries ({}/{})", progress, subjects.size())); + // store filenames + auto groomed_filenames = subject->get_groomed_filenames(); + groomed_filenames[first_domain] = first_filename; + groomed_filenames[second_domain] = second_filename; + if (shared_domain_index >= groomed_filenames.size()) { + groomed_filenames.resize(shared_domain_index + 1); + } + if (shared_boundary_index >= groomed_filenames.size()) { + groomed_filenames.resize(shared_boundary_index + 1); + } + groomed_filenames[shared_domain_index] = shared_surface_filename; + groomed_filenames[shared_boundary_index] = shared_boundary_filename; + subject->set_groomed_filenames(groomed_filenames); + subject->set_number_of_domains(domain_names.size()); + + // also store the two new ones as original + auto original_filenames = subject->get_original_filenames(); + if (original_filenames.size() <= shared_domain_index) { + original_filenames.resize(shared_domain_index + 1); + } + if (original_filenames.size() <= shared_boundary_index) { + original_filenames.resize(shared_boundary_index + 1); + } + original_filenames[shared_domain_index] = shared_surface_filename; + original_filenames[shared_boundary_index] = shared_boundary_filename; + subject->set_original_filenames(original_filenames); + + { + // lock + std::scoped_lock lock(progress_mutex); + progress++; + progress_ = static_cast(progress) / static_cast(subjects.size()) * 100.0; + SW_PROGRESS(progress_, fmt::format("Processing shared boundaries ({}/{})", progress, subjects.size())); + } } - } - }); - + }); + } return true; } diff --git a/Libs/Groom/GroomParameters.cpp b/Libs/Groom/GroomParameters.cpp index 0dfa3f645b..b4d0923cb3 100644 --- a/Libs/Groom/GroomParameters.cpp +++ b/Libs/Groom/GroomParameters.cpp @@ -62,6 +62,7 @@ const std::string SHARED_BOUNDARY = "shared_boundary"; const std::string SHARED_BOUNDARY_FIRST_DOMAIN = "shared_boundary_first_domain"; const std::string SHARED_BOUNDARY_SECOND_DOMAIN = "shared_boundary_second_domain"; const std::string SHARED_BOUNDARY_TOLERANCE = "shared_boundary_tolerance"; +const std::string SHARED_BOUNDARIES = "shared_boundaries"; } // namespace Keys @@ -107,59 +108,90 @@ const bool shared_boundary = false; const std::string shared_boundary_first_domain = ""; const std::string shared_boundary_second_domain = ""; const double shared_boundary_tolerance = 1e-1; +const std::string shared_boundaries = ""; // Empty string = no boundaries } // namespace Defaults +// Helper methods for SharedBoundary struct +//--------------------------------------------------------------------------- +std::string GroomParameters::SharedBoundary::to_string() const { + return first_domain + "|" + second_domain + "|" + std::to_string(tolerance); +} + +//--------------------------------------------------------------------------- +GroomParameters::SharedBoundary GroomParameters::SharedBoundary::from_string(const std::string& str) { + std::vector parts; + std::stringstream ss(str); + std::string item; + + while (std::getline(ss, item, '|')) { + parts.push_back(item); + } + + if (parts.size() != 3) { + throw std::runtime_error("Invalid shared boundary string format"); + } + + SharedBoundary boundary; + boundary.first_domain = parts[0]; + boundary.second_domain = parts[1]; + boundary.tolerance = std::stod(parts[2]); + return boundary; +} + //--------------------------------------------------------------------------- GroomParameters::GroomParameters(ProjectHandle project, std::string domain_name) : project_(project), domain_name_(domain_name) { params_ = project_->get_parameters(Parameters::GROOM_PARAMS, domain_name_); - std::vector all_params = {Keys::CROP, - Keys::REFLECT, - Keys::REFLECT_COLUMN, - Keys::REFLECT_CHOICE, - Keys::REFLECT_AXIS, - Keys::RESAMPLE, - Keys::ISOTROPIC, - Keys::ISO_SPACING, - Keys::SPACING, - Keys::CONVERT_MESH, - Keys::FILL_MESH_HOLES, - Keys::FILL_HOLES, - Keys::ISOLATE, - Keys::PAD, - Keys::PAD_VALUE, - Keys::ANTIALIAS, - Keys::ANTIALIAS_AMOUNT, - Keys::BLUR, - Keys::BLUR_SIGMA, - Keys::FASTMARCHING, - Keys::MESH_SMOOTH, - Keys::MESH_SMOOTHING_METHOD, - Keys::MESH_SMOOTHING_VTK_LAPLACIAN_ITERATIONS, - Keys::MESH_SMOOTHING_VTK_LAPLACIAN_RELAXATION, - Keys::MESH_SMOOTHING_VTK_WINDOWED_SINC_ITERATIONS, - Keys::MESH_SMOOTHING_VTK_WINDOWED_SINC_PASSBAND, - Keys::ALIGNMENT_METHOD, - Keys::ALIGNMENT_ENABLED, - Keys::ALIGNMENT_REFERENCE, - Keys::ALIGNMENT_REFERENCE_CHOSEN, - Keys::ALIGNMENT_SUBSET_SIZE, - Keys::GROOM_OUTPUT_PREFIX, - Keys::REMESH, - Keys::REMESH_PERCENT_MODE, - Keys::REMESH_PERCENT, - Keys::REMESH_NUM_VERTICES, - Keys::REMESH_GRADATION, - Keys::GROOM_ALL_DOMAINS_THE_SAME, - Keys::SKIP_GROOMING, - Keys::CENTER, - Keys::ICP, - Keys::SHARED_BOUNDARY, - Keys::SHARED_BOUNDARY_FIRST_DOMAIN, - Keys::SHARED_BOUNDARY_SECOND_DOMAIN, - Keys::SHARED_BOUNDARY_TOLERANCE}; + std::vector all_params = { + Keys::CROP, + Keys::REFLECT, + Keys::REFLECT_COLUMN, + Keys::REFLECT_CHOICE, + Keys::REFLECT_AXIS, + Keys::RESAMPLE, + Keys::ISOTROPIC, + Keys::ISO_SPACING, + Keys::SPACING, + Keys::CONVERT_MESH, + Keys::FILL_MESH_HOLES, + Keys::FILL_HOLES, + Keys::ISOLATE, + Keys::PAD, + Keys::PAD_VALUE, + Keys::ANTIALIAS, + Keys::ANTIALIAS_AMOUNT, + Keys::BLUR, + Keys::BLUR_SIGMA, + Keys::FASTMARCHING, + Keys::MESH_SMOOTH, + Keys::MESH_SMOOTHING_METHOD, + Keys::MESH_SMOOTHING_VTK_LAPLACIAN_ITERATIONS, + Keys::MESH_SMOOTHING_VTK_LAPLACIAN_RELAXATION, + Keys::MESH_SMOOTHING_VTK_WINDOWED_SINC_ITERATIONS, + Keys::MESH_SMOOTHING_VTK_WINDOWED_SINC_PASSBAND, + Keys::ALIGNMENT_METHOD, + Keys::ALIGNMENT_ENABLED, + Keys::ALIGNMENT_REFERENCE, + Keys::ALIGNMENT_REFERENCE_CHOSEN, + Keys::ALIGNMENT_SUBSET_SIZE, + Keys::GROOM_OUTPUT_PREFIX, + Keys::REMESH, + Keys::REMESH_PERCENT_MODE, + Keys::REMESH_PERCENT, + Keys::REMESH_NUM_VERTICES, + Keys::REMESH_GRADATION, + Keys::GROOM_ALL_DOMAINS_THE_SAME, + Keys::SKIP_GROOMING, + Keys::CENTER, + Keys::ICP, + Keys::SHARED_BOUNDARY, + Keys::SHARED_BOUNDARY_FIRST_DOMAIN, + Keys::SHARED_BOUNDARY_SECOND_DOMAIN, + Keys::SHARED_BOUNDARY_TOLERANCE, + Keys::SHARED_BOUNDARIES, + }; std::vector to_remove; @@ -497,40 +529,93 @@ bool GroomParameters::get_skip_grooming() { return params_.get(Keys::SKIP_GROOMI void GroomParameters::set_skip_grooming(bool skip) { params_.set(Keys::SKIP_GROOMING, skip); } //--------------------------------------------------------------------------- -bool GroomParameters::get_shared_boundary() { return params_.get(Keys::SHARED_BOUNDARY, Defaults::shared_boundary); } +bool GroomParameters::get_shared_boundaries_enabled() { + return params_.get(Keys::SHARED_BOUNDARY, Defaults::shared_boundary); +} //--------------------------------------------------------------------------- -void GroomParameters::set_shared_boundary(bool shared_boundary) { params_.set(Keys::SHARED_BOUNDARY, shared_boundary); } +void GroomParameters::set_shared_boundaries_enabled(bool enabled) { params_.set(Keys::SHARED_BOUNDARY, enabled); } //--------------------------------------------------------------------------- -std::string GroomParameters::get_shared_boundary_first_domain() { - return params_.get(Keys::SHARED_BOUNDARY_FIRST_DOMAIN, Defaults::shared_boundary_first_domain); -} +std::vector GroomParameters::get_shared_boundaries() { + std::vector boundaries; -//--------------------------------------------------------------------------- -void GroomParameters::set_shared_boundary_first_domain(const std::string& domain_name) { - params_.set(Keys::SHARED_BOUNDARY_FIRST_DOMAIN, domain_name); + // Try new format first + if (params_.key_exists(Keys::SHARED_BOUNDARIES)) { + std::string boundaries_str = params_.get(Keys::SHARED_BOUNDARIES, Defaults::shared_boundaries); + if (!boundaries_str.empty()) { + std::stringstream ss(boundaries_str); + std::string boundary_str; + + while (std::getline(ss, boundary_str, ';')) { + if (!boundary_str.empty()) { + try { + boundaries.push_back(SharedBoundary::from_string(boundary_str)); + } catch (const std::exception& e) { + SW_WARN("Failed to parse shared boundary: " + boundary_str); + } + } + } + } + return boundaries; + } + + // Migration: convert old single boundary format + if (params_.get(Keys::SHARED_BOUNDARY, Defaults::shared_boundary)) { + SharedBoundary boundary; + boundary.first_domain = + std::string(params_.get(Keys::SHARED_BOUNDARY_FIRST_DOMAIN, Defaults::shared_boundary_first_domain)); + boundary.second_domain = + std::string(params_.get(Keys::SHARED_BOUNDARY_SECOND_DOMAIN, Defaults::shared_boundary_second_domain)); + boundary.tolerance = params_.get(Keys::SHARED_BOUNDARY_TOLERANCE, Defaults::shared_boundary_tolerance); + boundaries.push_back(boundary); + + // Migrate to new format + set_shared_boundaries(boundaries); + } + + return boundaries; } //--------------------------------------------------------------------------- -std::string GroomParameters::get_shared_boundary_second_domain() { - return params_.get(Keys::SHARED_BOUNDARY_SECOND_DOMAIN, Defaults::shared_boundary_second_domain); +void GroomParameters::set_shared_boundaries(const std::vector& boundaries) { + if (boundaries.empty()) { + params_.set(Keys::SHARED_BOUNDARIES, ""); + } else { + std::string boundaries_str; + for (size_t i = 0; i < boundaries.size(); ++i) { + if (i > 0) boundaries_str += ";"; + boundaries_str += boundaries[i].to_string(); + } + params_.set(Keys::SHARED_BOUNDARIES, boundaries_str); + } + + // Clear old format keys to avoid confusion + params_.remove_entry(Keys::SHARED_BOUNDARY); + params_.remove_entry(Keys::SHARED_BOUNDARY_FIRST_DOMAIN); + params_.remove_entry(Keys::SHARED_BOUNDARY_SECOND_DOMAIN); + params_.remove_entry(Keys::SHARED_BOUNDARY_TOLERANCE); } //--------------------------------------------------------------------------- -void GroomParameters::set_shared_boundary_second_domain(const std::string& domain_name) { - params_.set(Keys::SHARED_BOUNDARY_SECOND_DOMAIN, domain_name); +void GroomParameters::add_shared_boundary(const std::string& first_domain, const std::string& second_domain, + double tolerance) { + auto boundaries = get_shared_boundaries(); + boundaries.push_back({first_domain, second_domain, tolerance}); + set_shared_boundaries(boundaries); } //--------------------------------------------------------------------------- -double GroomParameters::get_shared_boundary_tolerance() { - return params_.get(Keys::SHARED_BOUNDARY_TOLERANCE, Defaults::shared_boundary_tolerance); +void GroomParameters::remove_shared_boundary(size_t index) { + auto boundaries = get_shared_boundaries(); + if (index < boundaries.size()) { + boundaries.erase(boundaries.begin() + index); + set_shared_boundaries(boundaries); + } } //--------------------------------------------------------------------------- -void GroomParameters::set_shared_boundary_tolerance(double tolerance) { - params_.set(Keys::SHARED_BOUNDARY_TOLERANCE, tolerance); -} +void GroomParameters::clear_shared_boundaries() { set_shared_boundaries({}); } //--------------------------------------------------------------------------- } // namespace shapeworks diff --git a/Libs/Groom/GroomParameters.h b/Libs/Groom/GroomParameters.h index 6161840de3..346d0f88f9 100644 --- a/Libs/Groom/GroomParameters.h +++ b/Libs/Groom/GroomParameters.h @@ -11,6 +11,17 @@ namespace shapeworks { * This class encapsulated processing of Groom parameters */ class GroomParameters { + public: + struct SharedBoundary { + std::string first_domain; + std::string second_domain; + double tolerance; + + // Helper methods for serialization + std::string to_string() const; + static SharedBoundary from_string(const std::string& str); + }; + enum class MeshSmoothingOption { laplacian, sinc }; enum class AlignmentOption { none, center, icp }; @@ -137,17 +148,14 @@ class GroomParameters { bool get_skip_grooming(); void set_skip_grooming(bool skip); - bool get_shared_boundary(); - void set_shared_boundary(bool shared_boundary); - - std::string get_shared_boundary_first_domain(); - void set_shared_boundary_first_domain(const std::string& domain_name); - - std::string get_shared_boundary_second_domain(); - void set_shared_boundary_second_domain(const std::string& domain_name); + bool get_shared_boundaries_enabled(); + void set_shared_boundaries_enabled(bool enabled); - double get_shared_boundary_tolerance(); - void set_shared_boundary_tolerance(double tolerance); + std::vector get_shared_boundaries(); + void set_shared_boundaries(const std::vector& boundaries); + void add_shared_boundary(const std::string& first_domain, const std::string& second_domain, double tolerance); + void remove_shared_boundary(size_t index); + void clear_shared_boundaries(); void restore_defaults(); diff --git a/Libs/Mesh/MeshUtils.cpp b/Libs/Mesh/MeshUtils.cpp index 718e51efa1..ed53a0acde 100644 --- a/Libs/Mesh/MeshUtils.cpp +++ b/Libs/Mesh/MeshUtils.cpp @@ -1,5 +1,6 @@ #include "MeshUtils.h" +#include #include #include #include @@ -62,10 +63,8 @@ static std::vector get_random_subset(int num, int max) { return random_subset; } -const vtkSmartPointer MeshUtils::createICPTransform(const Mesh source, - const Mesh target, - Mesh::AlignmentType align, - const unsigned iterations, +const vtkSmartPointer MeshUtils::createICPTransform(const Mesh source, const Mesh target, + Mesh::AlignmentType align, const unsigned iterations, bool meshTransform) { if (source.numPoints() == 0 || target.numPoints() == 0) { throw std::invalid_argument("empty mesh passed to MeshUtils::createICPTransform"); @@ -97,9 +96,8 @@ const vtkSmartPointer MeshUtils::createICPTransform(const Mesh sou if (meshTransform) { m = icp->GetMatrix(); } else { - vtkMatrix4x4::Invert( - icp->GetMatrix(), - m); // It's inversed because when an image is transformed, + vtkMatrix4x4::Invert(icp->GetMatrix(), + m); // It's inversed because when an image is transformed, // a new image is created in the target space and samples through the transform back to the original space } @@ -108,14 +106,14 @@ const vtkSmartPointer MeshUtils::createICPTransform(const Mesh sou Mesh MeshUtils::create_mesh_from_file(std::string filename, double iso_value) { bool is_mesh = false; - for (auto& type: shapeworks::Mesh::getSupportedTypes()) { + for (auto& type : shapeworks::Mesh::getSupportedTypes()) { if (StringUtils::hasSuffix(filename, type)) { is_mesh = true; } } bool is_image = false; - for (auto& type: shapeworks::Image::getSupportedTypes()) { + for (auto& type : shapeworks::Image::getSupportedTypes()) { if (StringUtils::hasSuffix(filename, type)) { is_image = true; } @@ -147,7 +145,7 @@ PhysicalRegion MeshUtils::boundingBox(const std::vector& filenames, Mesh mesh(filenames[0]); PhysicalRegion bbox(mesh.boundingBox()); - for (auto& filename: filenames) { + for (auto& filename : filenames) { Mesh mesh(filename); bbox.expand(mesh.boundingBox()); } @@ -160,7 +158,7 @@ PhysicalRegion MeshUtils::boundingBox(const std::vector& meshes, bool cent PhysicalRegion bbox(meshes[0].boundingBox()); - for (auto& mesh: meshes) { + for (auto& mesh : meshes) { bbox.expand(mesh.boundingBox()); } @@ -206,32 +204,32 @@ int MeshUtils::findReferenceMesh(std::vector& meshes, int random_subset_si // mutex for access to results std::mutex mutex; - tbb::parallel_for(tbb::blocked_range{0, pairs.size()}, - [&](const tbb::blocked_range& r) { - for (size_t i = r.begin(); i < r.end(); ++i) { - auto pair = pairs[i]; - - vtkSmartPointer poly_data1 = vtkSmartPointer::New(); - poly_data1->DeepCopy(meshes[pair.first].getVTKMesh()); - vtkSmartPointer poly_data2 = vtkSmartPointer::New(); - poly_data2->DeepCopy(meshes[pair.second].getVTKMesh()); - - // register the two meshes - auto matrix = MeshUtils::createICPTransform(poly_data1, poly_data2, Mesh::Rigid, 10, true); - // transform - auto transform = createMeshTransform(matrix); - Mesh transformed = poly_data1; - transformed.applyTransform(transform); - - // compute distance - auto distances = transformed.distance(poly_data2)[0]; - double distance = mean(distances); { - // lock and store results - std::scoped_lock lock(mutex); - results[i] = distance; - } - } - }); + tbb::parallel_for(tbb::blocked_range{0, pairs.size()}, [&](const tbb::blocked_range& r) { + for (size_t i = r.begin(); i < r.end(); ++i) { + auto pair = pairs[i]; + + vtkSmartPointer poly_data1 = vtkSmartPointer::New(); + poly_data1->DeepCopy(meshes[pair.first].getVTKMesh()); + vtkSmartPointer poly_data2 = vtkSmartPointer::New(); + poly_data2->DeepCopy(meshes[pair.second].getVTKMesh()); + + // register the two meshes + auto matrix = MeshUtils::createICPTransform(poly_data1, poly_data2, Mesh::Rigid, 10, true); + // transform + auto transform = createMeshTransform(matrix); + Mesh transformed = poly_data1; + transformed.applyTransform(transform); + + // compute distance + auto distances = transformed.distance(poly_data2)[0]; + double distance = mean(distances); + { + // lock and store results + std::scoped_lock lock(mutex); + results[i] = distance; + } + } + }); std::vector means(meshes.size(), 0); @@ -255,9 +253,10 @@ int MeshUtils::findReferenceMesh(std::vector& meshes, int random_subset_si * Given a .ply mesh, extract the boundary loop and export the boundary loop as a VTK .vtp file */ +//--------------------------------------------------------------------------- static bool is_clockwise(const Eigen::MatrixXd& V, const Eigen::MatrixXi& F, const std::vector& loop) { Eigen::RowVector3d centroid{0.0, 0.0, 0.0}; - for (const auto& i: loop) { + for (const auto& i : loop) { centroid += V.row(i); } centroid /= loop.size(); @@ -271,6 +270,7 @@ static bool is_clockwise(const Eigen::MatrixXd& V, const Eigen::MatrixXi& F, con return angle0 > angle1; } +//--------------------------------------------------------------------------- Mesh MeshUtils::boundaryLoopExtractor(Mesh mesh) { Eigen::MatrixXd V = mesh.points(); Eigen::MatrixXi F = mesh.faces(); @@ -283,7 +283,7 @@ Mesh MeshUtils::boundaryLoopExtractor(Mesh mesh) { const auto is_cw = is_clockwise(V, F, loop); auto pts = vtkSmartPointer::New(); - for (const auto& i: loop) { + for (const auto& i : loop) { pts->InsertNextPoint(V(i, 0), V(i, 1), V(i, 2)); } @@ -318,6 +318,7 @@ Mesh MeshUtils::boundaryLoopExtractor(Mesh mesh) { * defines the threshold for two surfaces to be "close" */ +//--------------------------------------------------------------------------- static std::tuple rem_into_eigen_mesh(const std::vector& faces, const Eigen::MatrixXd& src_V, const Eigen::MatrixXi& src_F) { @@ -339,6 +340,7 @@ static std::tuple rem_into_eigen_mesh(const st return std::make_tuple(V, F); } +//--------------------------------------------------------------------------- static std::tuple shared_into_eigen_mesh(const std::vector& faces, const Eigen::MatrixXd& src_V, const Eigen::MatrixXi& src_F) { @@ -359,14 +361,13 @@ static std::tuple shared_into_eigen_mesh(const return std::make_tuple(V, F); } +//--------------------------------------------------------------------------- static bool is_empty(const Eigen::MatrixXd& V, const Eigen::MatrixXi& F) { return V.size() == 0 || F.size() == 0; } +//--------------------------------------------------------------------------- static std::tuple find_shared_surface( - const Eigen::MatrixXd& src_V, - const Eigen::MatrixXi& src_F, - const Eigen::MatrixXd& other_V, - const Eigen::MatrixXi& other_F, - double tol = 1e-3) { + const Eigen::MatrixXd& src_V, const Eigen::MatrixXi& src_F, const Eigen::MatrixXd& other_V, + const Eigen::MatrixXi& other_F, double tol = 1e-3) { Eigen::MatrixXd out_V; Eigen::MatrixXi out_F; Eigen::MatrixXd rem_V; @@ -406,18 +407,18 @@ static std::tuple > src_loops, shared_loops; igl::boundary_loop(src_F, src_loops); igl::boundary_loop(shared_F, shared_loops); - assert(src_loops.size() == 1); - assert(shared_loops.size() == 1); + if (src_loops.size() != 1 || shared_loops.size() != 1) { + SW_DEBUG("src_loops size: {}, shared_loops size: {}", src_loops.size(), shared_loops.size()); + throw std::runtime_error("Expected exactly one boundary loop in each mesh"); + } const auto& src_loop = src_loops[0]; const auto& shared_loop = shared_loops[0]; @@ -448,6 +449,7 @@ static void move_to_boundary(const Eigen::MatrixXd& src_V, } } +//--------------------------------------------------------------------------- std::array MeshUtils::sharedBoundaryExtractor(const Mesh& mesh_l, const Mesh& mesh_r, double tol) { Eigen::MatrixXd V_l, V_r; Eigen::MatrixXi F_l, F_r; @@ -462,7 +464,7 @@ std::array MeshUtils::sharedBoundaryExtractor(const Mesh& mesh_l, const std::tie(shared_V_r, shared_F_r, rem_V_r, rem_F_r) = find_shared_surface(V_r, F_r, V_l, F_l, tol); if (is_empty(shared_V_l, shared_F_l) || is_empty(shared_V_r, shared_F_r) || is_empty(rem_V_l, rem_F_l) || - is_empty(rem_V_r, rem_F_r)) { + is_empty(rem_V_r, rem_F_r)) { // todo this should return a status code to the caller so that it can be displayed or handled based on the // downstream task throw std::runtime_error("No shared surface detected. Please check the input meshes and/or increase the tolerance"); @@ -481,6 +483,7 @@ std::array MeshUtils::sharedBoundaryExtractor(const Mesh& mesh_l, const return std::array{std::move(out_l), std::move(out_r), std::move(out_s)}; } +//--------------------------------------------------------------------------- void MeshUtils::generateNormals(const std::vector >& meshes, bool forceRegen) { if (meshes.empty()) throw std::invalid_argument("No meshes provided to compute average normals"); @@ -496,29 +499,31 @@ void MeshUtils::generateNormals(const std::vector > } } +//--------------------------------------------------------------------------- Field MeshUtils::computeMeanNormals(const std::vector& filenames, bool autoGenerateNormals) { if (filenames.empty()) throw std::invalid_argument("No filenames provided to compute mean normals"); std::vector meshes; - meshes.reserve(filenames.size()); // create a vector large enough for all the meshes that will be loaded - for (auto& filename: filenames) { + meshes.reserve(filenames.size()); // create a vector large enough for all the meshes that will be loaded + for (auto& filename : filenames) { meshes.push_back(Mesh(filename)); } std::vector > rmeshes; rmeshes.reserve(meshes.size()); - for (Mesh& mesh: meshes) rmeshes.push_back(std::reference_wrapper(mesh)); + for (Mesh& mesh : meshes) rmeshes.push_back(std::reference_wrapper(mesh)); if (autoGenerateNormals) { MeshUtils::generateNormals(rmeshes); } std::vector > cmeshes; - for (Mesh& mesh: meshes) cmeshes.push_back(std::reference_wrapper(mesh)); + for (Mesh& mesh : meshes) cmeshes.push_back(std::reference_wrapper(mesh)); return computeMeanNormals(cmeshes); } +//--------------------------------------------------------------------------- Field MeshUtils::computeMeanNormals(const std::vector >& meshes) { if (meshes.empty()) throw std::invalid_argument("No meshes provided to compute average normals"); @@ -535,7 +540,7 @@ Field MeshUtils::computeMeanNormals(const std::vectorGetNumberOfTuples()) throw std::invalid_argument( - "Expected a normal for every point in mesh. Please call generateNormals to accomplish this"); + "Expected a normal for every point in mesh. Please call generateNormals to accomplish this"); for (int i = 0; i < num_normals; i++) { auto n = normals->GetTuple3(i); @@ -577,6 +582,7 @@ Field MeshUtils::computeMeanNormals(const std::vector mesh) { // std::cout << "VTK rendering" << std::endl; @@ -621,8 +627,7 @@ void MeshUtils::visualizeVectorFieldForFFCs(std::shared_ptr mesh) { // Computes grad vec for each face for (size_t i = 0; i < F.rows(); i++) { - const Eigen::Vector3d vert_dists(mesh->getFieldValue("value", F(i, 0)), - mesh->getFieldValue("value", F(i, 1)), + const Eigen::Vector3d vert_dists(mesh->getFieldValue("value", F(i, 0)), mesh->getFieldValue("value", F(i, 1)), mesh->getFieldValue("value", F(i, 2))); Eigen::Vector3d v1(V(F(i, 0), 0), V(F(i, 0), 1), V(F(i, 0), 2)); @@ -656,12 +661,14 @@ void MeshUtils::visualizeVectorFieldForFFCs(std::shared_ptr mesh) { renderWindowInteractor->Start(); } +//--------------------------------------------------------------------------- void CreateArrowMidPoint(double midPoint[], double startPoint[], double endPoint[]) { midPoint[0] = (startPoint[0] + endPoint[0]) / 3; midPoint[1] = (startPoint[1] + endPoint[1]) / 3; midPoint[2] = (startPoint[2] + endPoint[2]) / 3; } +//--------------------------------------------------------------------------- void PrintArray(std::string arrayName, double* arr, int size) { std::cout << arrayName << std::endl; @@ -672,22 +679,25 @@ void PrintArray(std::string arrayName, double* arr, int size) { std::cout << std::endl; } +//--------------------------------------------------------------------------- void PrintTransformParams(vtkSmartPointer transform) { PrintArray("Transform orientation: ", transform->GetOrientation(), 3); PrintArray("Transform position: ", transform->GetPosition(), 3); PrintArray("Transform scale: ", transform->GetScale(), 3); } +//--------------------------------------------------------------------------- void Print4x4Matrix(std::string matrixName, vtkSmartPointer matrix) { std::cout << matrixName << std::endl; for (unsigned int i = 0; i < sizeof(matrix->Element) / sizeof(matrix->Element[0]); i++) { std::cout << matrix->GetElement(i, 0) << " " << matrix->GetElement(i, 1) << " " << matrix->GetElement(i, 2) << " " - << matrix->GetElement(i, 3) << std::endl; + << matrix->GetElement(i, 3) << std::endl; } std::cout << std::endl; } +//--------------------------------------------------------------------------- vtkSmartPointer MeshUtils::getArrow(Eigen::Vector3d start, Eigen::Vector3d end) { // Create an arrow source vtkSmartPointer arrowSource = vtkSmartPointer::New(); @@ -778,14 +788,8 @@ vtkSmartPointer MeshUtils::getArrow(Eigen::Vector3d start, Eigen::Vect //------------------------------------------------------------------------- //! From vtkTriangle::EvaluatePosition -int MeshUtils::evaluate_triangle_position(const double x[3], - double closestPoint[3], - int& subId, - double pcoords[3], - double& dist2, - double weights[], - double pt3[3], - double pt1[3], +int MeshUtils::evaluate_triangle_position(const double x[3], double closestPoint[3], int& subId, double pcoords[3], + double& dist2, double weights[], double pt3[3], double pt1[3], double pt2[3]) { int i, j; double n[3], fabsn; @@ -794,7 +798,7 @@ int MeshUtils::evaluate_triangle_position(const double x[3], double maxComponent; int idx = 0, indices[2]; double dist2Point, dist2Line1, dist2Line2; - double* closest, closestPoint1[3], closestPoint2[3], cp[3]; + double *closest, closestPoint1[3], closestPoint2[3], cp[3]; subId = 0; pcoords[2] = 0.0; @@ -848,8 +852,8 @@ int MeshUtils::evaluate_triangle_position(const double x[3], weights[1] = pcoords[0]; weights[2] = pcoords[1]; - if (weights[0] >= 0.0 && weights[0] <= 1.0 && weights[1] >= 0.0 && weights[1] <= 1.0 && - weights[2] >= 0.0 && weights[2] <= 1.0) { + if (weights[0] >= 0.0 && weights[0] <= 1.0 && weights[1] >= 0.0 && weights[1] <= 1.0 && weights[2] >= 0.0 && + weights[2] <= 1.0) { // projection distance if (closestPoint) { dist2 = vtkMath::Distance2BetweenPoints(cp, x); @@ -926,4 +930,4 @@ int MeshUtils::evaluate_triangle_position(const double x[3], return 0; } } -} // namespace shapeworks +} // namespace shapeworks diff --git a/Studio/Groom/GroomTool.cpp b/Studio/Groom/GroomTool.cpp index aca6c9f2fd..71526f032a 100644 --- a/Studio/Groom/GroomTool.cpp +++ b/Studio/Groom/GroomTool.cpp @@ -38,6 +38,9 @@ GroomTool::GroomTool(Preferences& prefs, Telemetry& telemetry) : preferences_(pr connect(ui_->skip_grooming, &QCheckBox::toggled, this, &GroomTool::skip_grooming_toggled); + connect(ui_->add_shared_boundary, &QPushButton::clicked, this, &GroomTool::add_shared_boundary_clicked); + connect(ui_->delete_shared_boundary, &QPushButton::clicked, this, &GroomTool::delete_shared_boundary_clicked); + ui_->image_label->setAttribute(Qt::WA_TransparentForMouseEvents); ui_->mesh_label->setAttribute(Qt::WA_TransparentForMouseEvents); ui_->alignment_label->setAttribute(Qt::WA_TransparentForMouseEvents); @@ -174,6 +177,63 @@ void GroomTool::handle_error(QString msg) { enable_actions(); } +//--------------------------------------------------------------------------- +void GroomTool::add_shared_boundary_clicked() { + auto params = GroomParameters(session_->get_project(), ""); + + auto shared_boundaries = params.get_shared_boundaries(); + GroomParameters::SharedBoundary boundary; + boundary.first_domain = ui_->shared_boundary_first_domain->currentText().toStdString(); + boundary.second_domain = ui_->shared_boundary_second_domain->currentText().toStdString(); + boundary.tolerance = ui_->shared_boundary_tolerance->value(); + if (boundary.first_domain == boundary.second_domain) { + SW_ERROR("Cannot add a shared boundary between the same domain"); + return; + } + + for (const auto& existing_boundary : shared_boundaries) { + if (existing_boundary.first_domain == boundary.first_domain && + existing_boundary.second_domain == boundary.second_domain) { + SW_ERROR("Shared boundary already exists"); + return; + } + if (existing_boundary.first_domain == boundary.second_domain && + existing_boundary.second_domain == boundary.first_domain) { + SW_ERROR("Shared boundary already exists"); + return; + } + } + shared_boundaries.push_back(boundary); + params.set_shared_boundaries(shared_boundaries); + params.save_to_project(); + update_shared_boundary_table(); +} + +//--------------------------------------------------------------------------- +void GroomTool::delete_shared_boundary_clicked() { + auto params = GroomParameters(session_->get_project(), ""); + auto shared_boundaries = params.get_shared_boundaries(); + + QModelIndexList selected_rows = ui_->shared_boundary_table->selectionModel()->selectedRows(); + if (selected_rows.empty()) { + SW_ERROR("No shared boundary selected to delete"); + return; + } + + for (int i = selected_rows.size() - 1; i >= 0; i--) { + std::cerr << "erasing shared boundary: " << shared_boundaries[selected_rows[i].row()].first_domain << " <-> " + << shared_boundaries[selected_rows[i].row()].second_domain << " with tolerance " + << shared_boundaries[selected_rows[i].row()].tolerance << "\n"; + shared_boundaries.erase(shared_boundaries.begin() + selected_rows[i].row()); + } + + std::cerr << "There are now " << shared_boundaries.size() << " shared boundaries left.\n"; + + params.set_shared_boundaries(shared_boundaries); + params.save_to_project(); + update_shared_boundary_table(); +} + //--------------------------------------------------------------------------- void GroomTool::handle_progress(int val) { if (groom_is_running_) { @@ -233,8 +293,10 @@ void GroomTool::update_domain_box() { ui_->domain_box->clear(); } - StudioUtils::update_domain_combobox(ui_->shared_boundary_first_domain, session_); - StudioUtils::update_domain_combobox(ui_->shared_boundary_second_domain, session_); + StudioUtils::update_domain_combobox(ui_->shared_boundary_first_domain, session_, + {"shared_boundary", "shared_surface"}); + StudioUtils::update_domain_combobox(ui_->shared_boundary_second_domain, session_, + {"shared_boundary", "shared_surface"}); if (domain_names.size() > 1) { if (ui_->shared_boundary_first_domain->currentIndex() == ui_->shared_boundary_second_domain->currentIndex()) { ui_->shared_boundary_second_domain->setCurrentIndex(ui_->shared_boundary_second_domain->currentIndex() + 1); @@ -259,6 +321,44 @@ void GroomTool::apply_to_all_domains_changed() { update_domain_box(); } +//--------------------------------------------------------------------------- +void GroomTool::update_shared_boundary_table() { + auto params = GroomParameters(session_->get_project(), ""); + auto shared_boundaries = params.get_shared_boundaries(); + + std::cerr << "There are " << shared_boundaries.size() << " shared boundaries.\n"; + + QStringList table_headers; + table_headers << "First Domain"; + table_headers << "Second Domain"; + table_headers << "Tolerance"; + + QTableWidget* table = ui_->shared_boundary_table; + table->clear(); + table->setRowCount(shared_boundaries.size()); + table->setColumnCount(table_headers.size()); + + table->setHorizontalHeaderLabels(table_headers); + table->horizontalHeader()->setSectionResizeMode(QHeaderView::Stretch); + table->setSelectionBehavior(QAbstractItemView::SelectRows); + + int row = 0; + for (const auto& boundary : shared_boundaries) { + table->setItem(row, 0, new QTableWidgetItem(QString::fromStdString(boundary.first_domain))); + table->setItem(row, 1, new QTableWidgetItem(QString::fromStdString(boundary.second_domain))); + table->setItem(row, 2, new QTableWidgetItem(QString::number(boundary.tolerance))); + row++; + } + if (table->rowCount() < 1) { + table->setMaximumHeight(ui_->label->height() * 2); + } else { + table->setMaximumHeight(16777215); // default + } + table->resizeColumnsToContents(); + table->horizontalHeader()->setStretchLastSection(false); + table->setSelectionBehavior(QAbstractItemView::SelectRows); +} + //--------------------------------------------------------------------------- void GroomTool::update_reflect_columns() { if (!session_ || !session_->get_project()) { @@ -318,6 +418,7 @@ void GroomTool::load_params() { update_page(); update_reflect_columns(); + update_shared_boundary_table(); } //--------------------------------------------------------------------------- @@ -411,11 +512,10 @@ void GroomTool::set_ui_from_params(GroomParameters params) { ui_->remesh_gradation_slider->setValue(params.get_remesh_gradation() * 50.0); ui_->remesh_gradation_spinbox->setValue(params.get_remesh_gradation()); - ui_->shared_boundary->setChecked(params.get_shared_boundary()); - ui_->shared_boundary_first_domain->setCurrentText(QString::fromStdString(params.get_shared_boundary_first_domain())); - ui_->shared_boundary_second_domain->setCurrentText( - QString::fromStdString(params.get_shared_boundary_second_domain())); - ui_->shared_boundary_tolerance->setValue(params.get_shared_boundary_tolerance()); + ui_->shared_boundary->setChecked(params.get_shared_boundaries_enabled()); + ui_->shared_boundary_first_domain->setCurrentText(""); + ui_->shared_boundary_second_domain->setCurrentText(""); + ui_->shared_boundary_tolerance->setValue(1e-1); auto subjects = session_->get_project()->get_subjects(); int domain_id = std::max(ui_->domain_box->currentIndex(), 0); @@ -507,10 +607,7 @@ void GroomTool::store_params() { params.set_remesh_num_vertices(ui_->remesh_num_vertices->text().toInt()); params.set_remesh_gradation(ui_->remesh_gradation_spinbox->value()); - params.set_shared_boundary(ui_->shared_boundary->isChecked()); - params.set_shared_boundary_first_domain(ui_->shared_boundary_first_domain->currentText().toStdString()); - params.set_shared_boundary_second_domain(ui_->shared_boundary_second_domain->currentText().toStdString()); - params.set_shared_boundary_tolerance(ui_->shared_boundary_tolerance->value()); + params.set_shared_boundaries_enabled(ui_->shared_boundary->isChecked()); params.set_skip_grooming(ui_->skip_grooming->isChecked()); params.save_to_project(); @@ -531,10 +628,7 @@ void GroomTool::store_params() { params.set_groom_all_domains_the_same(ui_->apply_to_all_domains->isChecked()); params.set_skip_grooming(ui_->skip_grooming->isChecked()); - params.set_shared_boundary(ui_->shared_boundary->isChecked()); - params.set_shared_boundary_first_domain(ui_->shared_boundary_first_domain->currentText().toStdString()); - params.set_shared_boundary_second_domain(ui_->shared_boundary_second_domain->currentText().toStdString()); - params.set_shared_boundary_tolerance(ui_->shared_boundary_tolerance->value()); + params.set_shared_boundaries_enabled(ui_->shared_boundary->isChecked()); params.save_to_project(); } diff --git a/Studio/Groom/GroomTool.h b/Studio/Groom/GroomTool.h index 690eb220cb..d817cb4ded 100644 --- a/Studio/Groom/GroomTool.h +++ b/Studio/Groom/GroomTool.h @@ -78,12 +78,16 @@ class GroomTool : public QWidget { void handle_progress(int val); void handle_error(QString msg); + void add_shared_boundary_clicked(); + void delete_shared_boundary_clicked(); + private: void set_ui_from_params(GroomParameters params); void update_page(); void update_domain_box(); void apply_to_all_domains_changed(); + void update_shared_boundary_table(); void update_reflect_columns(); void update_reflect_choices(); diff --git a/Studio/Groom/GroomTool.ui b/Studio/Groom/GroomTool.ui index 07cf119084..69eba12305 100644 --- a/Studio/Groom/GroomTool.ui +++ b/Studio/Groom/GroomTool.ui @@ -7,7 +7,7 @@ 0 0 500 - 1361 + 1599 @@ -355,7 +355,7 @@ QWidget#domain_panel { - + @@ -389,6 +389,12 @@ QWidget#domain_panel { 24 + + Qt::AlignCenter + + + 3 + 99999.000000000000000 @@ -398,14 +404,14 @@ QWidget#domain_panel { - + Delete Selected - + Add diff --git a/Studio/Utils/StudioUtils.cpp b/Studio/Utils/StudioUtils.cpp index f92f6dfe5d..a3707d7709 100644 --- a/Studio/Utils/StudioUtils.cpp +++ b/Studio/Utils/StudioUtils.cpp @@ -178,13 +178,23 @@ void StudioUtils::window_width_level_to_brightness_contrast(double window_width, } //--------------------------------------------------------------------------- -void StudioUtils::update_domain_combobox(QComboBox* combobox, QSharedPointer session) { +void StudioUtils::update_domain_combobox(QComboBox* combobox, QSharedPointer session, + const std::vector& filters) { auto domain_names = session->get_project()->get_domain_names(); int currentIndex = combobox->currentIndex(); if (domain_names.size() != combobox->count()) { combobox->clear(); for (auto&& item : domain_names) { + bool skip = false; + for (const auto& filter : filters) { + if (QString::fromStdString(item).contains(filter, Qt::CaseInsensitive)) { + skip = true; + } + } + if (skip) { + continue; + } combobox->addItem(QString::fromStdString(item)); } } diff --git a/Studio/Utils/StudioUtils.h b/Studio/Utils/StudioUtils.h index ef86eb7326..9f1e12477d 100644 --- a/Studio/Utils/StudioUtils.h +++ b/Studio/Utils/StudioUtils.h @@ -50,7 +50,8 @@ class StudioUtils { double max_intensity, double& brightness, double& contrast); //! update a combobox with domain names - static void update_domain_combobox(QComboBox* combobox, QSharedPointer session); + static void update_domain_combobox(QComboBox* combobox, QSharedPointer session, + const std::vector& filters = {}); }; } // namespace shapeworks diff --git a/Studio/Visualization/PlaneWidget.cpp b/Studio/Visualization/PlaneWidget.cpp index 6fc6969a3d..56dde7472f 100644 --- a/Studio/Visualization/PlaneWidget.cpp +++ b/Studio/Visualization/PlaneWidget.cpp @@ -208,7 +208,7 @@ void PlaneWidget::update_planes() { double size = 100; auto meshes = viewer_->get_meshes(); - if (meshes.valid()) { + if (meshes.valid() && domain_id < meshes.meshes().size()) { size = meshes.meshes()[domain_id]->get_largest_dimension_size(); } else { continue; diff --git a/Studio/Visualization/Viewer.cpp b/Studio/Visualization/Viewer.cpp index 45ffb98e9b..fc3de04bb4 100644 --- a/Studio/Visualization/Viewer.cpp +++ b/Studio/Visualization/Viewer.cpp @@ -1088,6 +1088,9 @@ void Viewer::update_actors() { if (show_surface_ && meshes_.valid()) { // if (true && meshes_.valid()) { for (int i = 0; i < number_of_domains_; i++) { + if (meshes_.meshes().size() <= i) { + continue; + } renderer_->AddActor(unclipped_surface_actors_[i]); renderer_->AddActor(surface_actors_[i]); From a812e0c79e116b38b173f3203dd52c8ffd32a99d Mon Sep 17 00:00:00 2001 From: Alan Morris Date: Sat, 23 Aug 2025 12:41:05 -0600 Subject: [PATCH 03/18] Fixing error states of shared boundary (not leave project/session in inconsistent state) --- Libs/Analyze/Shape.cpp | 2 +- Libs/Groom/Groom.cpp | 87 +++++++++++++++++++++--- Libs/Mesh/MeshUtils.cpp | 26 +++++-- Libs/Mesh/MeshUtils.h | 4 +- Libs/Project/Project.h | 2 +- Libs/Python/ShapeworksPython.cpp | 4 +- Libs/Utils/Utils.cpp | 13 ++++ Libs/Utils/Utils.h | 3 + Studio/Groom/GroomTool.cpp | 7 -- Studio/Interface/ShapeWorksStudioApp.cpp | 10 ++- Testing/MeshTests/MeshTests.cpp | 2 +- 11 files changed, 130 insertions(+), 30 deletions(-) diff --git a/Libs/Analyze/Shape.cpp b/Libs/Analyze/Shape.cpp index 08a49f17b5..c029b419ce 100644 --- a/Libs/Analyze/Shape.cpp +++ b/Libs/Analyze/Shape.cpp @@ -100,7 +100,7 @@ MeshGroup Shape::get_original_meshes(bool wait) { return original_meshes_; } - if (!original_meshes_.valid() || original_meshes_.meshes().size() != subject_->get_number_of_domains()) { + if (!original_meshes_.valid() || original_meshes_.meshes().size() != subject_->get_original_filenames().size()) { original_meshes_ = MeshGroup(subject_->get_number_of_domains()); generate_meshes(subject_->get_original_filenames(), original_meshes_, true, wait); } diff --git a/Libs/Groom/Groom.cpp b/Libs/Groom/Groom.cpp index 2347643e87..42008f2128 100644 --- a/Libs/Groom/Groom.cpp +++ b/Libs/Groom/Groom.cpp @@ -7,6 +7,7 @@ #include #include #include +#include #include #include #include @@ -37,6 +38,8 @@ bool Groom::run() { auto subjects = project_->get_subjects(); + Project original_project = *project_; + if (subjects.empty()) { throw std::invalid_argument("No subjects to groom"); } @@ -93,16 +96,29 @@ bool Groom::run() { } }); - if (!run_shared_boundaries()) { + if (success && !run_shared_boundaries()) { success = false; } - if (!run_alignment()) { + if (success && !run_alignment()) { success = false; } + increment_progress(10); // alignment complete + if (!success) { + SW_DEBUG("Grooming failed, restoring original project state"); + + SW_LOG("1. Number of domains at the end: {}", project_->get_number_of_domains_per_subject()); + + // restore original project state + //*project_ = original_project; + + SW_LOG("2. Number of domains at the end: {}", project_->get_number_of_domains_per_subject()); + } project_->update_subjects(); + + SW_LOG("Number of domains at the end: {}", project_->get_number_of_domains_per_subject()); return success; } @@ -516,6 +532,7 @@ bool Groom::get_aborted() { return abort_; } //--------------------------------------------------------------------------- bool Groom::run_alignment() { size_t num_domains = project_->get_number_of_domains_per_subject(); + SW_DEBUG("Running alignment, number of domains = {}", num_domains); auto subjects = project_->get_subjects(); if (subjects.empty()) { @@ -708,7 +725,9 @@ bool Groom::run_shared_boundaries() { auto original_domain_types = project_->get_original_domain_types(); auto groomed_domain_types = project_->get_groomed_domain_types(); - for (auto shared_boundary : params.get_shared_boundaries()) { + std::atomic error_encountered = false; + + for (auto& shared_boundary : params.get_shared_boundaries()) { std::string first_domain_name = shared_boundary.first_domain; std::string second_domain_name = shared_boundary.second_domain; std::string shared_surface_name = "shared_surface_" + first_domain_name + "_" + second_domain_name; @@ -766,13 +785,28 @@ bool Groom::run_shared_boundaries() { Mesh first_mesh = get_mesh(i, first_domain, false); Mesh second_mesh = get_mesh(i, second_domain, false); - auto [extracted_l, extracted_r, extracted_s] = - MeshUtils::sharedBoundaryExtractor(first_mesh, second_mesh, shared_boundary.tolerance); + // create an empty vtk mesh + auto empty_mesh = vtkSmartPointer::New(); + + Mesh extracted_l(empty_mesh); + Mesh extracted_r(empty_mesh); + Mesh extracted_s(empty_mesh); + Mesh output_contour(empty_mesh); + try { + // returns left remainder, right remainder, and shared_surface + std::tie(extracted_l, extracted_r, extracted_s) = + MeshUtils::shared_boundary_extractor(first_mesh, second_mesh, shared_boundary.tolerance); - extracted_l.remeshPercent(.99, 1.0); - extracted_r.remeshPercent(.99, 1.0); + extracted_l.remeshPercent(.99, 1.0); + extracted_r.remeshPercent(.99, 1.0); - auto output_contour = MeshUtils::boundaryLoopExtractor(extracted_s); + output_contour = MeshUtils::extract_boundary_loop(extracted_s); + + } catch (const std::exception& e) { + SW_ERROR("Error extracting shared boundary for subject '{}': {}", subject->get_display_name(), e.what()); + // need to continue to write the empty meshes out or we will have a broken project that can no longer be saved + error_encountered = true; + } // overwrite groomed filename for first and second with extracted_l and extracted_r auto first_filename = subject->get_groomed_filenames()[first_domain]; @@ -826,6 +860,43 @@ bool Groom::run_shared_boundaries() { } }); } + + if (error_encountered) { + SW_ERROR("Errors encountered while processing shared boundaries. Please check the log for details."); + + // remove shared_surface and shared_boundary from domain_names + auto domain_names = project_->get_domain_names(); + domain_names.erase(std::remove_if(domain_names.begin(), domain_names.end(), + [](const std::string& name) { + return name.find("shared_surface_") != std::string::npos || + name.find("shared_boundary_") != std::string::npos; + }), + domain_names.end()); + project_->set_domain_names(domain_names); + + // go back and delete all the shared surfaces and boundaries that were created for all subjects + auto subjects = project_->get_subjects(); + for (auto& subject : subjects) { + auto original_filenames = subject->get_original_filenames(); + + for (auto it = original_filenames.begin(); it != original_filenames.end();) { + if (it->find("shared_surface_") != std::string::npos || it->find("shared_boundary_") != std::string::npos) { + SW_DEBUG("Removing original filename: {}", *it); + Utils::quiet_delete_file(*it); + it = original_filenames.erase(it); + } else { + ++it; + } + } + subject->set_original_filenames(original_filenames); + + subject->set_groomed_filenames({}); + subject->set_number_of_domains(domain_names.size()); + } + + return false; + } + return true; } diff --git a/Libs/Mesh/MeshUtils.cpp b/Libs/Mesh/MeshUtils.cpp index ed53a0acde..0a4b164706 100644 --- a/Libs/Mesh/MeshUtils.cpp +++ b/Libs/Mesh/MeshUtils.cpp @@ -271,13 +271,15 @@ static bool is_clockwise(const Eigen::MatrixXd& V, const Eigen::MatrixXi& F, con } //--------------------------------------------------------------------------- -Mesh MeshUtils::boundaryLoopExtractor(Mesh mesh) { +Mesh MeshUtils::extract_boundary_loop(Mesh mesh) { Eigen::MatrixXd V = mesh.points(); Eigen::MatrixXi F = mesh.faces(); std::vector > loops; igl::boundary_loop(F, loops); - assert(loops.size() == 1); + if (loops.size() != 1) { + throw std::runtime_error("Expected exactly one boundary loop in the mesh"); + } const auto& loop = loops[0]; const auto is_cw = is_clockwise(V, F, loop); @@ -365,6 +367,9 @@ static std::tuple shared_into_eigen_mesh(const static bool is_empty(const Eigen::MatrixXd& V, const Eigen::MatrixXi& F) { return V.size() == 0 || F.size() == 0; } //--------------------------------------------------------------------------- +// Identifies faces from the source mesh that lie within tolerance distance of the other mesh surface. +// Uses spatial queries to determine which faces are "shared" between the two meshes. +// Returns both the shared faces as one mesh and the remaining non-shared faces as another mesh. static std::tuple find_shared_surface( const Eigen::MatrixXd& src_V, const Eigen::MatrixXi& src_F, const Eigen::MatrixXd& other_V, const Eigen::MatrixXi& other_F, double tol = 1e-3) { @@ -408,9 +413,13 @@ static std::tuple > src_loops, shared_loops; igl::boundary_loop(src_F, src_loops); igl::boundary_loop(shared_F, shared_loops); @@ -450,7 +459,10 @@ static void move_to_boundary(const Eigen::MatrixXd& src_V, const Eigen::MatrixXi } //--------------------------------------------------------------------------- -std::array MeshUtils::sharedBoundaryExtractor(const Mesh& mesh_l, const Mesh& mesh_r, double tol) { +// Extracts shared boundary regions between two input meshes and prepares them for operations like stitching. +// Identifies shared surfaces on both meshes, then aligns the boundary of the remaining geometry to ensure +// compatibility. Returns three meshes: left remainder with aligned boundary, right remainder, and shared surface. +std::array MeshUtils::shared_boundary_extractor(const Mesh& mesh_l, const Mesh& mesh_r, double tol) { Eigen::MatrixXd V_l, V_r; Eigen::MatrixXi F_l, F_r; V_l = mesh_l.points(); @@ -473,7 +485,7 @@ std::array MeshUtils::sharedBoundaryExtractor(const Mesh& mesh_l, const Eigen::MatrixXd bridge_V; Eigen::MatrixXi bridge_F; - move_to_boundary(rem_V_l, rem_F_l, shared_V_r, shared_F_r, bridge_V, bridge_F); + snap_boundary_to_mesh(rem_V_l, rem_F_l, shared_V_r, shared_F_r, bridge_V, bridge_F); Mesh out_l(bridge_V, bridge_F); Mesh out_r(rem_V_r, rem_F_r); diff --git a/Libs/Mesh/MeshUtils.h b/Libs/Mesh/MeshUtils.h index 26e302d59a..ab991ae8d6 100644 --- a/Libs/Mesh/MeshUtils.h +++ b/Libs/Mesh/MeshUtils.h @@ -47,10 +47,10 @@ class MeshUtils { static int findReferenceMesh(std::vector& meshes, int random_subset_size = -1); /// boundary loop extractor for a given mesh - static Mesh boundaryLoopExtractor(Mesh mesh); + static Mesh extract_boundary_loop(Mesh mesh); /// shared boundary extractor for the left and right mesh - static std::array sharedBoundaryExtractor(const Mesh& mesh_l, const Mesh& mesh_r, double tol); + static std::array shared_boundary_extractor(const Mesh& mesh_l, const Mesh& mesh_r, double tol); /// generates and adds normals for points and faces for each mesh in given set of meshes static void generateNormals(const std::vector >& meshes, bool forceRegen = false); diff --git a/Libs/Project/Project.h b/Libs/Project/Project.h index f773dcc5cd..a2bd755695 100644 --- a/Libs/Project/Project.h +++ b/Libs/Project/Project.h @@ -196,7 +196,7 @@ class Project { // map of type (e.g. groom, optimize) to map (domain->Parameters) std::map> parameters_; - const int supported_version_{2}; + static constexpr int supported_version_{2}; int version_{2}; }; } // namespace shapeworks diff --git a/Libs/Python/ShapeworksPython.cpp b/Libs/Python/ShapeworksPython.cpp index 6ce5b7e81a..a68df59df6 100644 --- a/Libs/Python/ShapeworksPython.cpp +++ b/Libs/Python/ShapeworksPython.cpp @@ -1187,13 +1187,13 @@ PYBIND11_MODULE(shapeworks_py, m) { .def_static("findReferenceMesh", &MeshUtils::findReferenceMesh, "find reference mesh from a set of meshes", "meshes"_a, "random_subset"_a = -1) - .def_static("boundaryLoopExtractor", &MeshUtils::boundaryLoopExtractor, + .def_static("boundaryLoopExtractor", &MeshUtils::extract_boundary_loop, "for a mesh extracts the boundary loop and export the boundary loop as a contour .vtp file", "mesh"_a) .def_static( "sharedBoundaryExtractor", [](const Mesh& mesh_l, const Mesh& mesh_r, float tol) -> decltype(auto) { - std::array output = MeshUtils::sharedBoundaryExtractor(mesh_l, mesh_r, tol); + std::array output = MeshUtils::shared_boundary_extractor(mesh_l, mesh_r, tol); // std::move passes ownership to Python return py::make_tuple(std::move(output[0]), std::move(output[1]), std::move(output[2])); diff --git a/Libs/Utils/Utils.cpp b/Libs/Utils/Utils.cpp index 813b79c811..fa6a0d7cab 100644 --- a/Libs/Utils/Utils.cpp +++ b/Libs/Utils/Utils.cpp @@ -181,6 +181,19 @@ void Utils::writeParticleIds(char* filename, std::vector ids) ofs.close(); } +//--------------------------------------------------------------------- +// quietly delete a file, if it doesn't exist, no error, if it fails, no error +void Utils::quiet_delete_file(const std::string &filename) +{ + try { + if (std::remove(filename.c_str()) != 0) { + // file deletion failed, but we don't care + } + } catch (...) { + // no error message, just quietly ignore it. + } +} + //--------------- point cloud queries -------------------------------- void Utils::computeCenterOfMassForShapeEnsemble (std::vector< std::vector< itk::Point< double, 3 > > > points_list, itk::Point< double, 3 > & center) { diff --git a/Libs/Utils/Utils.h b/Libs/Utils/Utils.h index afa6f53790..5835639df8 100644 --- a/Libs/Utils/Utils.h +++ b/Libs/Utils/Utils.h @@ -66,6 +66,9 @@ class Utils{ static std::vector readParticleIds(char* filename); static void writeParticleIds(char* filename, std::vector ids); + + static void quiet_delete_file(const std::string& filename); + //--------------- point cloud queries -------------------------------- static void computeCenterOfMassForShapeEnsemble (std::vector< std::vector< itk::Point< double, 3 > > > points_list, itk::Point< double, 3 > & center); static void computeCenterOfMassForShape (std::vector< itk::Point< double, 3 > > points, itk::Point< double, 3 > & center); diff --git a/Studio/Groom/GroomTool.cpp b/Studio/Groom/GroomTool.cpp index 71526f032a..9369387ec0 100644 --- a/Studio/Groom/GroomTool.cpp +++ b/Studio/Groom/GroomTool.cpp @@ -221,14 +221,9 @@ void GroomTool::delete_shared_boundary_clicked() { } for (int i = selected_rows.size() - 1; i >= 0; i--) { - std::cerr << "erasing shared boundary: " << shared_boundaries[selected_rows[i].row()].first_domain << " <-> " - << shared_boundaries[selected_rows[i].row()].second_domain << " with tolerance " - << shared_boundaries[selected_rows[i].row()].tolerance << "\n"; shared_boundaries.erase(shared_boundaries.begin() + selected_rows[i].row()); } - std::cerr << "There are now " << shared_boundaries.size() << " shared boundaries left.\n"; - params.set_shared_boundaries(shared_boundaries); params.save_to_project(); update_shared_boundary_table(); @@ -326,8 +321,6 @@ void GroomTool::update_shared_boundary_table() { auto params = GroomParameters(session_->get_project(), ""); auto shared_boundaries = params.get_shared_boundaries(); - std::cerr << "There are " << shared_boundaries.size() << " shared boundaries.\n"; - QStringList table_headers; table_headers << "First Domain"; table_headers << "Second Domain"; diff --git a/Studio/Interface/ShapeWorksStudioApp.cpp b/Studio/Interface/ShapeWorksStudioApp.cpp index 046eb7f884..9c5d4426f6 100644 --- a/Studio/Interface/ShapeWorksStudioApp.cpp +++ b/Studio/Interface/ShapeWorksStudioApp.cpp @@ -1318,8 +1318,16 @@ void ShapeWorksStudioApp::handle_groom_start() { //--------------------------------------------------------------------------- void ShapeWorksStudioApp::handle_groom_complete() { update_view_combo(); - ui_->view_mode_combobox->setCurrentIndex(DisplayMode::Groomed); + + if (!session_->groomed_present()) { + // grooming may have failed, if so, we don't want to switch to groomed that don't exist + session_->set_display_mode(DisplayMode::Original); + } else { + session_->set_display_mode(DisplayMode::Groomed); + } + session_->handle_clear_cache(); + update_display(true); visualizer_->reset_camera(); enable_possible_actions(); diff --git a/Testing/MeshTests/MeshTests.cpp b/Testing/MeshTests/MeshTests.cpp index 07bd92c9d2..355776c25d 100644 --- a/Testing/MeshTests/MeshTests.cpp +++ b/Testing/MeshTests/MeshTests.cpp @@ -799,7 +799,7 @@ TEST(MeshTests, sharedBoundaryExtractor) { TEST(MeshTests, boundaryLoopExtractor) { Mesh ground_truth(std::string(TEST_DATA_DIR) + "/shared_boundary/00_out_c.vtp"); Mesh mesh(std::string(TEST_DATA_DIR) + "/shared_boundary/00_out_s.vtk"); - Mesh loop = MeshUtils::boundaryLoopExtractor(mesh); + Mesh loop = MeshUtils::extract_boundary_loop(mesh); ASSERT_TRUE(loop == ground_truth); } From 7fe69e38a44c0e94843e0aa6cb8ae5ebe8cce589 Mon Sep 17 00:00:00 2001 From: Alan Morris Date: Mon, 25 Aug 2025 10:49:05 -0600 Subject: [PATCH 04/18] Replace asserts with exceptions for release handling --- Libs/Groom/Groom.cpp | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/Libs/Groom/Groom.cpp b/Libs/Groom/Groom.cpp index 42008f2128..a2693c07bf 100644 --- a/Libs/Groom/Groom.cpp +++ b/Libs/Groom/Groom.cpp @@ -1001,8 +1001,13 @@ std::string Groom::get_output_filename(std::string input, DomainType domain_type //--------------------------------------------------------------------------- Mesh Groom::get_mesh(int subject, int domain, bool transformed) { auto subjects = project_->get_subjects(); - assert(subject < subjects.size()); - assert(domain < subjects[subject]->get_original_filenames().size()); + if (subject >= subjects.size()) { + throw std::out_of_range("subject index out of range"); + } + if (domain >= subjects[subject]->get_original_filenames().size()) { + throw std::out_of_range("domain index out of range"); + } + auto path = subjects[subject]->get_original_filenames()[domain]; auto constraint_filename = subjects[subject]->get_constraints_filenames(); From 58a8b71f8b040b462d3f253f9e3d61d309f05aff Mon Sep 17 00:00:00 2001 From: Alan Morris Date: Mon, 25 Aug 2025 10:50:33 -0600 Subject: [PATCH 05/18] Multiple shared boundaries handling --- Libs/Groom/Groom.cpp | 256 +++++++++++++++++++++++++++++++++- Libs/Groom/Groom.h | 3 + Libs/Mesh/MeshUtils.cpp | 8 +- Libs/Project/Project.cpp | 10 +- Libs/Project/ProjectUtils.cpp | 2 +- 5 files changed, 266 insertions(+), 13 deletions(-) diff --git a/Libs/Groom/Groom.cpp b/Libs/Groom/Groom.cpp index a2693c07bf..a22be35dd5 100644 --- a/Libs/Groom/Groom.cpp +++ b/Libs/Groom/Groom.cpp @@ -19,14 +19,16 @@ using namespace shapeworks; +namespace { // for concurrent access static std::mutex mutex; +} // namespace using PixelType = float; using ImageType = itk::Image; //--------------------------------------------------------------------------- -Groom::Groom(ProjectHandle project) { project_ = project; } +Groom::Groom(ProjectHandle project) : project_{project} {} //--------------------------------------------------------------------------- bool Groom::run() { @@ -715,6 +717,10 @@ bool Groom::run_alignment() { bool Groom::run_shared_boundaries() { // only need to check on one domain auto params = GroomParameters(project_, project_->get_domain_names()[0]); + + // first, we must remove any existing shared surfaces and boundaries from the project + clear_unused_shared_boundaries(); + if (!params.get_shared_boundaries_enabled()) { return true; } @@ -727,7 +733,7 @@ bool Groom::run_shared_boundaries() { std::atomic error_encountered = false; - for (auto& shared_boundary : params.get_shared_boundaries()) { + for (const auto& shared_boundary : params.get_shared_boundaries()) { std::string first_domain_name = shared_boundary.first_domain; std::string second_domain_name = shared_boundary.second_domain; std::string shared_surface_name = "shared_surface_" + first_domain_name + "_" + second_domain_name; @@ -900,6 +906,252 @@ bool Groom::run_shared_boundaries() { return true; } +/* +//--------------------------------------------------------------------------- +void Groom::clear_unused_shared_boundaries() { + auto params = GroomParameters(project_, project_->get_domain_names()[0]); + + std::vector shared_surfaces_to_keep; + std::vector shared_boundaries_to_keep; + + for (const auto& shared_boundary : params.get_shared_boundaries()) { + std::string first_domain_name = shared_boundary.first_domain; + std::string second_domain_name = shared_boundary.second_domain; + std::string shared_surface_name = "shared_surface_" + first_domain_name + "_" + second_domain_name; + std::string shared_boundary_name = "shared_boundary_" + first_domain_name + "_" + second_domain_name; + shared_surfaces_to_keep.push_back(shared_surface_name); + shared_boundaries_to_keep.push_back(shared_boundary_name); + } + + // Lambda to check if a name is a shared domain that should be removed + auto is_shared_domain_to_remove = [&shared_surfaces_to_keep, &shared_boundaries_to_keep](const std::string& name) { + bool is_shared = + name.find("shared_surface_") != std::string::npos || name.find("shared_boundary_") != std::string::npos; + + if (!is_shared) { + return false; // Not a shared domain at all + } + + // Check if it's in either keep list + bool should_keep = std::find(shared_surfaces_to_keep.begin(), shared_surfaces_to_keep.end(), name) != + shared_surfaces_to_keep.end() || + std::find(shared_boundaries_to_keep.begin(), shared_boundaries_to_keep.end(), name) != + shared_boundaries_to_keep.end(); + + return !should_keep; // Remove if it's shared but not in keep lists + }; + + // Lambda to remove shared files from a vector and delete them from disk + auto remove_shared_files = [&is_shared_domain_to_remove](std::vector& filenames) { + for (auto it = filenames.begin(); it != filenames.end();) { + if (is_shared_domain_to_remove(*it)) { + Utils::quiet_delete_file(*it); + it = filenames.erase(it); + } else { + ++it; + } + } + }; + + auto domain_names = project_->get_domain_names(); + + // Remove shared_surface and shared_boundary from domain_names + domain_names.erase(std::remove_if(domain_names.begin(), domain_names.end(), is_shared_domain_to_remove), + domain_names.end()); + + // Find indices to remove (in reverse order to avoid index shifting issues) + std::vector indices_to_remove; + const auto& original_domain_names = project_->get_domain_names(); + for (size_t i = 0; i < original_domain_names.size(); ++i) { + if (is_shared_domain_to_remove(original_domain_names[i])) { + indices_to_remove.push_back(i); + std::cerr << "planning to remove domain index: " << i << " with name: " << original_domain_names[i] << "\n"; + } + } + + // Remove corresponding entries from domain types arrays + auto original_domain_types = project_->get_original_domain_types(); + auto groomed_domain_types = project_->get_groomed_domain_types(); + + for (auto it = indices_to_remove.rbegin(); it != indices_to_remove.rend(); ++it) { + if (*it < original_domain_types.size()) { + original_domain_types.erase(original_domain_types.begin() + *it); + } + if (*it < groomed_domain_types.size()) { + groomed_domain_types.erase(groomed_domain_types.begin() + *it); + } + } + + project_->set_domain_names(domain_names); + project_->set_original_domain_types(original_domain_types); + project_->set_groomed_domain_types(groomed_domain_types); + + // Clean up files from all subjects + auto subjects = project_->get_subjects(); + for (auto& subject : subjects) { + // Remove shared surface/boundary files using lambda + auto original_filenames = subject->get_original_filenames(); + remove_shared_files(original_filenames); + subject->set_original_filenames(original_filenames); + + auto groomed_filenames = subject->get_groomed_filenames(); + remove_shared_files(groomed_filenames); + subject->set_groomed_filenames(groomed_filenames); + + // Remove alignment transforms for removed domains + auto transforms = subject->get_groomed_transforms(); + for (auto it = indices_to_remove.rbegin(); it != indices_to_remove.rend(); ++it) { + if (*it < transforms.size()) { + transforms.erase(transforms.begin() + *it); + } + } + subject->set_groomed_transforms(transforms); + + // Clean up particle files using lambda + auto local_particles = subject->get_local_particle_filenames(); + remove_shared_files(local_particles); + subject->set_local_particle_filenames(local_particles); + + auto world_particles = subject->get_world_particle_filenames(); + remove_shared_files(world_particles); + subject->set_world_particle_filenames(world_particles); + + // procrustes + auto procrustes = subject->get_procrustes_transforms(); + for (auto it = indices_to_remove.rbegin(); it != indices_to_remove.rend(); ++it) { + if (*it < procrustes.size()) { + procrustes.erase(procrustes.begin() + *it); + } + } + subject->set_procrustes_transforms(procrustes); + + subject->set_number_of_domains(domain_names.size()); + } +} +*/ + +//--------------------------------------------------------------------------- +void Groom::clear_unused_shared_boundaries() { + auto params = GroomParameters(project_, project_->get_domain_names()[0]); + + std::vector shared_surfaces_to_keep; + std::vector shared_boundaries_to_keep; + + if (params.get_shared_boundaries_enabled()) { + for (const auto& shared_boundary : params.get_shared_boundaries()) { + std::string first_domain_name = shared_boundary.first_domain; + std::string second_domain_name = shared_boundary.second_domain; + std::string shared_surface_name = "shared_surface_" + first_domain_name + "_" + second_domain_name; + std::string shared_boundary_name = "shared_boundary_" + first_domain_name + "_" + second_domain_name; + shared_surfaces_to_keep.push_back(shared_surface_name); + shared_boundaries_to_keep.push_back(shared_boundary_name); + } + } + + // Lambda to check if a name is a shared domain that should be removed + auto is_shared_domain_to_remove = [&shared_surfaces_to_keep, &shared_boundaries_to_keep](const std::string& name) { + // name may be dir/file and we only care about file + std::string filename = StringUtils::getFilename(name); + + bool is_shared = + filename.find("shared_surface_") != std::string::npos || filename.find("shared_boundary_") != std::string::npos; + + if (!is_shared) { + return false; // Not a shared domain at all + } + + // Check if it's in either keep list + bool should_keep = std::find(shared_surfaces_to_keep.begin(), shared_surfaces_to_keep.end(), name) != + shared_surfaces_to_keep.end() || + std::find(shared_boundaries_to_keep.begin(), shared_boundaries_to_keep.end(), name) != + shared_boundaries_to_keep.end(); + + return !should_keep; // Remove if it's shared but not in keep lists + }; + + // Lambda to remove shared files from a vector and delete them from disk + auto remove_shared_files = [&is_shared_domain_to_remove](std::vector& filenames) { + for (auto it = filenames.begin(); it != filenames.end();) { + if (is_shared_domain_to_remove(*it)) { + std::cerr << "Removing file: " << *it << "\n"; + Utils::quiet_delete_file(*it); + it = filenames.erase(it); + } else { + ++it; + } + } + }; + + // Lambda to remove elements from vector by matching domain names + auto remove_by_domain_names = [&is_shared_domain_to_remove](auto& vec, const std::vector& domain_names) { + for (auto it = vec.begin(); it != vec.end();) { + size_t index = std::distance(vec.begin(), it); + if (index < domain_names.size() && is_shared_domain_to_remove(domain_names[index])) { + it = vec.erase(it); + } else { + ++it; + } + } + }; + + auto domain_names = project_->get_domain_names(); + auto original_domain_names = domain_names; // Save original for reference + + // Remove shared domains from domain_names + domain_names.erase(std::remove_if(domain_names.begin(), domain_names.end(), is_shared_domain_to_remove), + domain_names.end()); + + // Remove corresponding entries from domain types arrays using original domain names as reference + auto original_domain_types = project_->get_original_domain_types(); + auto groomed_domain_types = project_->get_groomed_domain_types(); + + remove_by_domain_names(original_domain_types, original_domain_names); + remove_by_domain_names(groomed_domain_types, original_domain_names); + + project_->set_domain_names(domain_names); + project_->set_original_domain_types(original_domain_types); + project_->set_groomed_domain_types(groomed_domain_types); + + // Clean up files from all subjects + auto subjects = project_->get_subjects(); + for (auto& subject : subjects) { + // Remove shared surface/boundary files using lambda + auto original_filenames = subject->get_original_filenames(); + remove_shared_files(original_filenames); + subject->set_original_filenames(original_filenames); + + auto groomed_filenames = subject->get_groomed_filenames(); + remove_shared_files(groomed_filenames); + subject->set_groomed_filenames(groomed_filenames); + + // Remove alignment transforms for removed domains + auto transforms = subject->get_groomed_transforms(); + remove_by_domain_names(transforms, original_domain_names); + subject->set_groomed_transforms(transforms); + + // Clean up particle files using lambda + auto local_particles = subject->get_local_particle_filenames(); + remove_shared_files(local_particles); + subject->set_local_particle_filenames(local_particles); + + auto world_particles = subject->get_world_particle_filenames(); + remove_shared_files(world_particles); + subject->set_world_particle_filenames(world_particles); + + // procrustes transforms + auto procrustes = subject->get_procrustes_transforms(); + remove_by_domain_names(procrustes, original_domain_names); + subject->set_procrustes_transforms(procrustes); + + subject->set_number_of_domains(domain_names.size()); + } + + std::cerr << "Final domain names after cleanup:\n"; + for (const auto& name : domain_names) { + std::cerr << " - " << name << "\n"; + } + project_->set_domain_names(domain_names); +} //--------------------------------------------------------------------------- void Groom::assign_transforms(std::vector> transforms, int domain, bool global) { auto subjects = project_->get_subjects(); diff --git a/Libs/Groom/Groom.h b/Libs/Groom/Groom.h index 22e3e6882a..6d73225e4e 100644 --- a/Libs/Groom/Groom.h +++ b/Libs/Groom/Groom.h @@ -66,6 +66,9 @@ class Groom { //! Create shared boundary surface and contour if requested bool run_shared_boundaries(); + //! Remove all shared boundaries + void clear_unused_shared_boundaries(); + void assign_transforms(std::vector> transforms, int domain, bool global = false); static std::vector> get_icp_transforms(const std::vector meshes, Mesh reference); diff --git a/Libs/Mesh/MeshUtils.cpp b/Libs/Mesh/MeshUtils.cpp index 0a4b164706..f7e35ed87e 100644 --- a/Libs/Mesh/MeshUtils.cpp +++ b/Libs/Mesh/MeshUtils.cpp @@ -277,8 +277,8 @@ Mesh MeshUtils::extract_boundary_loop(Mesh mesh) { std::vector > loops; igl::boundary_loop(F, loops); - if (loops.size() != 1) { - throw std::runtime_error("Expected exactly one boundary loop in the mesh"); + if (loops.size() < 1) { + throw std::runtime_error("Expected at least one boundary loop in the mesh"); } const auto& loop = loops[0]; @@ -424,9 +424,9 @@ static void snap_boundary_to_mesh(const Eigen::MatrixXd& src_V, const Eigen::Mat igl::boundary_loop(src_F, src_loops); igl::boundary_loop(shared_F, shared_loops); - if (src_loops.size() != 1 || shared_loops.size() != 1) { + if (src_loops.size() < 1 || shared_loops.size() < 1) { SW_DEBUG("src_loops size: {}, shared_loops size: {}", src_loops.size(), shared_loops.size()); - throw std::runtime_error("Expected exactly one boundary loop in each mesh"); + throw std::runtime_error("Expected at least one boundary loop in each mesh"); } const auto& src_loop = src_loops[0]; diff --git a/Libs/Project/Project.cpp b/Libs/Project/Project.cpp index 848c1b1157..e963cdd999 100644 --- a/Libs/Project/Project.cpp +++ b/Libs/Project/Project.cpp @@ -192,16 +192,14 @@ void Project::set_subjects(const std::vector>& subjects //--------------------------------------------------------------------------- void Project::update_subjects() { + originals_present_ = false; + groomed_present_ = false; + particles_present_ = false; + if (subjects_.empty()) { - originals_present_ = false; - groomed_present_ = false; - particles_present_ = false; return; } - originals_present_ = false; - groomed_present_ = false; - particles_present_ = false; images_present_ = false; auto groomed_subject = subjects_[0]; diff --git a/Libs/Project/ProjectUtils.cpp b/Libs/Project/ProjectUtils.cpp index d858abddb2..5b4d4ed97c 100644 --- a/Libs/Project/ProjectUtils.cpp +++ b/Libs/Project/ProjectUtils.cpp @@ -318,7 +318,7 @@ static void assign_transforms(StringMap& j, std::string prefix, std::vector Date: Tue, 26 Aug 2025 12:44:59 -0600 Subject: [PATCH 06/18] More fixes for multiple shared boundaries --- Libs/Groom/Groom.cpp | 205 ++++++++----------------------------- Libs/Groom/Groom.h | 5 +- Studio/Groom/GroomTool.cpp | 35 +++++-- Studio/Groom/GroomTool.ui | 35 +++---- 4 files changed, 82 insertions(+), 198 deletions(-) diff --git a/Libs/Groom/Groom.cpp b/Libs/Groom/Groom.cpp index a22be35dd5..e2f05502ad 100644 --- a/Libs/Groom/Groom.cpp +++ b/Libs/Groom/Groom.cpp @@ -40,12 +40,15 @@ bool Groom::run() { auto subjects = project_->get_subjects(); - Project original_project = *project_; - if (subjects.empty()) { throw std::invalid_argument("No subjects to groom"); } + // clear alignment transforms + for (auto& subject : subjects) { + subject->set_groomed_transforms({}); + } + total_ops_ = get_total_ops(); std::atomic success = true; @@ -57,9 +60,13 @@ bool Groom::run() { success = false; continue; } + + auto domain_name = project_->get_domain_names()[domain]; // skip "shared_surface" and "shared_boundary" - if (project_->get_domain_names()[domain] == "shared_surface" || - project_->get_domain_names()[domain] == "shared_boundary") { + bool is_shared_boundary = domain_name.find("shared_surface_") != std::string::npos || + domain_name.find("shared_boundary_") != std::string::npos; + if (is_shared_boundary) { + SW_DEBUG("Skipping shared boundary domain '{}' during grooming", domain_name); continue; } @@ -108,19 +115,8 @@ bool Groom::run() { increment_progress(10); // alignment complete - if (!success) { - SW_DEBUG("Grooming failed, restoring original project state"); - - SW_LOG("1. Number of domains at the end: {}", project_->get_number_of_domains_per_subject()); - - // restore original project state - //*project_ = original_project; - - SW_LOG("2. Number of domains at the end: {}", project_->get_number_of_domains_per_subject()); - } project_->update_subjects(); - SW_LOG("Number of domains at the end: {}", project_->get_number_of_domains_per_subject()); return success; } @@ -131,6 +127,8 @@ bool Groom::image_pipeline(std::shared_ptr subject, size_t domain) { auto original = subject->get_original_filenames()[domain]; + SW_DEBUG("Grooming image: {}", original); + // load the image Image image(original); @@ -312,6 +310,7 @@ bool Groom::mesh_pipeline(std::shared_ptr subject, size_t domain) { auto params = GroomParameters(project_, project_->get_domain_names()[domain]); auto original = subject->get_original_filenames()[domain]; + SW_DEBUG("Grooming mesh: {}", original); // groomed mesh name std::string groom_name = get_output_filename(original, DomainType::Mesh); @@ -375,6 +374,11 @@ bool Groom::run_mesh_pipeline(Mesh& mesh, GroomParameters params) { 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(); @@ -788,8 +792,8 @@ bool Groom::run_shared_boundaries() { auto subject = subjects[i]; - Mesh first_mesh = get_mesh(i, first_domain, false); - Mesh second_mesh = get_mesh(i, second_domain, false); + Mesh first_mesh = get_mesh(i, first_domain, false, MeshSource::Groomed); + Mesh second_mesh = get_mesh(i, second_domain, false, MeshSource::Groomed); // create an empty vtk mesh auto empty_mesh = vtkSmartPointer::New(); @@ -823,8 +827,8 @@ bool Groom::run_shared_boundaries() { MeshUtils::threadSafeWriteMesh(first_filename, extracted_l); MeshUtils::threadSafeWriteMesh(second_filename, extracted_r); - auto shared_surface_filename = get_output_filename("shared_surface", DomainType::Mesh); - auto shared_boundary_filename = get_output_filename("shared_boundary", DomainType::Contour); + auto shared_surface_filename = get_output_filename(shared_surface_name, DomainType::Mesh); + auto shared_boundary_filename = get_output_filename(shared_boundary_name, DomainType::Contour); extracted_s.write(shared_surface_filename); output_contour.write(shared_boundary_filename); @@ -906,130 +910,6 @@ bool Groom::run_shared_boundaries() { return true; } -/* -//--------------------------------------------------------------------------- -void Groom::clear_unused_shared_boundaries() { - auto params = GroomParameters(project_, project_->get_domain_names()[0]); - - std::vector shared_surfaces_to_keep; - std::vector shared_boundaries_to_keep; - - for (const auto& shared_boundary : params.get_shared_boundaries()) { - std::string first_domain_name = shared_boundary.first_domain; - std::string second_domain_name = shared_boundary.second_domain; - std::string shared_surface_name = "shared_surface_" + first_domain_name + "_" + second_domain_name; - std::string shared_boundary_name = "shared_boundary_" + first_domain_name + "_" + second_domain_name; - shared_surfaces_to_keep.push_back(shared_surface_name); - shared_boundaries_to_keep.push_back(shared_boundary_name); - } - - // Lambda to check if a name is a shared domain that should be removed - auto is_shared_domain_to_remove = [&shared_surfaces_to_keep, &shared_boundaries_to_keep](const std::string& name) { - bool is_shared = - name.find("shared_surface_") != std::string::npos || name.find("shared_boundary_") != std::string::npos; - - if (!is_shared) { - return false; // Not a shared domain at all - } - - // Check if it's in either keep list - bool should_keep = std::find(shared_surfaces_to_keep.begin(), shared_surfaces_to_keep.end(), name) != - shared_surfaces_to_keep.end() || - std::find(shared_boundaries_to_keep.begin(), shared_boundaries_to_keep.end(), name) != - shared_boundaries_to_keep.end(); - - return !should_keep; // Remove if it's shared but not in keep lists - }; - - // Lambda to remove shared files from a vector and delete them from disk - auto remove_shared_files = [&is_shared_domain_to_remove](std::vector& filenames) { - for (auto it = filenames.begin(); it != filenames.end();) { - if (is_shared_domain_to_remove(*it)) { - Utils::quiet_delete_file(*it); - it = filenames.erase(it); - } else { - ++it; - } - } - }; - - auto domain_names = project_->get_domain_names(); - - // Remove shared_surface and shared_boundary from domain_names - domain_names.erase(std::remove_if(domain_names.begin(), domain_names.end(), is_shared_domain_to_remove), - domain_names.end()); - - // Find indices to remove (in reverse order to avoid index shifting issues) - std::vector indices_to_remove; - const auto& original_domain_names = project_->get_domain_names(); - for (size_t i = 0; i < original_domain_names.size(); ++i) { - if (is_shared_domain_to_remove(original_domain_names[i])) { - indices_to_remove.push_back(i); - std::cerr << "planning to remove domain index: " << i << " with name: " << original_domain_names[i] << "\n"; - } - } - - // Remove corresponding entries from domain types arrays - auto original_domain_types = project_->get_original_domain_types(); - auto groomed_domain_types = project_->get_groomed_domain_types(); - - for (auto it = indices_to_remove.rbegin(); it != indices_to_remove.rend(); ++it) { - if (*it < original_domain_types.size()) { - original_domain_types.erase(original_domain_types.begin() + *it); - } - if (*it < groomed_domain_types.size()) { - groomed_domain_types.erase(groomed_domain_types.begin() + *it); - } - } - - project_->set_domain_names(domain_names); - project_->set_original_domain_types(original_domain_types); - project_->set_groomed_domain_types(groomed_domain_types); - - // Clean up files from all subjects - auto subjects = project_->get_subjects(); - for (auto& subject : subjects) { - // Remove shared surface/boundary files using lambda - auto original_filenames = subject->get_original_filenames(); - remove_shared_files(original_filenames); - subject->set_original_filenames(original_filenames); - - auto groomed_filenames = subject->get_groomed_filenames(); - remove_shared_files(groomed_filenames); - subject->set_groomed_filenames(groomed_filenames); - - // Remove alignment transforms for removed domains - auto transforms = subject->get_groomed_transforms(); - for (auto it = indices_to_remove.rbegin(); it != indices_to_remove.rend(); ++it) { - if (*it < transforms.size()) { - transforms.erase(transforms.begin() + *it); - } - } - subject->set_groomed_transforms(transforms); - - // Clean up particle files using lambda - auto local_particles = subject->get_local_particle_filenames(); - remove_shared_files(local_particles); - subject->set_local_particle_filenames(local_particles); - - auto world_particles = subject->get_world_particle_filenames(); - remove_shared_files(world_particles); - subject->set_world_particle_filenames(world_particles); - - // procrustes - auto procrustes = subject->get_procrustes_transforms(); - for (auto it = indices_to_remove.rbegin(); it != indices_to_remove.rend(); ++it) { - if (*it < procrustes.size()) { - procrustes.erase(procrustes.begin() + *it); - } - } - subject->set_procrustes_transforms(procrustes); - - subject->set_number_of_domains(domain_names.size()); - } -} -*/ - //--------------------------------------------------------------------------- void Groom::clear_unused_shared_boundaries() { auto params = GroomParameters(project_, project_->get_domain_names()[0]); @@ -1129,27 +1009,13 @@ void Groom::clear_unused_shared_boundaries() { remove_by_domain_names(transforms, original_domain_names); subject->set_groomed_transforms(transforms); - // Clean up particle files using lambda - auto local_particles = subject->get_local_particle_filenames(); - remove_shared_files(local_particles); - subject->set_local_particle_filenames(local_particles); - - auto world_particles = subject->get_world_particle_filenames(); - remove_shared_files(world_particles); - subject->set_world_particle_filenames(world_particles); - - // procrustes transforms - auto procrustes = subject->get_procrustes_transforms(); - remove_by_domain_names(procrustes, original_domain_names); - subject->set_procrustes_transforms(procrustes); + subject->set_local_particle_filenames({}); + subject->set_world_particle_filenames({}); + subject->set_procrustes_transforms({}); subject->set_number_of_domains(domain_names.size()); } - std::cerr << "Final domain names after cleanup:\n"; - for (const auto& name : domain_names) { - std::cerr << " - " << name << "\n"; - } project_->set_domain_names(domain_names); } //--------------------------------------------------------------------------- @@ -1251,16 +1117,25 @@ std::string Groom::get_output_filename(std::string input, DomainType domain_type } //--------------------------------------------------------------------------- -Mesh Groom::get_mesh(int subject, int domain, bool transformed) { +Mesh Groom::get_mesh(int subject, int domain, bool transformed, MeshSource source) { auto subjects = project_->get_subjects(); if (subject >= subjects.size()) { throw std::out_of_range("subject index out of range"); } - if (domain >= subjects[subject]->get_original_filenames().size()) { - throw std::out_of_range("domain index out of range"); - } - auto path = subjects[subject]->get_original_filenames()[domain]; + std::string path; + + if (source == MeshSource::Original) { + if (domain >= subjects[subject]->get_original_filenames().size()) { + throw std::out_of_range("domain index out of range"); + } + path = subjects[subject]->get_original_filenames()[domain]; + } else { + if (domain >= subjects[subject]->get_groomed_filenames().size()) { + throw std::out_of_range("domain index out of range"); + } + path = subjects[subject]->get_groomed_filenames()[domain]; + } auto constraint_filename = subjects[subject]->get_constraints_filenames(); Constraints constraint; diff --git a/Libs/Groom/Groom.h b/Libs/Groom/Groom.h index 6d73225e4e..f1f4174d86 100644 --- a/Libs/Groom/Groom.h +++ b/Libs/Groom/Groom.h @@ -39,6 +39,8 @@ class Groom { std::atomic progress_counter_ = 0; private: + enum class MeshSource { Original, Groomed }; + //! Return the number of operations that will be performed int get_total_ops(); @@ -83,7 +85,8 @@ class Groom { std::vector> get_combined_points(); - Mesh get_mesh(int subject, int domain, bool transformed = false); + Mesh get_mesh(int subject, int domain, bool transformed = false, MeshSource source = MeshSource::Original); + vtkSmartPointer get_landmarks(int subject, int domain); diff --git a/Studio/Groom/GroomTool.cpp b/Studio/Groom/GroomTool.cpp index 9369387ec0..98ec69854b 100644 --- a/Studio/Groom/GroomTool.cpp +++ b/Studio/Groom/GroomTool.cpp @@ -5,6 +5,7 @@ #include #include #include +#include #include #include @@ -116,13 +117,13 @@ GroomTool::GroomTool(Preferences& prefs, Telemetry& telemetry) : preferences_(pr QIntValidator* above_zero = new QIntValidator(1, std::numeric_limits::max(), this); QIntValidator* zero_and_up = new QIntValidator(0, std::numeric_limits::max(), this); - QDoubleValidator* double_validator = new QDoubleValidator(0, std::numeric_limits::max(), 1000, this); ui_->laplacian_iterations->setValidator(zero_and_up); ui_->laplacian_relaxation->setValidator(double_validator); ui_->sinc_iterations->setValidator(zero_and_up); ui_->sinc_passband->setValidator(double_validator); + ui_->shared_boundary_tolerance->setValidator(double_validator); auto line_edits = findChildren(); for (auto line_edit : line_edits) { @@ -145,6 +146,9 @@ GroomTool::GroomTool(Preferences& prefs, Telemetry& telemetry) : preferences_(pr connect(combo_box, qOverload(&QComboBox::currentIndexChanged), this, &GroomTool::set_session_modified); } + // read only + ui_->shared_boundary_table->setEditTriggers(QAbstractItemView::NoEditTriggers); + update_ui(); } @@ -185,7 +189,7 @@ void GroomTool::add_shared_boundary_clicked() { GroomParameters::SharedBoundary boundary; boundary.first_domain = ui_->shared_boundary_first_domain->currentText().toStdString(); boundary.second_domain = ui_->shared_boundary_second_domain->currentText().toStdString(); - boundary.tolerance = ui_->shared_boundary_tolerance->value(); + boundary.tolerance = ui_->shared_boundary_tolerance->text().toDouble(); if (boundary.first_domain == boundary.second_domain) { SW_ERROR("Cannot add a shared boundary between the same domain"); return; @@ -320,21 +324,17 @@ void GroomTool::apply_to_all_domains_changed() { void GroomTool::update_shared_boundary_table() { auto params = GroomParameters(session_->get_project(), ""); auto shared_boundaries = params.get_shared_boundaries(); - QStringList table_headers; table_headers << "First Domain"; table_headers << "Second Domain"; table_headers << "Tolerance"; - QTableWidget* table = ui_->shared_boundary_table; table->clear(); table->setRowCount(shared_boundaries.size()); table->setColumnCount(table_headers.size()); - table->setHorizontalHeaderLabels(table_headers); table->horizontalHeader()->setSectionResizeMode(QHeaderView::Stretch); table->setSelectionBehavior(QAbstractItemView::SelectRows); - int row = 0; for (const auto& boundary : shared_boundaries) { table->setItem(row, 0, new QTableWidgetItem(QString::fromStdString(boundary.first_domain))); @@ -342,11 +342,28 @@ void GroomTool::update_shared_boundary_table() { table->setItem(row, 2, new QTableWidgetItem(QString::number(boundary.tolerance))); row++; } + + // Calculate the height to fit exactly the number of rows + int header_height = table->horizontalHeader()->height(); + int row_height = table->rowHeight(0); // Get height of first row (all rows should be same height) + int frame_width = table->frameWidth() * 2; // Top and bottom frame + int scrollbar_height = 0; + + // Check if horizontal scrollbar might be visible + if (table->horizontalScrollBar()->isVisible()) { + scrollbar_height = table->horizontalScrollBar()->height(); + } + + int total_height; if (table->rowCount() < 1) { - table->setMaximumHeight(ui_->label->height() * 2); + total_height = ui_->label->height() * 2; // Keep your original fallback } else { - table->setMaximumHeight(16777215); // default + total_height = header_height + (row_height * table->rowCount()) + frame_width + scrollbar_height; } + + table->setMaximumHeight(total_height); + table->setMinimumHeight(total_height); // Also set minimum to prevent shrinking + table->resizeColumnsToContents(); table->horizontalHeader()->setStretchLastSection(false); table->setSelectionBehavior(QAbstractItemView::SelectRows); @@ -508,7 +525,7 @@ void GroomTool::set_ui_from_params(GroomParameters params) { ui_->shared_boundary->setChecked(params.get_shared_boundaries_enabled()); ui_->shared_boundary_first_domain->setCurrentText(""); ui_->shared_boundary_second_domain->setCurrentText(""); - ui_->shared_boundary_tolerance->setValue(1e-1); + ui_->shared_boundary_tolerance->setText("0.01"); auto subjects = session_->get_project()->get_subjects(); int domain_id = std::max(ui_->domain_box->currentIndex(), 0); diff --git a/Studio/Groom/GroomTool.ui b/Studio/Groom/GroomTool.ui index 69eba12305..d3e79ecb53 100644 --- a/Studio/Groom/GroomTool.ui +++ b/Studio/Groom/GroomTool.ui @@ -355,7 +355,11 @@ QWidget#domain_panel { - + + + true + + @@ -381,28 +385,6 @@ QWidget#domain_panel { - - - - - 0 - 24 - - - - Qt::AlignCenter - - - 3 - - - 99999.000000000000000 - - - 0.100000000000000 - - - @@ -417,6 +399,13 @@ QWidget#domain_panel { + + + + Qt::AlignCenter + + + From d00bb26e98b50812236cb6fbf5bcb93a60fa2769 Mon Sep 17 00:00:00 2001 From: Alan Morris Date: Tue, 26 Aug 2025 12:59:59 -0600 Subject: [PATCH 07/18] Fix Viewer members when number of domains decreases (e.g. remove shared boundary) --- Studio/Visualization/Viewer.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Studio/Visualization/Viewer.cpp b/Studio/Visualization/Viewer.cpp index fc3de04bb4..4e2dec8f8f 100644 --- a/Studio/Visualization/Viewer.cpp +++ b/Studio/Visualization/Viewer.cpp @@ -1382,7 +1382,7 @@ std::shared_ptr Viewer::get_shape() { return shape_; } //----------------------------------------------------------------------------- void Viewer::initialize_surfaces() { - if (number_of_domains_ > surface_mappers_.size()) { + if (number_of_domains_ != surface_mappers_.size()) { surface_mappers_.resize(number_of_domains_); surface_actors_.resize(number_of_domains_); unclipped_surface_mappers_.resize(number_of_domains_); From 84ffae534cd66c20a4ae5d997f7ee73d6ebbc431 Mon Sep 17 00:00:00 2001 From: Alan Morris Date: Tue, 26 Aug 2025 13:00:14 -0600 Subject: [PATCH 08/18] Trigger save at end of groom so that project is consistent. --- Studio/Groom/GroomTool.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/Studio/Groom/GroomTool.cpp b/Studio/Groom/GroomTool.cpp index 98ec69854b..68ea36c6ea 100644 --- a/Studio/Groom/GroomTool.cpp +++ b/Studio/Groom/GroomTool.cpp @@ -730,6 +730,8 @@ void GroomTool::on_run_groom_button_clicked() { void GroomTool::handle_thread_complete() { Q_EMIT progress(95); + session_->trigger_save(); + QString duration = QString::number(timer_.elapsed() / 1000.0, 'f', 1); SW_LOG("Groom Complete. Duration: {} seconds", duration); From 286ad98eb900c8e793cd03dee8ad0f09f50d8ee0 Mon Sep 17 00:00:00 2001 From: Alan Morris Date: Thu, 28 Aug 2025 23:44:44 -0600 Subject: [PATCH 09/18] Fix domain type check in groom --- Libs/Groom/Groom.cpp | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/Libs/Groom/Groom.cpp b/Libs/Groom/Groom.cpp index e2f05502ad..c4a9de5b7a 100644 --- a/Libs/Groom/Groom.cpp +++ b/Libs/Groom/Groom.cpp @@ -1125,16 +1125,20 @@ Mesh Groom::get_mesh(int subject, int domain, bool transformed, MeshSource sourc std::string path; + DomainType domain_type; + if (source == MeshSource::Original) { if (domain >= subjects[subject]->get_original_filenames().size()) { throw std::out_of_range("domain index out of range"); } path = subjects[subject]->get_original_filenames()[domain]; + domain_type = project_->get_original_domain_types()[domain]; } else { if (domain >= subjects[subject]->get_groomed_filenames().size()) { throw std::out_of_range("domain index out of range"); } path = subjects[subject]->get_groomed_filenames()[domain]; + domain_type = project_->get_groomed_domain_types()[domain]; } auto constraint_filename = subjects[subject]->get_constraints_filenames(); @@ -1143,7 +1147,7 @@ Mesh Groom::get_mesh(int subject, int domain, bool transformed, MeshSource sourc constraint.read(constraint_filename[domain]); } - if (project_->get_original_domain_types()[domain] == DomainType::Contour) { + if (domain_type == DomainType::Contour) { Mesh mesh = MeshUtils::threadSafeReadMesh(path); mesh.set_id(subject); return mesh; @@ -1151,14 +1155,14 @@ Mesh Groom::get_mesh(int subject, int domain, bool transformed, MeshSource sourc Mesh mesh = vtkSmartPointer::New(); - if (project_->get_original_domain_types()[domain] == DomainType::Image) { + if (domain_type == DomainType::Image) { Image image(path); mesh = image.toMesh(0.5); constraint.clipMesh(mesh); - } else if (project_->get_original_domain_types()[domain] == DomainType::Mesh) { + } else if (domain_type == DomainType::Mesh) { mesh = MeshUtils::threadSafeReadMesh(path); constraint.clipMesh(mesh); - } else if (project_->get_original_domain_types()[domain] == DomainType::Contour) { + } else if (domain_type == DomainType::Contour) { mesh = MeshUtils::threadSafeReadMesh(path); } else { throw std::invalid_argument("invalid domain type"); From f272de8317e974d1ebac32270e44261d53b1500a Mon Sep 17 00:00:00 2001 From: Alan Morris Date: Thu, 28 Aug 2025 23:56:07 -0600 Subject: [PATCH 10/18] debugging logs --- Libs/Project/ProjectUtils.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/Libs/Project/ProjectUtils.cpp b/Libs/Project/ProjectUtils.cpp index 5b4d4ed97c..022ef6fc8d 100644 --- a/Libs/Project/ProjectUtils.cpp +++ b/Libs/Project/ProjectUtils.cpp @@ -294,10 +294,13 @@ static void assign_keys(StringMap& j, std::vector prefixes, std::ve auto prefix = prefixes[0]; if (filenames.size() != domains.size()) { SW_DEBUG("filename size: " + std::to_string(filenames.size())); + SW_DEBUG("domains size: " + std::to_string(domains.size())); for (auto& filename : filenames) { SW_DEBUG("filename: " + filename); } - + for (auto& domain : domains) { + SW_DEBUG("domain: " + domain); + } throw std::invalid_argument(prefix + " filenames and number of domains mismatch (" + std::to_string(filenames.size()) + " vs " + std::to_string(domains.size()) + ")"); } From a4658a1b269df106ec20a0577f93094fe3bf52dd Mon Sep 17 00:00:00 2001 From: Alan Morris Date: Thu, 28 Aug 2025 23:56:17 -0600 Subject: [PATCH 11/18] Adding biharmonic debugging --- Libs/Mesh/MeshWarper.cpp | 143 +++++++++++++++++++++++++++++++++++---- Libs/Mesh/MeshWarper.h | 5 +- 2 files changed, 132 insertions(+), 16 deletions(-) diff --git a/Libs/Mesh/MeshWarper.cpp b/Libs/Mesh/MeshWarper.cpp index f6eed42654..70991bd013 100644 --- a/Libs/Mesh/MeshWarper.cpp +++ b/Libs/Mesh/MeshWarper.cpp @@ -2,6 +2,9 @@ #include #include #include +#include +#include +#include #include #include #include @@ -23,11 +26,11 @@ static std::mutex mutex; //--------------------------------------------------------------------------- vtkSmartPointer MeshWarper::build_mesh(const Eigen::MatrixXd& particles) { - if (!this->warp_available_) { + if (!warp_available_) { return nullptr; } - if (!this->check_warp_ready()) { + if (!check_warp_ready()) { return nullptr; } @@ -38,7 +41,7 @@ vtkSmartPointer MeshWarper::build_mesh(const Eigen::MatrixXd& parti return blank; } - auto points = this->remove_bad_particles(particles); + auto points = remove_bad_particles(particles); vtkSmartPointer poly_data = MeshWarper::warp_mesh(points); @@ -59,9 +62,9 @@ void MeshWarper::set_reference_mesh(vtkSmartPointer reference_mesh, // lock so that we don't swap out the reference mesh while we are using it std::scoped_lock lock(mutex); - if (this->incoming_reference_mesh_ == reference_mesh) { - if (this->reference_particles_.size() == reference_particles.size()) { - if (this->reference_particles_ == reference_particles && landmarks_points_ == landmarks) { + if (incoming_reference_mesh_ == reference_mesh) { + if (reference_particles_.size() == reference_particles.size()) { + if (reference_particles_ == reference_particles && landmarks_points_ == landmarks) { // we can skip, nothing has changed return; } @@ -72,18 +75,18 @@ void MeshWarper::set_reference_mesh(vtkSmartPointer reference_mesh, incoming_reference_mesh_ = vtkSmartPointer::New(); incoming_reference_mesh_->DeepCopy(reference_mesh); - this->landmarks_points_ = landmarks; - this->reference_particles_ = reference_particles; + landmarks_points_ = landmarks; + reference_particles_ = reference_particles; // mark that the warp needs to be generated - this->needs_warp_ = true; + needs_warp_ = true; - this->warp_available_ = true; + warp_available_ = true; // TODO This is temporary for detecting if contour until the contour type is fully supported if (reference_mesh->GetNumberOfCells() > 0 && reference_mesh->GetCell(0)->GetNumberOfPoints() == 2) { - this->is_contour_ = true; - this->update_progress(1.0); + is_contour_ = true; + update_progress(1.0); } } @@ -410,7 +413,7 @@ vtkSmartPointer MeshWarper::remove_zero_area_triangles(vtkSmartPoin filtered_mesh->SetPoints(points); // Create a new cell array to store triangles - auto newTriangles = vtkSmartPointer::New(); + auto new_triangles = vtkSmartPointer::New(); // Get the triangle cells from the input mesh vtkCellArray* triangles = mesh->GetPolys(); @@ -436,13 +439,13 @@ vtkSmartPointer MeshWarper::remove_zero_area_triangles(vtkSmartPoin const double EPSILON = 1e-10; if (area > EPSILON) { // Add this triangle to the new cell array - newTriangles->InsertNextCell(cell_point_ids); + new_triangles->InsertNextCell(cell_point_ids); } } } // Set the triangles in the filtered mesh - filtered_mesh->SetPolys(newTriangles); + filtered_mesh->SetPolys(new_triangles); // Copy the point data filtered_mesh->GetPointData()->ShallowCopy(mesh->GetPointData()); @@ -558,6 +561,7 @@ bool MeshWarper::generate_warp_matrix(Eigen::MatrixXd TV, Eigen::MatrixXi TF, co const int k = 2; if (!igl::biharmonic_coordinates(TV, TF, S, k, W)) { SW_ERROR("Mesh Warp Error: igl:biharmonic_coordinates failed"); + diagnose_biharmonic_failure(TV, TF, S, k); return false; } // Throw away interior tet-vertices, keep weights and indices of boundary @@ -608,6 +612,115 @@ vtkSmartPointer MeshWarper::warp_mesh(const Eigen::MatrixXd& points return poly_data; } +//--------------------------------------------------------------------------- +void MeshWarper::diagnose_biharmonic_failure(const Eigen::MatrixXd& TV, const Eigen::MatrixXi& TF, + const std::vector>& S, int k) { + std::cerr << "\n=== Biharmonic Coordinates Failure Diagnostics ===\n"; + + // 1. Check matrix dimensions and consistency + std::cerr << "Vertices: " << TV.rows() << "x" << TV.cols() << std::endl; + std::cerr << "Faces: " << TF.rows() << "x" << TF.cols() << std::endl; + std::cerr << "Handles: " << S.size() << std::endl; + std::cerr << "k parameter: " << k << std::endl; + + // 2. Check handle indices validity + for (size_t i = 0; i < S.size(); ++i) { + for (int handle_idx : S[i]) { + if (handle_idx < 0 || handle_idx >= TV.rows()) { + std::cerr << "ERROR: Handle " << i << " contains invalid vertex index: " << handle_idx << " (valid range: 0-" + << (TV.rows() - 1) << ")\n"; + } + } + } + + // 3. Check face indices validity + int max_face_idx = TF.maxCoeff(); + int min_face_idx = TF.minCoeff(); + if (min_face_idx < 0 || max_face_idx >= TV.rows()) { + std::cerr << "ERROR: Face indices out of range. Min: " << min_face_idx << ", Max: " << max_face_idx + << " (valid range: 0-" << (TV.rows() - 1) << ")\n"; + } + + // 4. Check mesh connectivity + Eigen::VectorXi components; + int num_components = igl::facet_components(TF, components); + std::cerr << "Connected components: " << num_components << std::endl; + if (num_components > 1) { + std::cerr << "WARNING: Mesh has multiple disconnected components\n"; + } + + // 5. Check for boundary vertices (important for biharmonic coords) + Eigen::VectorXi boundary; + igl::boundary_loop(TF, boundary); + std::cerr << "Boundary vertices: " << boundary.rows() << std::endl; + + // 6. Check mesh genus/Euler characteristic + int V = TV.rows(); + int F = TF.rows(); + Eigen::MatrixXi E; + igl::edges(TF, E); + int euler_char = V - E.rows() + F; + std::cerr << "Euler characteristic: " << euler_char << " (genus ≈ " << (2 - euler_char) / 2 << ")\n"; + + // 7. Check for numerical issues + double min_coord = TV.minCoeff(); + double max_coord = TV.maxCoeff(); + double coord_range = max_coord - min_coord; + std::cerr << "Coordinate range: [" << min_coord << ", " << max_coord << "] (span: " << coord_range << ")\n"; + + if (coord_range < 1e-10) { + std::cerr << "ERROR: Coordinates too small - numerical precision issues likely\n"; + } + if (coord_range > 1e6) { + std::cerr << "WARNING: Large coordinate range - consider normalizing mesh\n"; + } + + // 8. Check triangle quality + Eigen::VectorXd areas; + igl::doublearea(TV, TF, areas); + areas *= 0.5; + + double min_area = areas.minCoeff(); + double max_area = areas.maxCoeff(); + std::cerr << "Triangle areas - Min: " << min_area << ", Max: " << max_area; + if (min_area > 0) { + std::cerr << ", Ratio: " << (max_area / min_area); + } + std::cerr << std::endl; + + if (min_area < 1e-15) { + std::cerr << "ERROR: Extremely small triangles detected - likely degenerate\n"; + } + + // 9. Test if we can build the Laplacian matrix (common failure point) + Eigen::SparseMatrix L; + try { + igl::cotmatrix(TV, TF, L); + std::cerr << "Cotangent Laplacian: " << L.rows() << "x" << L.cols() << " (" << L.nonZeros() << " non-zeros)\n"; + + // Check for NaN or inf in Laplacian + bool has_invalid = false; + for (int k = 0; k < L.outerSize(); ++k) { + for (Eigen::SparseMatrix::InnerIterator it(L, k); it; ++it) { + if (!std::isfinite(it.value())) { + has_invalid = true; + break; + } + } + if (has_invalid) break; + } + + if (has_invalid) { + std::cerr << "ERROR: Laplacian matrix contains invalid values (NaN/inf)\n"; + } + + } catch (const std::exception& e) { + std::cerr << "ERROR: Failed to build Laplacian matrix: " << e.what() << std::endl; + } + + std::cerr << "=== End Diagnostics ===\n\n"; +} + //--------------------------------------------------------------------------- Eigen::MatrixXd MeshWarper::extract_landmarks(vtkSmartPointer warped_mesh) { Eigen::MatrixXd landmarks; diff --git a/Libs/Mesh/MeshWarper.h b/Libs/Mesh/MeshWarper.h index 7ece8735c4..6b8db3aa01 100644 --- a/Libs/Mesh/MeshWarper.h +++ b/Libs/Mesh/MeshWarper.h @@ -111,6 +111,8 @@ class MeshWarper { //! Return the number of bad particles size_t bad_particle_count() const { return size_t(reference_particles_.rows()) - good_particles_.size(); } + void diagnose_biharmonic_failure(const Eigen::MatrixXd& TV, const Eigen::MatrixXi& TF, + const std::vector>& S, int k); // Members Eigen::MatrixXi faces_; Eigen::MatrixXd vertices_; @@ -123,7 +125,8 @@ class MeshWarper { bool warp_available_ = false; - std::map landmarks_map_; // map the landmarks id (Key) to the vertex(point) id (Value) belonging to the clean Reference mesh + std::map landmarks_map_; // map the landmarks id (Key) to the vertex(point) id (Value) belonging to the + // clean Reference mesh //! Reference mesh as it was given to us vtkSmartPointer incoming_reference_mesh_; //! Processed reference mesh From 9c14bf1bd7cc2c9f534c14794e85c2ec5d5c6622 Mon Sep 17 00:00:00 2001 From: Alan Morris Date: Fri, 29 Aug 2025 00:03:54 -0600 Subject: [PATCH 12/18] Rename vars for readability --- Libs/Mesh/MeshWarper.cpp | 49 +++++++++++++++++++--------------------- Libs/Mesh/MeshWarper.h | 3 ++- 2 files changed, 25 insertions(+), 27 deletions(-) diff --git a/Libs/Mesh/MeshWarper.cpp b/Libs/Mesh/MeshWarper.cpp index 70991bd013..398e615bbf 100644 --- a/Libs/Mesh/MeshWarper.cpp +++ b/Libs/Mesh/MeshWarper.cpp @@ -537,40 +537,37 @@ bool MeshWarper::generate_warp() { } //--------------------------------------------------------------------------- -bool MeshWarper::generate_warp_matrix(Eigen::MatrixXd TV, Eigen::MatrixXi TF, const Eigen::MatrixXd& Vref, - Eigen::MatrixXd& W) { - Eigen::VectorXi b; - { - Eigen::VectorXi J = Eigen::VectorXi::LinSpaced(TV.rows(), 0, TV.rows() - 1); - Eigen::VectorXd sqrD; - Eigen::MatrixXd _2; - // using J which is N by 1 instead of a matrix that represents faces of N by 3 - // so that we will find the closest vertices instead of closest point on the face - // so far the two meshes are not separated. So what we are really doing here - // is computing handles from low resolution and use that for the high resolution one - igl::point_mesh_squared_distance(Vref, TV, J, sqrD, b, _2); - // assert(sqrD.maxCoeff() < 1e-7 && "Particles must exist on vertices"); - } - - // list of points --> list of singleton lists - std::vector> S; - igl::matrix_to_list(b, S); +bool MeshWarper::generate_warp_matrix(Eigen::MatrixXd target_vertices, Eigen::MatrixXi target_faces, + const Eigen::MatrixXd& references_vertices, Eigen::MatrixXd& warp) { + Eigen::VectorXi closest_vertex_indices; + + Eigen::VectorXi vertex_indices = Eigen::VectorXi::LinSpaced(target_vertices.rows(), 0, target_vertices.rows() - 1); + Eigen::VectorXd squared_distances; + Eigen::MatrixXd unused_closest_points; + igl::point_mesh_squared_distance(references_vertices, target_vertices, vertex_indices, squared_distances, + closest_vertex_indices, unused_closest_points); + + std::vector> handle_sets; + igl::matrix_to_list(closest_vertex_indices, handle_sets); // Technically k should equal 3 for smooth interpolation in 3d, but 2 is // faster and looks OK const int k = 2; - if (!igl::biharmonic_coordinates(TV, TF, S, k, W)) { + if (!igl::biharmonic_coordinates(target_vertices, target_faces, handle_sets, k, warp)) { SW_ERROR("Mesh Warp Error: igl:biharmonic_coordinates failed"); - diagnose_biharmonic_failure(TV, TF, S, k); + diagnose_biharmonic_failure(target_vertices, target_faces, handle_sets, k); return false; } + // Throw away interior tet-vertices, keep weights and indices of boundary - Eigen::VectorXi I, J; - igl::remove_unreferenced(TV.rows(), TF, I, J); - std::for_each(TF.data(), TF.data() + TF.size(), [&I](int& a) { a = I(a); }); - std::for_each(b.data(), b.data() + b.size(), [&I](int& a) { a = I(a); }); - igl::slice(Eigen::MatrixXd(TV), J, 1, TV); - igl::slice(Eigen::MatrixXd(W), J, 1, W); + Eigen::VectorXi old_to_new_indices, referenced_vertices; + igl::remove_unreferenced(target_vertices.rows(), target_faces, old_to_new_indices, referenced_vertices); + std::for_each(target_faces.data(), target_faces.data() + target_faces.size(), + [&old_to_new_indices](int& a) { a = old_to_new_indices(a); }); + std::for_each(closest_vertex_indices.data(), closest_vertex_indices.data() + closest_vertex_indices.size(), + [&old_to_new_indices](int& a) { a = old_to_new_indices(a); }); + igl::slice(Eigen::MatrixXd(target_vertices), referenced_vertices, 1, target_vertices); + igl::slice(Eigen::MatrixXd(warp), referenced_vertices, 1, warp); return true; } diff --git a/Libs/Mesh/MeshWarper.h b/Libs/Mesh/MeshWarper.h index 6b8db3aa01..6630ce6dbb 100644 --- a/Libs/Mesh/MeshWarper.h +++ b/Libs/Mesh/MeshWarper.h @@ -103,7 +103,8 @@ class MeshWarper { static vtkSmartPointer recreate_mesh(vtkSmartPointer mesh); //! Generate the warp matrix - bool generate_warp_matrix(Eigen::MatrixXd TV, Eigen::MatrixXi TF, const Eigen::MatrixXd& Vref, Eigen::MatrixXd& W); + bool generate_warp_matrix(Eigen::MatrixXd target_vertices, Eigen::MatrixXi target_faces, + const Eigen::MatrixXd& references_vertices, Eigen::MatrixXd& warp); //! Generate a polydata from a set of points (e.g. warp the reference mesh) vtkSmartPointer warp_mesh(const Eigen::MatrixXd& points); From 863a511a58d2952e6e06f5cec7ccefe599533e21 Mon Sep 17 00:00:00 2001 From: Alan Morris Date: Fri, 29 Aug 2025 00:05:01 -0600 Subject: [PATCH 13/18] Remove "this->" from MeshWarper.cpp --- Libs/Mesh/MeshWarper.cpp | 112 +++++++++++++++++++-------------------- 1 file changed, 56 insertions(+), 56 deletions(-) diff --git a/Libs/Mesh/MeshWarper.cpp b/Libs/Mesh/MeshWarper.cpp index 398e615bbf..792d4a5756 100644 --- a/Libs/Mesh/MeshWarper.cpp +++ b/Libs/Mesh/MeshWarper.cpp @@ -48,7 +48,7 @@ vtkSmartPointer MeshWarper::build_mesh(const Eigen::MatrixXd& parti for (int i = 0; i < poly_data->GetNumberOfPoints(); i++) { double* p = poly_data->GetPoint(i); if (std::isnan(p[0]) || std::isnan(p[1]) || std::isnan(p[2])) { - this->warp_available_ = false; // failed + warp_available_ = false; // failed SW_ERROR("Reconstruction failed. NaN detected in mesh."); return nullptr; } @@ -91,18 +91,18 @@ void MeshWarper::set_reference_mesh(vtkSmartPointer reference_mesh, } //--------------------------------------------------------------------------- -bool MeshWarper::get_warp_available() { return this->warp_available_; } +bool MeshWarper::get_warp_available() { return warp_available_; } //--------------------------------------------------------------------------- bool MeshWarper::check_warp_ready() { std::scoped_lock lock(mutex); - if (!this->needs_warp_) { + if (!needs_warp_) { // warp already done return true; } - return this->generate_warp(); + return generate_warp(); } //--------------------------------------------------------------------------- @@ -122,15 +122,15 @@ void MeshWarper::add_particle_vertices(Eigen::MatrixXd& vertices) { while (not_all_done) { not_all_done = false; - this->reference_mesh_->BuildLinks(); + reference_mesh_->BuildLinks(); auto locator = vtkSmartPointer::New(); locator->SetCacheCellBounds(true); - locator->SetDataSet(this->reference_mesh_); + locator->SetDataSet(reference_mesh_); locator->BuildLocator(); std::vector> new_triangles; vtkSmartPointer new_points = vtkSmartPointer::New(); - int next_point = this->reference_mesh_->GetNumberOfPoints(); + int next_point = reference_mesh_->GetNumberOfPoints(); for (int i = 0; i < vertices.rows(); i++) { if (done[i]) { @@ -149,7 +149,7 @@ void MeshWarper::add_particle_vertices(Eigen::MatrixXd& vertices) { pt[2] = closest_point[2]; // grab the closest cell - vtkCell* cell = this->reference_mesh_->GetCell(cell_id); + vtkCell* cell = reference_mesh_->GetCell(cell_id); if (cell->GetCellType() == VTK_EMPTY_CELL) { // this one was deleted not_all_done = true; @@ -159,14 +159,14 @@ void MeshWarper::add_particle_vertices(Eigen::MatrixXd& vertices) { // mark this particle as complete done[i] = true; - this->update_progress(count++ / vertices.rows()); + update_progress(count++ / vertices.rows()); double point[3] = {pt[0], pt[1], pt[2]}; double closest[3]; double pcoords[3]; double dist2; double weights[3]; - this->reference_mesh_->GetCell(cell_id)->EvaluatePosition(point, closest, sub_id, pcoords, dist2, weights); + reference_mesh_->GetCell(cell_id)->EvaluatePosition(point, closest, sub_id, pcoords, dist2, weights); if (weights[0] > 0.99 || weights[1] > 0.99 || weights[2] > 0.99) { // close enough to a vertex @@ -195,11 +195,11 @@ void MeshWarper::add_particle_vertices(Eigen::MatrixXd& vertices) { } auto neighbors = vtkSmartPointer::New(); - this->reference_mesh_->GetCellEdgeNeighbors(cell_id, v0_index, v1_index, neighbors); + reference_mesh_->GetCellEdgeNeighbors(cell_id, v0_index, v1_index, neighbors); if (on_edge) { // first confirm that the other cell hasn't been deleted already if (neighbors->GetNumberOfIds() == 1) { // could be 0 for the boundary of an open mesh - vtkCell* other_cell = this->reference_mesh_->GetCell(neighbors->GetId(0)); + vtkCell* other_cell = reference_mesh_->GetCell(neighbors->GetId(0)); if (other_cell->GetCellType() == VTK_EMPTY_CELL) { // this one was deleted not_all_done = true; continue; @@ -214,11 +214,11 @@ void MeshWarper::add_particle_vertices(Eigen::MatrixXd& vertices) { if (on_edge) { if (neighbors->GetNumberOfIds() == 1) { // could be 0 for the boundary of an open mesh // split the neighbor cell into two triangles as well - this->split_cell_on_edge(neighbors->GetId(0), new_vertex, v0_index, v1_index, new_triangles); + split_cell_on_edge(neighbors->GetId(0), new_vertex, v0_index, v1_index, new_triangles); } // split the current cell into two triangles - this->split_cell_on_edge(cell_id, new_vertex, v0_index, v1_index, new_triangles); + split_cell_on_edge(cell_id, new_vertex, v0_index, v1_index, new_triangles); } else { // interior to a current triangle vtkSmartPointer list = vtkSmartPointer::New(); list->SetNumberOfIds(3); @@ -242,22 +242,22 @@ void MeshWarper::add_particle_vertices(Eigen::MatrixXd& vertices) { list->SetId(2, cell->GetPointId(2)); new_triangles.push_back(list); - this->reference_mesh_->DeleteCell(cell_id); + reference_mesh_->DeleteCell(cell_id); } } // add new points for (int i = 0; i < new_points->GetNumberOfPoints(); i++) { - this->reference_mesh_->GetPoints()->InsertNextPoint(new_points->GetPoint(i)); + reference_mesh_->GetPoints()->InsertNextPoint(new_points->GetPoint(i)); } // add new triangles for (int i = 0; i < new_triangles.size(); i++) { - this->reference_mesh_->InsertNextCell(VTK_TRIANGLE, new_triangles[i]); + reference_mesh_->InsertNextCell(VTK_TRIANGLE, new_triangles[i]); } // remove deleted triangles and clean up - this->reference_mesh_ = MeshWarper::recreate_mesh(this->reference_mesh_); + reference_mesh_ = MeshWarper::recreate_mesh(reference_mesh_); } } @@ -265,19 +265,19 @@ void MeshWarper::add_particle_vertices(Eigen::MatrixXd& vertices) { void MeshWarper::find_good_particles() { // Create the tree auto tree = vtkSmartPointer::New(); - tree->SetDataSet(this->reference_mesh_); + tree->SetDataSet(reference_mesh_); tree->BuildLocator(); std::vector ids; - for (int i = 0; i < this->vertices_.rows(); i++) { - double p[3]{this->vertices_(i, 0), this->vertices_(i, 1), this->vertices_(i, 2)}; + for (int i = 0; i < vertices_.rows(); i++) { + double p[3]{vertices_(i, 0), vertices_(i, 1), vertices_(i, 2)}; int id = tree->FindClosestPoint(p); ids.push_back(id); } std::set set; // initially store in set to avoid duplicates - for (int i = 0; i < this->vertices_.rows(); i++) { - for (int j = i + 1; j < this->vertices_.rows(); j++) { + for (int i = 0; i < vertices_.rows(); i++) { + for (int j = i + 1; j < vertices_.rows(); j++) { if (ids[i] == ids[j]) { set.insert(i); set.insert(j); @@ -285,19 +285,19 @@ void MeshWarper::find_good_particles() { } } - this->good_particles_.clear(); - for (int i = 0; i < this->vertices_.rows(); i++) { + good_particles_.clear(); + for (int i = 0; i < vertices_.rows(); i++) { if (set.find(i) == set.end()) { - this->good_particles_.push_back(i); + good_particles_.push_back(i); } } } //--------------------------------------------------------------------------- Eigen::MatrixXd MeshWarper::remove_bad_particles(const Eigen::MatrixXd& particles) { - Eigen::MatrixXd new_particles(this->good_particles_.size(), 3); - for (int i = 0; i < this->good_particles_.size(); i++) { - int id = this->good_particles_[i]; + Eigen::MatrixXd new_particles(good_particles_.size(), 3); + for (int i = 0; i < good_particles_.size(); i++) { + int id = good_particles_[i]; new_particles(i, 0) = particles(id, 0); new_particles(i, 1) = particles(id, 1); new_particles(i, 2) = particles(id, 2); @@ -309,7 +309,7 @@ Eigen::MatrixXd MeshWarper::remove_bad_particles(const Eigen::MatrixXd& particle //--------------------------------------------------------------------------- void MeshWarper::split_cell_on_edge(int cell_id, int new_vertex, int v0, int v1, std::vector>& new_triangles) { - vtkCell* cell = this->reference_mesh_->GetCell(cell_id); + vtkCell* cell = reference_mesh_->GetCell(cell_id); int p0 = cell->GetPointId(0); int p1 = cell->GetPointId(1); int p2 = cell->GetPointId(2); @@ -366,7 +366,7 @@ void MeshWarper::split_cell_on_edge(int cell_id, int new_vertex, int v0, int v1, new_triangles.push_back(list); } - this->reference_mesh_->DeleteCell(cell_id); + reference_mesh_->DeleteCell(cell_id); } //--------------------------------------------------------------------------- @@ -494,14 +494,14 @@ vtkSmartPointer MeshWarper::recreate_mesh(vtkSmartPointerupdate_progress(1.0); - this->needs_warp_ = false; + update_progress(1.0); + needs_warp_ = false; return true; } - this->reference_mesh_ = MeshWarper::prep_mesh(this->incoming_reference_mesh_); + reference_mesh_ = MeshWarper::prep_mesh(incoming_reference_mesh_); // prep points - this->vertices_ = this->reference_particles_; + vertices_ = reference_particles_; if (reference_mesh_->GetNumberOfPoints() == 0) { SW_ERROR("Unable to warp mesh, no points in surface mesh"); @@ -510,29 +510,29 @@ bool MeshWarper::generate_warp() { return false; } - this->add_particle_vertices(this->vertices_); + add_particle_vertices(vertices_); - if (this->landmarks_points_.size() > 0) { + if (landmarks_points_.size() > 0) { // to ensure that the landmark points sit on the ref mesh vertices - this->add_particle_vertices(this->landmarks_points_); - this->find_landmarks_vertices_on_ref_mesh(); + add_particle_vertices(landmarks_points_); + find_landmarks_vertices_on_ref_mesh(); } - this->find_good_particles(); - this->vertices_ = this->remove_bad_particles(this->vertices_); + find_good_particles(); + vertices_ = remove_bad_particles(vertices_); Mesh referenceMesh(reference_mesh_); Eigen::MatrixXd vertices = referenceMesh.points(); - this->faces_ = referenceMesh.faces(); + faces_ = referenceMesh.faces(); // perform warp - if (!MeshWarper::generate_warp_matrix(vertices, this->faces_, this->vertices_, this->warp_)) { - this->update_progress(1.0); - this->warp_available_ = false; + if (!MeshWarper::generate_warp_matrix(vertices, faces_, vertices_, warp_)) { + update_progress(1.0); + warp_available_ = false; return false; } - this->update_progress(1.0); - this->needs_warp_ = false; + update_progress(1.0); + needs_warp_ = false; return true; } @@ -574,10 +574,10 @@ bool MeshWarper::generate_warp_matrix(Eigen::MatrixXd target_vertices, Eigen::Ma //--------------------------------------------------------------------------- bool MeshWarper::find_landmarks_vertices_on_ref_mesh() { auto tree = vtkSmartPointer::New(); - tree->SetDataSet(this->reference_mesh_); + tree->SetDataSet(reference_mesh_); tree->BuildLocator(); - for (int i = 0; i < this->landmarks_points_.rows(); i++) { - double p[3]{this->landmarks_points_(i, 0), this->landmarks_points_(i, 1), this->landmarks_points_(i, 2)}; + for (int i = 0; i < landmarks_points_.rows(); i++) { + double p[3]{landmarks_points_(i, 0), landmarks_points_(i, 1), landmarks_points_(i, 2)}; int id = tree->FindClosestPoint(p); landmarks_map_.insert({i, id}); // std::cout << "Landmark id: " << i << " Vertex id: " << id << std::endl; @@ -587,9 +587,9 @@ bool MeshWarper::find_landmarks_vertices_on_ref_mesh() { //--------------------------------------------------------------------------- vtkSmartPointer MeshWarper::warp_mesh(const Eigen::MatrixXd& points) { - int num_vertices = this->warp_.rows(); - int num_faces = this->faces_.rows(); - Eigen::MatrixXd v_out = this->warp_ * (points.rowwise() + Eigen::RowVector3d(0, 0, 0)); + int num_vertices = warp_.rows(); + int num_faces = faces_.rows(); + Eigen::MatrixXd v_out = warp_ * (points.rowwise() + Eigen::RowVector3d(0, 0, 0)); vtkSmartPointer poly_data = vtkSmartPointer::New(); vtkSmartPointer out_points = vtkSmartPointer::New(); out_points->SetNumberOfPoints(num_vertices); @@ -599,9 +599,9 @@ vtkSmartPointer MeshWarper::warp_mesh(const Eigen::MatrixXd& points vtkSmartPointer polys = vtkSmartPointer::New(); for (vtkIdType i = 0; i < num_faces; i++) { polys->InsertNextCell(3); - polys->InsertCellPoint(this->faces_(i, 0)); - polys->InsertCellPoint(this->faces_(i, 1)); - polys->InsertCellPoint(this->faces_(i, 2)); + polys->InsertCellPoint(faces_(i, 0)); + polys->InsertCellPoint(faces_(i, 1)); + polys->InsertCellPoint(faces_(i, 2)); } poly_data->SetPoints(out_points); poly_data->SetPolys(polys); From a8bc2bf4c137ca49fec0842d8eb48f8bc9c8505b Mon Sep 17 00:00:00 2001 From: Alan Morris Date: Fri, 29 Aug 2025 00:48:32 -0600 Subject: [PATCH 14/18] When igl::biharmonic fails, we retry again with remeshing. --- Libs/Mesh/MeshWarper.cpp | 137 +++++++++++++++++++++------------------ Libs/Mesh/MeshWarper.h | 2 + 2 files changed, 75 insertions(+), 64 deletions(-) diff --git a/Libs/Mesh/MeshWarper.cpp b/Libs/Mesh/MeshWarper.cpp index 792d4a5756..d4b7200ddc 100644 --- a/Libs/Mesh/MeshWarper.cpp +++ b/Libs/Mesh/MeshWarper.cpp @@ -371,12 +371,12 @@ void MeshWarper::split_cell_on_edge(int cell_id, int new_vertex, int v0, int v1, //--------------------------------------------------------------------------- vtkSmartPointer MeshWarper::prep_mesh(vtkSmartPointer mesh) { - vtkSmartPointer triangle_filter = vtkSmartPointer::New(); + auto triangle_filter = vtkSmartPointer::New(); triangle_filter->SetInputData(mesh); triangle_filter->PassLinesOff(); triangle_filter->Update(); - vtkSmartPointer connectivity = vtkSmartPointer::New(); + auto connectivity = vtkSmartPointer::New(); connectivity->SetInputConnection(triangle_filter->GetOutputPort()); connectivity->SetExtractionModeToLargestRegion(); connectivity->Update(); @@ -388,7 +388,9 @@ vtkSmartPointer MeshWarper::prep_mesh(vtkSmartPointer // remove deleted triangles and clean up auto recreated = MeshWarper::recreate_mesh(fixed.getVTKMesh()); - return MeshWarper::remove_zero_area_triangles(recreated); + auto final = MeshWarper::remove_zero_area_triangles(recreated); + + return final; } //--------------------------------------------------------------------------- @@ -498,42 +500,57 @@ bool MeshWarper::generate_warp() { needs_warp_ = false; return true; } - reference_mesh_ = MeshWarper::prep_mesh(incoming_reference_mesh_); - // prep points - vertices_ = reference_particles_; + const int MAX_ATTEMPTS = 2; - if (reference_mesh_->GetNumberOfPoints() == 0) { - SW_ERROR("Unable to warp mesh, no points in surface mesh"); - update_progress(1.0); - warp_available_ = false; - return false; - } + for (int attempt = 0; attempt < MAX_ATTEMPTS; ++attempt) { + // Prepare mesh (remesh on retry) + if (attempt == 0) { + reference_mesh_ = MeshWarper::prep_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()); + } + + // Prep points + vertices_ = reference_particles_; + if (reference_mesh_->GetNumberOfPoints() == 0) { + SW_ERROR("Unable to warp mesh, no points in surface mesh"); + update_progress(1.0); + warp_available_ = false; + return false; + } - add_particle_vertices(vertices_); + add_particle_vertices(vertices_); + if (landmarks_points_.size() > 0) { + // to ensure that the landmark points sit on the ref mesh vertices + add_particle_vertices(landmarks_points_); + find_landmarks_vertices_on_ref_mesh(); + } - if (landmarks_points_.size() > 0) { - // to ensure that the landmark points sit on the ref mesh vertices - add_particle_vertices(landmarks_points_); - find_landmarks_vertices_on_ref_mesh(); + find_good_particles(); + vertices_ = remove_bad_particles(vertices_); + Mesh reference_mesh(reference_mesh_); + Eigen::MatrixXd vertices = reference_mesh.points(); + faces_ = reference_mesh.faces(); + + // Perform warp + if (MeshWarper::generate_warp_matrix(vertices, faces_, vertices_, warp_)) { + // Success! + update_progress(1.0); + needs_warp_ = false; + return true; + } } - find_good_particles(); - vertices_ = remove_bad_particles(vertices_); + SW_ERROR("Mesh Warp Error: igl:biharmonic_coordinates failed"); - Mesh referenceMesh(reference_mesh_); - Eigen::MatrixXd vertices = referenceMesh.points(); - faces_ = referenceMesh.faces(); - - // perform warp - if (!MeshWarper::generate_warp_matrix(vertices, faces_, vertices_, warp_)) { - update_progress(1.0); - warp_available_ = false; - return false; - } + // All attempts failed update_progress(1.0); - needs_warp_ = false; - return true; + warp_available_ = false; + return false; } //--------------------------------------------------------------------------- @@ -554,7 +571,6 @@ bool MeshWarper::generate_warp_matrix(Eigen::MatrixXd target_vertices, Eigen::Ma // faster and looks OK const int k = 2; if (!igl::biharmonic_coordinates(target_vertices, target_faces, handle_sets, k, warp)) { - SW_ERROR("Mesh Warp Error: igl:biharmonic_coordinates failed"); diagnose_biharmonic_failure(target_vertices, target_faces, handle_sets, k); return false; } @@ -609,23 +625,23 @@ vtkSmartPointer MeshWarper::warp_mesh(const Eigen::MatrixXd& points return poly_data; } +//--------------------------------------------------------------------------- //--------------------------------------------------------------------------- void MeshWarper::diagnose_biharmonic_failure(const Eigen::MatrixXd& TV, const Eigen::MatrixXi& TF, const std::vector>& S, int k) { - std::cerr << "\n=== Biharmonic Coordinates Failure Diagnostics ===\n"; - + SW_LOG_ONLY("=== Biharmonic Coordinates Failure Diagnostics ==="); // 1. Check matrix dimensions and consistency - std::cerr << "Vertices: " << TV.rows() << "x" << TV.cols() << std::endl; - std::cerr << "Faces: " << TF.rows() << "x" << TF.cols() << std::endl; - std::cerr << "Handles: " << S.size() << std::endl; - std::cerr << "k parameter: " << k << std::endl; + SW_LOG_ONLY("Vertices: {}x{}", TV.rows(), TV.cols()); + SW_LOG_ONLY("Faces: {}x{}", TF.rows(), TF.cols()); + SW_LOG_ONLY("Handles: {}", S.size()); + SW_LOG_ONLY("k parameter: {}", k); // 2. Check handle indices validity for (size_t i = 0; i < S.size(); ++i) { for (int handle_idx : S[i]) { if (handle_idx < 0 || handle_idx >= TV.rows()) { - std::cerr << "ERROR: Handle " << i << " contains invalid vertex index: " << handle_idx << " (valid range: 0-" - << (TV.rows() - 1) << ")\n"; + SW_LOG_ONLY("ERROR: Handle {} contains invalid vertex index: {} (valid range: 0-{})", i, handle_idx, + TV.rows() - 1); } } } @@ -634,22 +650,22 @@ void MeshWarper::diagnose_biharmonic_failure(const Eigen::MatrixXd& TV, const Ei int max_face_idx = TF.maxCoeff(); int min_face_idx = TF.minCoeff(); if (min_face_idx < 0 || max_face_idx >= TV.rows()) { - std::cerr << "ERROR: Face indices out of range. Min: " << min_face_idx << ", Max: " << max_face_idx - << " (valid range: 0-" << (TV.rows() - 1) << ")\n"; + SW_LOG_ONLY("ERROR: Face indices out of range. Min: {}, Max: {} (valid range: 0-{})", min_face_idx, max_face_idx, + TV.rows() - 1); } // 4. Check mesh connectivity Eigen::VectorXi components; int num_components = igl::facet_components(TF, components); - std::cerr << "Connected components: " << num_components << std::endl; + SW_LOG_ONLY("Connected components: {}", num_components); if (num_components > 1) { - std::cerr << "WARNING: Mesh has multiple disconnected components\n"; + SW_LOG_ONLY("WARNING: Mesh has multiple disconnected components"); } // 5. Check for boundary vertices (important for biharmonic coords) Eigen::VectorXi boundary; igl::boundary_loop(TF, boundary); - std::cerr << "Boundary vertices: " << boundary.rows() << std::endl; + SW_LOG_ONLY("Boundary vertices: {}", boundary.rows()); // 6. Check mesh genus/Euler characteristic int V = TV.rows(); @@ -657,44 +673,41 @@ void MeshWarper::diagnose_biharmonic_failure(const Eigen::MatrixXd& TV, const Ei Eigen::MatrixXi E; igl::edges(TF, E); int euler_char = V - E.rows() + F; - std::cerr << "Euler characteristic: " << euler_char << " (genus ≈ " << (2 - euler_char) / 2 << ")\n"; + SW_LOG_ONLY("Euler characteristic: {} (genus ≈ {})", euler_char, (2 - euler_char) / 2); // 7. Check for numerical issues double min_coord = TV.minCoeff(); double max_coord = TV.maxCoeff(); double coord_range = max_coord - min_coord; - std::cerr << "Coordinate range: [" << min_coord << ", " << max_coord << "] (span: " << coord_range << ")\n"; - + SW_LOG_ONLY("Coordinate range: [{}, {}] (span: {})", min_coord, max_coord, coord_range); if (coord_range < 1e-10) { - std::cerr << "ERROR: Coordinates too small - numerical precision issues likely\n"; + SW_LOG_ONLY("ERROR: Coordinates too small - numerical precision issues likely"); } if (coord_range > 1e6) { - std::cerr << "WARNING: Large coordinate range - consider normalizing mesh\n"; + SW_LOG_ONLY("WARNING: Large coordinate range - consider normalizing mesh"); } // 8. Check triangle quality Eigen::VectorXd areas; igl::doublearea(TV, TF, areas); areas *= 0.5; - double min_area = areas.minCoeff(); double max_area = areas.maxCoeff(); - std::cerr << "Triangle areas - Min: " << min_area << ", Max: " << max_area; if (min_area > 0) { - std::cerr << ", Ratio: " << (max_area / min_area); + SW_LOG_ONLY("Triangle areas - Min: {}, Max: {}, Ratio: {}", min_area, max_area, max_area / min_area); + } else { + SW_LOG_ONLY("Triangle areas - Min: {}, Max: {}", min_area, max_area); } - std::cerr << std::endl; if (min_area < 1e-15) { - std::cerr << "ERROR: Extremely small triangles detected - likely degenerate\n"; + SW_LOG_ONLY("ERROR: Extremely small triangles detected - likely degenerate"); } // 9. Test if we can build the Laplacian matrix (common failure point) Eigen::SparseMatrix L; try { igl::cotmatrix(TV, TF, L); - std::cerr << "Cotangent Laplacian: " << L.rows() << "x" << L.cols() << " (" << L.nonZeros() << " non-zeros)\n"; - + SW_LOG_ONLY("Cotangent Laplacian: {}x{} ({} non-zeros)", L.rows(), L.cols(), L.nonZeros()); // Check for NaN or inf in Laplacian bool has_invalid = false; for (int k = 0; k < L.outerSize(); ++k) { @@ -706,18 +719,14 @@ void MeshWarper::diagnose_biharmonic_failure(const Eigen::MatrixXd& TV, const Ei } if (has_invalid) break; } - if (has_invalid) { - std::cerr << "ERROR: Laplacian matrix contains invalid values (NaN/inf)\n"; + SW_LOG_ONLY("ERROR: Laplacian matrix contains invalid values (NaN/inf)"); } - } catch (const std::exception& e) { - std::cerr << "ERROR: Failed to build Laplacian matrix: " << e.what() << std::endl; + SW_LOG_ONLY("ERROR: Failed to build Laplacian matrix: {}", e.what()); } - - std::cerr << "=== End Diagnostics ===\n\n"; + SW_LOG_ONLY("=== End Diagnostics ==="); } - //--------------------------------------------------------------------------- Eigen::MatrixXd MeshWarper::extract_landmarks(vtkSmartPointer warped_mesh) { Eigen::MatrixXd landmarks; diff --git a/Libs/Mesh/MeshWarper.h b/Libs/Mesh/MeshWarper.h index 6630ce6dbb..500e6deb71 100644 --- a/Libs/Mesh/MeshWarper.h +++ b/Libs/Mesh/MeshWarper.h @@ -114,6 +114,8 @@ class MeshWarper { void diagnose_biharmonic_failure(const Eigen::MatrixXd& TV, const Eigen::MatrixXi& TF, const std::vector>& S, int k); + + // Members Eigen::MatrixXi faces_; Eigen::MatrixXd vertices_; From 72d1edae204c69a659a15506be630b15429928a0 Mon Sep 17 00:00:00 2001 From: Alan Morris Date: Fri, 29 Aug 2025 00:52:42 -0600 Subject: [PATCH 15/18] Upgrade logging so that log-only messages become debug when debug is enabled --- Libs/Common/Logging.cpp | 8 +++++++- Libs/Common/Logging.h | 4 ++-- 2 files changed, 9 insertions(+), 3 deletions(-) diff --git a/Libs/Common/Logging.cpp b/Libs/Common/Logging.cpp index 7e527543c2..fdd92b975c 100644 --- a/Libs/Common/Logging.cpp +++ b/Libs/Common/Logging.cpp @@ -70,7 +70,13 @@ void Logging::log_message(const std::string& message, const int line, const char } //----------------------------------------------------------------------------- -void Logging::log_only(const std::string& message, const int line, const char* file) const { +void Logging::log_only(const std::string& message, const int line, const char* file, const char* function) const { + if (spd::get_level() == spd::level::debug) { + // when in debug mode, treat this as a debug message + log_debug(message, line, file, function); + return; + } + spd::info(message); if (log_open_) { spd::get("file")->info(message); diff --git a/Libs/Common/Logging.h b/Libs/Common/Logging.h index 1e3961097f..eaaeaf887f 100644 --- a/Libs/Common/Logging.h +++ b/Libs/Common/Logging.h @@ -103,7 +103,7 @@ class Logging { void log_message(const std::string& message, const int line, const char* file) const; //! Log a message, use SW_LOG_ONLY macro - void log_only(const std::string& message, const int line, const char* file) const; + void log_only(const std::string& message, const int line, const char* file, const char *function) const; //! Log a stack trace message, use SW_LOG_STACK macro void log_stack(const std::string& message) const; @@ -176,7 +176,7 @@ class Logging { //! Log only macro #define SW_LOG_ONLY(message, ...) \ - shapeworks::Logging::Instance().log_only(safe_format(message, ##__VA_ARGS__), __LINE__, __FILE__); + shapeworks::Logging::Instance().log_only(safe_format(message, ##__VA_ARGS__), __LINE__, __FILE__, __FUNCTION__); //! Log warning macro #define SW_WARN(message, ...) \ From 920d9fb728001b192d99cb5db35d8e0d5e8cb0f6 Mon Sep 17 00:00:00 2001 From: Alan Morris Date: Fri, 29 Aug 2025 10:12:14 -0600 Subject: [PATCH 16/18] Added some debugging, update DomainType enum --- Libs/Groom/Groom.cpp | 11 ++++++++--- Libs/Optimize/Domain/DomainType.h | 10 ++-------- 2 files changed, 10 insertions(+), 11 deletions(-) diff --git a/Libs/Groom/Groom.cpp b/Libs/Groom/Groom.cpp index c4a9de5b7a..8f2b74240c 100644 --- a/Libs/Groom/Groom.cpp +++ b/Libs/Groom/Groom.cpp @@ -1125,7 +1125,7 @@ Mesh Groom::get_mesh(int subject, int domain, bool transformed, MeshSource sourc std::string path; - DomainType domain_type; + DomainType domain_type = DomainType::Mesh; if (source == MeshSource::Original) { if (domain >= subjects[subject]->get_original_filenames().size()) { @@ -1133,12 +1133,17 @@ Mesh Groom::get_mesh(int subject, int domain, bool transformed, MeshSource sourc } path = subjects[subject]->get_original_filenames()[domain]; domain_type = project_->get_original_domain_types()[domain]; + SW_DEBUG("Getting original mesh for subject {}, domain {}: {}", subject, domain, path); + SW_DEBUG("Domain type: {}", static_cast(domain_type)); } else { if (domain >= subjects[subject]->get_groomed_filenames().size()) { throw std::out_of_range("domain index out of range"); } path = subjects[subject]->get_groomed_filenames()[domain]; - domain_type = project_->get_groomed_domain_types()[domain]; + + domain_type = ProjectUtils::determine_domain_type(path); + SW_DEBUG("Getting groomed mesh for subject {}, domain {}: {}", subject, domain, path); + SW_DEBUG("Domain type: {}", static_cast(domain_type)); } auto constraint_filename = subjects[subject]->get_constraints_filenames(); @@ -1165,7 +1170,7 @@ Mesh Groom::get_mesh(int subject, int domain, bool transformed, MeshSource sourc } else if (domain_type == DomainType::Contour) { mesh = MeshUtils::threadSafeReadMesh(path); } else { - throw std::invalid_argument("invalid domain type"); + throw std::invalid_argument("invalid domain type: " + std::to_string(static_cast(domain_type))); } if (transformed) { diff --git a/Libs/Optimize/Domain/DomainType.h b/Libs/Optimize/Domain/DomainType.h index 111e2908d9..c6b253401d 100644 --- a/Libs/Optimize/Domain/DomainType.h +++ b/Libs/Optimize/Domain/DomainType.h @@ -1,11 +1,5 @@ #pragma once -#define DIMENSION 3 - namespace shapeworks { - enum class DomainType : char { - Image = 'I', - Mesh = 'M', - Contour = 'C' - }; -} +enum class DomainType { Image, Mesh, Contour }; +} // namespace shapeworks From 5bb650358eb7e7e7f7b816fb19267fc3c258a4e7 Mon Sep 17 00:00:00 2001 From: Alan Morris Date: Fri, 29 Aug 2025 12:05:19 -0600 Subject: [PATCH 17/18] Fix compilation for gcc --- Libs/Groom/Groom.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/Libs/Groom/Groom.cpp b/Libs/Groom/Groom.cpp index 8f2b74240c..155ba085e2 100644 --- a/Libs/Groom/Groom.cpp +++ b/Libs/Groom/Groom.cpp @@ -804,8 +804,10 @@ bool Groom::run_shared_boundaries() { Mesh output_contour(empty_mesh); try { // returns left remainder, right remainder, and shared_surface - std::tie(extracted_l, extracted_r, extracted_s) = - MeshUtils::shared_boundary_extractor(first_mesh, second_mesh, shared_boundary.tolerance); + auto result = MeshUtils::shared_boundary_extractor(first_mesh, second_mesh, shared_boundary.tolerance); + extracted_l = result[0]; + extracted_r = result[1]; + extracted_s = result[2]; extracted_l.remeshPercent(.99, 1.0); extracted_r.remeshPercent(.99, 1.0); From 7b9b565374e3377ed13dce2ccb82b86cf5797630 Mon Sep 17 00:00:00 2001 From: Alan Morris Date: Sat, 30 Aug 2025 02:09:09 -0600 Subject: [PATCH 18/18] Only clean up groomed files for non-fixed domains. --- Libs/Groom/Groom.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/Libs/Groom/Groom.cpp b/Libs/Groom/Groom.cpp index 155ba085e2..005b40ec35 100644 --- a/Libs/Groom/Groom.cpp +++ b/Libs/Groom/Groom.cpp @@ -994,9 +994,12 @@ void Groom::clear_unused_shared_boundaries() { project_->set_original_domain_types(original_domain_types); project_->set_groomed_domain_types(groomed_domain_types); - // Clean up files from all subjects + // Clean up files from all non-fixed subjects auto subjects = project_->get_subjects(); for (auto& subject : subjects) { + if (subject->is_fixed()) { + continue; // Skip fixed subjects + } // Remove shared surface/boundary files using lambda auto original_filenames = subject->get_original_filenames(); remove_shared_files(original_filenames);