diff --git a/Applications/shapeworks/MeshCommands.cpp b/Applications/shapeworks/MeshCommands.cpp index 19a374ddf6..a6aee4b1a2 100644 --- a/Applications/shapeworks/MeshCommands.cpp +++ b/Applications/shapeworks/MeshCommands.cpp @@ -1505,7 +1505,7 @@ bool WarpMesh::execute(const optparse::Values &options, SharedCommandData &share Eigen::MatrixXd vertices_updated = output.points(); std::map landmarks_map = warper.get_landmarks_map(); - for (vtkIdType i = 0; i < warper.get_warp_matrix().rows(); i++) { + for (vtkIdType i = 0; i < vertices_updated.rows(); i++) { if (landmarks_map.count(i)) { // This vertex i corresponds to a landmark int landmark_index = landmarks_map.lower_bound(i)->second; diff --git a/Libs/Mesh/MeshWarper.cpp b/Libs/Mesh/MeshWarper.cpp index 0a2d18ba8d..cc066b71e9 100644 --- a/Libs/Mesh/MeshWarper.cpp +++ b/Libs/Mesh/MeshWarper.cpp @@ -2,10 +2,13 @@ #include #include #include +#include #include #include +#include #include #include +#include #include #include #include @@ -19,6 +22,7 @@ #include #include +#include namespace shapeworks { @@ -27,6 +31,7 @@ static std::mutex mutex; //--------------------------------------------------------------------------- vtkSmartPointer MeshWarper::build_mesh(const Eigen::MatrixXd& particles) { + TIME_SCOPE("MeshWarper::build_mesh"); if (!warp_available_) { return nullptr; } @@ -44,7 +49,8 @@ vtkSmartPointer MeshWarper::build_mesh(const Eigen::MatrixXd& parti auto points = remove_bad_particles(particles); - vtkSmartPointer poly_data = MeshWarper::warp_mesh(points); + vtkSmartPointer poly_data = + laplacian_ready_ ? MeshWarper::warp_mesh_laplacian(points) : MeshWarper::warp_mesh(points); for (int i = 0; i < poly_data->GetNumberOfPoints(); i++) { double* p = poly_data->GetPoint(i); @@ -81,6 +87,7 @@ void MeshWarper::set_reference_mesh(vtkSmartPointer reference_mesh, // mark that the warp needs to be generated needs_warp_ = true; + laplacian_ready_ = false; warp_available_ = true; @@ -372,6 +379,7 @@ void MeshWarper::split_cell_on_edge(int cell_id, int new_vertex, int v0, int v1, //--------------------------------------------------------------------------- bool MeshWarper::generate_warp() { + TIME_SCOPE("MeshWarper::generate_warp"); if (is_contour_) { update_progress(1.0); needs_warp_ = false; @@ -414,6 +422,15 @@ bool MeshWarper::generate_warp() { faces_ = reference_mesh.faces(); // Perform warp + if (warp_method_ == WarpMethod::Laplacian) { + if (generate_laplacian_warp(vertices, faces_, vertices_)) { + update_progress(1.0); + needs_warp_ = false; + return true; + } + SW_DEBUG("Laplacian warp failed, falling back to biharmonic"); + } + if (MeshWarper::generate_warp_matrix(vertices, faces_, vertices_, warp_)) { // Success! update_progress(1.0); @@ -422,7 +439,7 @@ bool MeshWarper::generate_warp() { } } - SW_ERROR("Mesh Warp Error: igl:biharmonic_coordinates failed"); + SW_ERROR("Mesh Warp Error: warp generation failed"); // All attempts failed update_progress(1.0); @@ -433,6 +450,7 @@ bool MeshWarper::generate_warp() { //--------------------------------------------------------------------------- bool MeshWarper::generate_warp_matrix(Eigen::MatrixXd target_vertices, Eigen::MatrixXi target_faces, const Eigen::MatrixXd& references_vertices, Eigen::MatrixXd& warp) { + TIME_SCOPE("MeshWarper::generate_warp_matrix (biharmonic)"); Eigen::VectorXi closest_vertex_indices; Eigen::VectorXi vertex_indices = Eigen::VectorXi::LinSpaced(target_vertices.rows(), 0, target_vertices.rows() - 1); @@ -480,6 +498,7 @@ bool MeshWarper::find_landmarks_vertices_on_ref_mesh() { //--------------------------------------------------------------------------- vtkSmartPointer MeshWarper::warp_mesh(const Eigen::MatrixXd& points) { + TIME_SCOPE("MeshWarper::warp_mesh (biharmonic)"); int num_vertices = warp_.rows(); int num_faces = faces_.rows(); Eigen::MatrixXd v_out = warp_ * (points.rowwise() + Eigen::RowVector3d(0, 0, 0)); @@ -502,6 +521,135 @@ vtkSmartPointer MeshWarper::warp_mesh(const Eigen::MatrixXd& points return poly_data; } +//--------------------------------------------------------------------------- +int MeshWarper::get_num_warp_vertices() const { + if (laplacian_ready_) { + return laplacian_vertices_.rows(); + } + return warp_.rows(); +} + +//--------------------------------------------------------------------------- +bool MeshWarper::generate_laplacian_warp(const Eigen::MatrixXd& V, const Eigen::MatrixXi& F, + const Eigen::MatrixXd& control_points) { + TIME_SCOPE("MeshWarper::generate_laplacian_warp"); + // Map control points to closest mesh vertices + Eigen::VectorXi vertex_indices = Eigen::VectorXi::LinSpaced(V.rows(), 0, V.rows() - 1); + Eigen::VectorXd squared_distances; + Eigen::MatrixXd unused_closest_points; + Eigen::VectorXi closest_vertex_indices; + igl::point_mesh_squared_distance(control_points, V, vertex_indices, squared_distances, closest_vertex_indices, + unused_closest_points); + + // Deduplicate handle indices + std::set handle_set(closest_vertex_indices.data(), + closest_vertex_indices.data() + closest_vertex_indices.size()); + laplacian_handles_.resize(handle_set.size()); + int idx = 0; + for (int h : handle_set) { + laplacian_handles_(idx++) = h; + } + + // Build cotangent Laplacian + Eigen::SparseMatrix L; + igl::cotmatrix(V, F, L); + + // Q = L^T * L (bilaplacian energy) + Eigen::SparseMatrix Q = L.transpose() * L; + + // Precompute factorization with fixed handle constraints + Eigen::SparseMatrix Aeq; // empty equality constraints + if (!igl::min_quad_with_fixed_precompute(Q, laplacian_handles_, Aeq, true, mqwf_data_)) { + SW_DEBUG("Laplacian precompute failed"); + return false; + } + + // Store constant RHS: B = -Q * V_ref (the linear term for the unconstrained part) + laplacian_rhs_ = -Q * V; + + // Store reference vertices + laplacian_vertices_ = V; + + // Precompute handle_map_: a sparse matrix (num_handles x num_controls) that maps + // control point positions to handle Y rows, averaging when multiple controls map to the same handle. + // Also precompute handle_default_: reference positions for handles with no control point. + int num_handles = laplacian_handles_.size(); + int num_controls = closest_vertex_indices.size(); + + // Build a map from vertex id to handle index + std::unordered_map vertex_to_handle; + for (int i = 0; i < num_handles; i++) { + vertex_to_handle[laplacian_handles_(i)] = i; + } + + // Count how many controls map to each handle + Eigen::VectorXi handle_counts = Eigen::VectorXi::Zero(num_handles); + for (int i = 0; i < num_controls; i++) { + auto it = vertex_to_handle.find(closest_vertex_indices(i)); + if (it != vertex_to_handle.end()) { + handle_counts(it->second)++; + } + } + + // Build sparse matrix with weights = 1/count for averaging + std::vector> triplets; + triplets.reserve(num_controls); + for (int i = 0; i < num_controls; i++) { + auto it = vertex_to_handle.find(closest_vertex_indices(i)); + if (it != vertex_to_handle.end()) { + int h = it->second; + triplets.emplace_back(h, i, 1.0 / handle_counts(h)); + } + } + handle_map_.resize(num_handles, num_controls); + handle_map_.setFromTriplets(triplets.begin(), triplets.end()); + + // Default Y: reference positions for handles with no control point + handle_default_ = Eigen::MatrixXd::Zero(num_handles, 3); + for (int i = 0; i < num_handles; i++) { + if (handle_counts(i) == 0) { + handle_default_.row(i) = V.row(laplacian_handles_(i)); + } + } + + laplacian_ready_ = true; + return true; +} + +//--------------------------------------------------------------------------- +vtkSmartPointer MeshWarper::warp_mesh_laplacian(const Eigen::MatrixXd& points) { + TIME_SCOPE("MeshWarper::warp_mesh_laplacian"); + // Build Y from precomputed sparse map: Y = handle_map_ * points + handle_default_ + Eigen::MatrixXd Y = handle_map_ * points + handle_default_; + + Eigen::MatrixXd Beq; // empty + Eigen::MatrixXd Z; + if (!igl::min_quad_with_fixed_solve(mqwf_data_, laplacian_rhs_, Y, Beq, Z)) { + SW_ERROR("Laplacian solve failed"); + return warp_mesh(points); // fallback to biharmonic + } + + // Build VTK polydata from solved vertices + int num_vertices = Z.rows(); + int num_faces = faces_.rows(); + vtkSmartPointer poly_data = vtkSmartPointer::New(); + vtkSmartPointer out_points = vtkSmartPointer::New(); + out_points->SetNumberOfPoints(num_vertices); + for (vtkIdType i = 0; i < num_vertices; i++) { + out_points->SetPoint(i, Z(i, 0), Z(i, 1), Z(i, 2)); + } + vtkSmartPointer polys = vtkSmartPointer::New(); + for (vtkIdType i = 0; i < num_faces; i++) { + polys->InsertNextCell(3); + polys->InsertCellPoint(faces_(i, 0)); + polys->InsertCellPoint(faces_(i, 1)); + polys->InsertCellPoint(faces_(i, 2)); + } + poly_data->SetPoints(out_points); + poly_data->SetPolys(polys); + return poly_data; +} + //--------------------------------------------------------------------------- //--------------------------------------------------------------------------- void MeshWarper::diagnose_biharmonic_failure(const Eigen::MatrixXd& TV, const Eigen::MatrixXi& TF, diff --git a/Libs/Mesh/MeshWarper.h b/Libs/Mesh/MeshWarper.h index 610cb2e95d..70ce57433e 100644 --- a/Libs/Mesh/MeshWarper.h +++ b/Libs/Mesh/MeshWarper.h @@ -7,6 +7,7 @@ * The MeshWarper provides an object to warp meshes for surface reconstruction */ +#include #include #include @@ -14,13 +15,16 @@ namespace shapeworks { +//! Warp method selection for MeshWarper +enum class WarpMethod { Laplacian, Biharmonic }; + /** * \class MeshWarper * \ingroup Group-Mesh * * This class implements mesh warping based on correspondence particles. - * Correspondence points are embedded into the mesh as new vertices (traingles split). Then a biharmonic deformation - * is used to warp the mesh to new sets of correspondence particles. + * Correspondence points are embedded into the mesh as new vertices (triangles split). Then a deformation + * (Laplacian or biharmonic) is used to warp the mesh to new sets of correspondence particles. * * It can optionally be used to warp landmarks along with the mesh by embedding them as vertices * @@ -66,6 +70,15 @@ class MeshWarper { //! Return the reference particles const Eigen::MatrixXd& get_reference_particles() const { return this->reference_particles_; } + //! Set the warp method (Laplacian or Biharmonic) + void set_warp_method(WarpMethod method) { warp_method_ = method; } + + //! Get the current warp method + WarpMethod get_warp_method() const { return warp_method_; } + + //! Return the number of vertices in the warped mesh + int get_num_warp_vertices() const; + protected: //! For overriding to handle progress updates virtual void update_progress(float p) {} @@ -103,6 +116,12 @@ class MeshWarper { void diagnose_biharmonic_failure(const Eigen::MatrixXd& TV, const Eigen::MatrixXi& TF, const std::vector>& S, int k); + //! Generate the Laplacian warp data + bool generate_laplacian_warp(const Eigen::MatrixXd& V, const Eigen::MatrixXi& F, + const Eigen::MatrixXd& control_points); + + //! Warp mesh using Laplacian deformation + vtkSmartPointer warp_mesh_laplacian(const Eigen::MatrixXd& points); // Members Eigen::MatrixXi faces_; @@ -126,5 +145,17 @@ class MeshWarper { Eigen::MatrixXd reference_particles_; //! Whether the reference is a contour bool is_contour_ = false; + + //! Selected warp method + WarpMethod warp_method_ = WarpMethod::Biharmonic; + + //! Laplacian solve data + igl::min_quad_with_fixed_data mqwf_data_; + Eigen::MatrixXd laplacian_rhs_; //!< Constant RHS term: -Q * V_ref (n x 3) + Eigen::MatrixXd laplacian_vertices_; //!< Reference surface vertices + Eigen::VectorXi laplacian_handles_; //!< Handle vertex indices + Eigen::SparseMatrix handle_map_; //!< Maps control points to handle Y rows (num_handles x num_controls) + Eigen::MatrixXd handle_default_; //!< Default Y for handles with no control point (num_handles x 3) + bool laplacian_ready_ = false; }; } // namespace shapeworks diff --git a/Libs/Python/ShapeworksPython.cpp b/Libs/Python/ShapeworksPython.cpp index ba606bcbf1..a63f18aa39 100644 --- a/Libs/Python/ShapeworksPython.cpp +++ b/Libs/Python/ShapeworksPython.cpp @@ -1112,6 +1112,12 @@ PYBIND11_MODULE(shapeworks_py, m) { ; + // WarpMethod enum + py::enum_(m, "WarpMethod") + .value("Laplacian", WarpMethod::Laplacian) + .value("Biharmonic", WarpMethod::Biharmonic) + .export_values(); + // MeshWarping py::class_(m, "MeshWarper") @@ -1177,7 +1183,19 @@ PYBIND11_MODULE(shapeworks_py, m) { [](MeshWarper& w, const Mesh& warped_mesh) -> decltype(auto) { return w.extract_landmarks(warped_mesh.getVTKMesh()); }, - "Extract the landmarks from the warped mesh and return the landmarks (matrix [Nx3])", "warped_mesh"_a); + "Extract the landmarks from the warped mesh and return the landmarks (matrix [Nx3])", "warped_mesh"_a) + + .def( + "setWarpMethod", [](MeshWarper& w, WarpMethod method) { w.set_warp_method(method); }, + "Set the warp method (WarpMethod.Laplacian or WarpMethod.Biharmonic)", "method"_a) + + .def( + "getWarpMethod", [](MeshWarper& w) -> decltype(auto) { return w.get_warp_method(); }, + "Return the current warp method.") + + .def( + "getNumWarpVertices", [](MeshWarper& w) -> decltype(auto) { return w.get_num_warp_vertices(); }, + "Return the number of vertices in the warped mesh."); // MeshUtils py::class_(m, "MeshUtils") diff --git a/Studio/Analysis/AnalysisTool.cpp b/Studio/Analysis/AnalysisTool.cpp index 0e9d6b52f9..98c7e7690a 100644 --- a/Studio/Analysis/AnalysisTool.cpp +++ b/Studio/Analysis/AnalysisTool.cpp @@ -23,6 +23,7 @@ #include #include +#include #include "ParticleAreaPanel.h" #include "ShapeScalarPanel.h" @@ -36,6 +37,7 @@ const std::string AnalysisTool::MODE_SINGLE_SAMPLE_C("single sample"); const std::string AnalysisTool::MODE_REGRESSION_C("regression"); constexpr auto MESH_WARP_TEMPLATE_INDEX = "mesh_warp_template_index"; +constexpr auto MESH_WARP_METHOD = "mesh_warp_method"; //--------------------------------------------------------------------------- //! Helper to extract x,y,z from x,y,z,scalar @@ -205,6 +207,8 @@ AnalysisTool::AnalysisTool(Preferences& prefs) : preferences_(prefs) { connect(ui_->mesh_warp_median_button, &QPushButton::clicked, this, &AnalysisTool::mesh_warp_median_clicked); connect(ui_->mesh_warp_sample_spinbox, qOverload(&QSpinBox::valueChanged), this, &AnalysisTool::mesh_warp_sample_changed); + connect(ui_->warp_method_combo, qOverload(&QComboBox::currentIndexChanged), this, + &AnalysisTool::mesh_warp_method_changed); connect(ui_->mesh_warp_run_button, &QPushButton::clicked, this, &AnalysisTool::mesh_warp_run_clicked); } @@ -914,6 +918,7 @@ void AnalysisTool::load_settings() { QString::number(static_cast(params.get("network_pvalue_threshold", 0.05)))); ui_->mesh_warp_sample_spinbox->setValue(params.get(MESH_WARP_TEMPLATE_INDEX, -1)); + ui_->warp_method_combo->setCurrentIndex(params.get(MESH_WARP_METHOD, preferences_.get_warp_method())); update_group_boxes(); ui_->group_combo->setCurrentText(QString::fromStdString(params.get("current_group", "-None-"))); @@ -1718,6 +1723,8 @@ void AnalysisTool::initialize_mesh_warper() { SW_ERROR("Unable to set reference mesh, groomed mesh is unavailable"); return; } + WarpMethod warp_method = ui_->warp_method_combo->currentIndex() == 0 ? WarpMethod::Laplacian : WarpMethod::Biharmonic; + auto meshes = mesh_group.meshes(); for (int i = 0; i < mesh_group.meshes().size(); i++) { Eigen::VectorXd particles = median_shape->get_particles().get_local_particles(i); @@ -1729,7 +1736,9 @@ void AnalysisTool::initialize_mesh_warper() { Mesh mesh(poly_data); median_shape->get_constraints(i).clipMesh(mesh); - session_->get_mesh_manager()->get_mesh_warper(i)->set_reference_mesh(mesh.getVTKMesh(), points); + auto warper = session_->get_mesh_manager()->get_mesh_warper(i); + warper->set_warp_method(warp_method); + warper->set_reference_mesh(mesh.getVTKMesh(), points); } } } @@ -1870,6 +1879,7 @@ void AnalysisTool::handle_lda_complete() { group_lda_job_->plot(ui_->lda_graph, left_group, right_group); ui_->lda_graph->setVisible(true); ui_->lda_hint_label->setVisible(true); + QTimer::singleShot(0, this, &AnalysisTool::resize_tab_to_current); } void AnalysisTool::handle_network_analysis_progress(int progress) { @@ -1909,15 +1919,17 @@ void AnalysisTool::group_analysis_combo_changed() { } if (index == 0) { // none ui_->group_analysis_stacked_widget->setVisible(false); - ui_->group_analysis_stacked_widget->setEnabled(false); + ui_->group_analysis_stacked_widget->setMaximumHeight(0); } else { ui_->group_analysis_stacked_widget->setCurrentIndex(index - 1); + ui_->group_analysis_stacked_widget->setMaximumHeight(QWIDGETSIZE_MAX); ui_->group_analysis_stacked_widget->setVisible(true); - ui_->group_analysis_stacked_widget->setEnabled(true); if (ui_->group_analysis_stacked_widget->currentWidget() == ui_->lda_page) { update_lda_graph(); } } + // Recalculate tab height since analysis content changed + QTimer::singleShot(0, this, &AnalysisTool::resize_tab_to_current); Q_EMIT update_view(); } @@ -2094,6 +2106,13 @@ void AnalysisTool::mesh_warp_sample_changed() { ui_->template_mesh_name_label->setText(QString::fromStdString(shapes[index]->get_subject()->get_display_name())); } +//--------------------------------------------------------------------------- +void AnalysisTool::mesh_warp_method_changed() { + auto params = session_->get_project()->get_parameters(Parameters::ANALYSIS_PARAMS); + params.set(MESH_WARP_METHOD, ui_->warp_method_combo->currentIndex()); + session_->get_project()->set_parameters(Parameters::ANALYSIS_PARAMS, params); +} + //--------------------------------------------------------------------------- void AnalysisTool::mesh_warp_run_clicked() { session_->handle_clear_cache(); @@ -2105,6 +2124,22 @@ void AnalysisTool::mesh_warp_run_clicked() { void AnalysisTool::handle_tab_changed() { stats_ready_ = false; compute_stats(); + + resize_tab_to_current(); +} + +void AnalysisTool::resize_tab_to_current() { + auto* current = ui_->tabWidget->currentWidget(); + if (!current) return; + + // Force Qt to recalculate layout so minimumSize reflects hidden/shown widgets + current->layout()->activate(); + int h = current->layout()->minimumSize().height(); + int tab_bar = ui_->tabWidget->tabBar()->sizeHint().height(); + // Add frame margins + int frame = ui_->tabWidget->height() - ui_->tabWidget->currentWidget()->height(); + if (frame < tab_bar) frame = tab_bar + 6; + ui_->tabWidget->setMaximumHeight(frame + h); } //--------------------------------------------------------------------------- @@ -2148,6 +2183,11 @@ void AnalysisTool::set_active(bool active) { } active_ = active; update_interface(); + if (active) { + handle_tab_changed(); + // Deferred recalculation after layout settles + QTimer::singleShot(0, this, &AnalysisTool::handle_tab_changed); + } } //--------------------------------------------------------------------------- diff --git a/Studio/Analysis/AnalysisTool.h b/Studio/Analysis/AnalysisTool.h index 33842f2fc7..c6f0dd8511 100644 --- a/Studio/Analysis/AnalysisTool.h +++ b/Studio/Analysis/AnalysisTool.h @@ -207,6 +207,7 @@ class AnalysisTool : public QWidget { // mesh warping options void mesh_warp_median_clicked(); void mesh_warp_sample_changed(); + void mesh_warp_method_changed(); void mesh_warp_run_clicked(); void handle_tab_changed(); @@ -221,6 +222,7 @@ class AnalysisTool : public QWidget { private: void compute_reconstructed_domain_transforms(); + void resize_tab_to_current(); bool active_ = false; diff --git a/Studio/Analysis/AnalysisTool.ui b/Studio/Analysis/AnalysisTool.ui index 2ac40def64..cb2237ba5a 100644 --- a/Studio/Analysis/AnalysisTool.ui +++ b/Studio/Analysis/AnalysisTool.ui @@ -1096,19 +1096,6 @@ QWidget#particles_panel { - - - - Qt::Vertical - - - - 20 - 40 - - - - @@ -1225,6 +1212,12 @@ QWidget#particles_panel { + + + 0 + 0 + + Groups @@ -1244,7 +1237,11 @@ QWidget#particles_panel { - + + + QAbstractScrollArea::AdjustToContents + + @@ -1678,19 +1675,6 @@ QWidget#particles_panel { - - - - Qt::Vertical - - - - 20 - 0 - - - - @@ -2202,7 +2186,7 @@ Reference Domain - + Surface reconstruction based on distance transforms and particles. @@ -2212,33 +2196,7 @@ Reference Domain - - - - Qt::Horizontal - - - - 40 - 20 - - - - - - - - Qt::Horizontal - - - - 40 - 20 - - - - - + Legacy surface reconstruction method that requires only particles @@ -2284,6 +2242,9 @@ Reference Domain 16777215 + + Qt::AlignCenter + 3 @@ -2316,6 +2277,9 @@ Reference Domain 16777215 + + Qt::AlignCenter + 90.000000000000000 @@ -2344,6 +2308,9 @@ Reference Domain 16777215 + + Qt::AlignCenter + 0.010000000000000 @@ -2385,7 +2352,7 @@ Reference Domain - + Surface reconstruction method that uses mesh warping from a median mesh and particle positions. @@ -2417,20 +2384,21 @@ Reference Domain - Template Sample: + Template: - - - - - - true + + + Qt::AlignCenter + + + + - Run Mesh Warp + Template Name: @@ -2441,10 +2409,34 @@ Reference Domain - - + + - Template Name: + Warp Method: + + + + + + + + Laplacian + + + + + Biharmonic + + + + + + + + true + + + Run Mesh Warp @@ -2599,6 +2591,9 @@ Reference Domain 16777215 + + Qt::AlignCenter + 180.000000000000000 diff --git a/Studio/Data/Preferences.cpp b/Studio/Data/Preferences.cpp index 8f1da93a96..63bba1c7ee 100644 --- a/Studio/Data/Preferences.cpp +++ b/Studio/Data/Preferences.cpp @@ -387,6 +387,12 @@ bool Preferences::get_telemetry_asked() { return settings_.value("General/teleme //----------------------------------------------------------------------------- void Preferences::set_telemetry_asked(bool asked) { settings_.setValue("General/telemetry_asked", asked); } +//----------------------------------------------------------------------------- +int Preferences::get_warp_method() { return settings_.value("Analysis/warp_method", 1).toInt(); } + +//----------------------------------------------------------------------------- +void Preferences::set_warp_method(int method) { settings_.setValue("Analysis/warp_method", method); } + //----------------------------------------------------------------------------- QStringList Preferences::get_pending_telemetry_events() { return settings_.value("Telemetry/pending_events").toStringList(); diff --git a/Studio/Data/Preferences.h b/Studio/Data/Preferences.h index ed60fb17da..6591bebf8d 100644 --- a/Studio/Data/Preferences.h +++ b/Studio/Data/Preferences.h @@ -130,6 +130,9 @@ class Preferences : public QObject { bool get_telemetry_asked(); void set_telemetry_asked(bool asked); + int get_warp_method(); + void set_warp_method(int method); + QStringList get_pending_telemetry_events(); void set_pending_telemetry_events(QStringList events); diff --git a/Studio/Data/PreferencesWindow.cpp b/Studio/Data/PreferencesWindow.cpp index 10ce75cf18..1729ab1ae4 100644 --- a/Studio/Data/PreferencesWindow.cpp +++ b/Studio/Data/PreferencesWindow.cpp @@ -54,6 +54,8 @@ PreferencesWindow::PreferencesWindow(QWidget* parent, Preferences& prefs) : pref connect(ui_->discrete_color_mode, &QCheckBox::toggled, this, &PreferencesWindow::save_to_preferences); connect(ui_->reverse_color_map, &QCheckBox::toggled, this, &PreferencesWindow::save_to_preferences); connect(ui_->auto_update_checkbox, &QCheckBox::toggled, this, &PreferencesWindow::save_to_preferences); + connect(ui_->warp_method, qOverload(&QComboBox::currentIndexChanged), this, + &PreferencesWindow::save_to_preferences); } //----------------------------------------------------------------------------- @@ -117,6 +119,7 @@ void PreferencesWindow::set_values_from_preferences() { ui_->auto_update_checkbox->setChecked(preferences_.get_auto_update_check()); ui_->telemetry_enabled->setChecked(preferences_.get_telemetry_enabled()); ui_->data_loader_num_workers->setText(QString::number(preferences_.get_dataloader_num_workers())); + ui_->warp_method->setCurrentIndex(preferences_.get_warp_method()); update_labels(); } @@ -149,6 +152,7 @@ void PreferencesWindow::save_to_preferences() { preferences_.set_auto_update_check(ui_->auto_update_checkbox->isChecked()); preferences_.set_telemetry_enabled(ui_->telemetry_enabled->isChecked()); preferences_.set_dataloader_num_workers(ui_->data_loader_num_workers->text().toInt()); + preferences_.set_warp_method(ui_->warp_method->currentIndex()); update_labels(); Q_EMIT update_view(); } diff --git a/Studio/Data/PreferencesWindow.ui b/Studio/Data/PreferencesWindow.ui index 9c79973764..a565aca98d 100644 --- a/Studio/Data/PreferencesWindow.ui +++ b/Studio/Data/PreferencesWindow.ui @@ -522,7 +522,28 @@ - + + + + Default Warp Method: + + + + + + + + Laplacian + + + + + Biharmonic + + + + + Qt::Vertical diff --git a/Testing/MeshTests/MeshTests.cpp b/Testing/MeshTests/MeshTests.cpp index d7d343e090..050fb81429 100644 --- a/Testing/MeshTests/MeshTests.cpp +++ b/Testing/MeshTests/MeshTests.cpp @@ -1,3 +1,4 @@ +#include #include #include "Image.h" @@ -715,6 +716,7 @@ void mesh_warp_test(std::string ref_mesh, std::string ref_particles, std::string movingPoints.transposeInPlace(); MeshWarper warper; + warper.set_warp_method(WarpMethod::Biharmonic); warper.set_reference_mesh(reference.getVTKMesh(), staticPoints); ASSERT_TRUE(warper.get_warp_available()); @@ -750,6 +752,149 @@ TEST(MeshTests, warpTest4) { // "/mesh_warp/lv_shared2_baseline.vtk"); //} +// Laplacian warp tests +TEST(MeshTests, laplacianWarpBasic) { + // Test that Laplacian warp succeeds and produces a valid mesh + Mesh reference(std::string(TEST_DATA_DIR) + "/ellipsoid_0.ply"); + + std::string staticPath = std::string(TEST_DATA_DIR) + "/ellipsoid_0.particles"; + std::string movingPath = std::string(TEST_DATA_DIR) + "/ellipsoid_1.particles"; + std::vector paths = {staticPath, movingPath}; + ParticleSystemEvaluation particlesystem(paths); + Eigen::MatrixXd allPts = particlesystem.get_matrix(); + Eigen::MatrixXd staticPoints = allPts.col(0); + Eigen::MatrixXd movingPoints = allPts.col(1); + + int numParticles = staticPoints.rows() / 3; + staticPoints.resize(3, numParticles); + staticPoints.transposeInPlace(); + movingPoints.resize(3, numParticles); + movingPoints.transposeInPlace(); + + MeshWarper warper; + warper.set_warp_method(WarpMethod::Laplacian); + warper.set_reference_mesh(reference.getVTKMesh(), staticPoints); + ASSERT_TRUE(warper.get_warp_available()); + + Mesh output = warper.build_mesh(movingPoints); + ASSERT_TRUE(output.numPoints() > 0); + ASSERT_TRUE(output.numFaces() > 0); + + // Verify no NaN values in output + Eigen::MatrixXd points = output.points(); + for (int i = 0; i < points.rows(); i++) { + ASSERT_FALSE(std::isnan(points(i, 0))); + ASSERT_FALSE(std::isnan(points(i, 1))); + ASSERT_FALSE(std::isnan(points(i, 2))); + } +} + +TEST(MeshTests, laplacianWarpSelfWarp) { + // Self-warp: warping with same particles should produce mesh close to reference + Mesh reference(std::string(TEST_DATA_DIR) + "/ellipsoid_0.ply"); + + std::string staticPath = std::string(TEST_DATA_DIR) + "/ellipsoid_0.particles"; + std::vector paths = {staticPath, staticPath}; + ParticleSystemEvaluation particlesystem(paths); + Eigen::MatrixXd allPts = particlesystem.get_matrix(); + Eigen::MatrixXd staticPoints = allPts.col(0); + + int numParticles = staticPoints.rows() / 3; + staticPoints.resize(3, numParticles); + staticPoints.transposeInPlace(); + + MeshWarper warper; + warper.set_warp_method(WarpMethod::Laplacian); + warper.set_reference_mesh(reference.getVTKMesh(), staticPoints); + ASSERT_TRUE(warper.get_warp_available()); + + // Self-warp: use same particles as reference + Mesh output = warper.build_mesh(staticPoints); + ASSERT_TRUE(output.numPoints() > 0); + + // Vertex count should match + ASSERT_EQ(output.numPoints(), warper.get_num_warp_vertices()); +} + +TEST(MeshTests, warpMethodBenchmark) { + // Benchmark: run both methods on the same data and print timing comparison + std::string ref_mesh = "/mesh_warp/mesh_warp2.vtk"; + std::string ref_particles = "/mesh_warp/mesh_warp2.particles"; + + Mesh reference(std::string(TEST_DATA_DIR) + ref_mesh); + + std::string staticPath = std::string(TEST_DATA_DIR) + ref_particles; + std::vector paths = {staticPath, staticPath}; + ParticleSystemEvaluation particlesystem(paths); + Eigen::MatrixXd allPts = particlesystem.get_matrix(); + Eigen::MatrixXd staticPoints = allPts.col(0); + int numParticles = staticPoints.rows() / 3; + staticPoints.resize(3, numParticles); + staticPoints.transposeInPlace(); + + const int num_warps = 10; + + // -- Laplacian -- + { + MeshWarper warper; + warper.set_warp_method(WarpMethod::Laplacian); + warper.set_reference_mesh(reference.getVTKMesh(), staticPoints); + + auto t0 = std::chrono::high_resolution_clock::now(); + warper.generate_warp(); + auto t1 = std::chrono::high_resolution_clock::now(); + double precompute_ms = std::chrono::duration(t1 - t0).count(); + + auto t2 = std::chrono::high_resolution_clock::now(); + for (int i = 0; i < num_warps; i++) { + warper.build_mesh(staticPoints); + } + auto t3 = std::chrono::high_resolution_clock::now(); + double warp_ms = std::chrono::duration(t3 - t2).count(); + + std::cout << "\n=== Warp Method Benchmark ===" << std::endl; + std::cout << "Laplacian: precompute=" << precompute_ms << "ms" + << " warp=" << warp_ms / num_warps << "ms/shape" + << " (" << num_warps << " warps)" << std::endl; + } + + // -- Biharmonic -- + { + MeshWarper warper; + warper.set_warp_method(WarpMethod::Biharmonic); + warper.set_reference_mesh(reference.getVTKMesh(), staticPoints); + + auto t0 = std::chrono::high_resolution_clock::now(); + warper.generate_warp(); + auto t1 = std::chrono::high_resolution_clock::now(); + double precompute_ms = std::chrono::duration(t1 - t0).count(); + + auto t2 = std::chrono::high_resolution_clock::now(); + for (int i = 0; i < num_warps; i++) { + warper.build_mesh(staticPoints); + } + auto t3 = std::chrono::high_resolution_clock::now(); + double warp_ms = std::chrono::duration(t3 - t2).count(); + + std::cout << "Biharmonic: precompute=" << precompute_ms << "ms" + << " warp=" << warp_ms / num_warps << "ms/shape" + << " (" << num_warps << " warps)" << std::endl; + std::cout << "==============================\n" << std::endl; + } +} + +TEST(MeshTests, laplacianWarpMethodSwitch) { + // Verify that warp method can be set and queried + MeshWarper warper; + ASSERT_EQ(warper.get_warp_method(), WarpMethod::Biharmonic); // default + + warper.set_warp_method(WarpMethod::Biharmonic); + ASSERT_EQ(warper.get_warp_method(), WarpMethod::Biharmonic); + + warper.set_warp_method(WarpMethod::Laplacian); + ASSERT_EQ(warper.get_warp_method(), WarpMethod::Laplacian); +} + TEST(MeshTests, findReferenceMeshTest) { std::vector meshes; meshes.push_back(Mesh(std::string(TEST_DATA_DIR) + "/m03_L_femur.ply")); diff --git a/docs/studio/images/surface_reconstruction_biharmonic.jpg b/docs/studio/images/surface_reconstruction_biharmonic.jpg new file mode 100644 index 0000000000..5e840b0a26 Binary files /dev/null and b/docs/studio/images/surface_reconstruction_biharmonic.jpg differ diff --git a/docs/studio/images/surface_reconstruction_laplacian.jpg b/docs/studio/images/surface_reconstruction_laplacian.jpg new file mode 100644 index 0000000000..fb5235f765 Binary files /dev/null and b/docs/studio/images/surface_reconstruction_laplacian.jpg differ diff --git a/docs/studio/images/surface_reconstruction_panel.png b/docs/studio/images/surface_reconstruction_panel.png new file mode 100644 index 0000000000..fc2a4a1aba Binary files /dev/null and b/docs/studio/images/surface_reconstruction_panel.png differ diff --git a/docs/studio/studio-analyze.md b/docs/studio/studio-analyze.md index 5540aa29ab..66727a7d22 100644 --- a/docs/studio/studio-analyze.md +++ b/docs/studio/studio-analyze.md @@ -101,19 +101,7 @@ See [Shape Model Evaluation](../new/ssm-eval.md) for more information about shap ## Surface Reconstruction Panel -The surface reconstruction panel provides options for the surface reconstruction method. There are three surface reconstruction methods available depending on the data you supply. - -`Legacy` - If an older XML file with only particle files is supplied, then only this option is available. This is the fallback option since it requires only the particles. - -`Mesh Warping Based` - Mesh warping based method that utilizes the mean mesh. You must have either meshes supplied or image based (distance transforms). - -`Distance Transform Based` - Surface reconstruction based on distance transforms. Project must have distance transforms. - -![ShapeWorks Studio Analysis Surface Reconstruction Panel](../img/studio/studio_analyze_surface_reconstruction.png) - -Below is an example of the difference in using this option. - -![ShapeWorks Studio Surface Reconstruction Result](../img/studio/studio_analyze_surface_reconstruction_result.png) +See the [Surface Reconstruction](surface-reconstruction.md) page for details on reconstruction methods, warp method options, and examples. ## Good/Bad Particles Panel ## diff --git a/docs/studio/surface-reconstruction.md b/docs/studio/surface-reconstruction.md new file mode 100644 index 0000000000..501479e064 --- /dev/null +++ b/docs/studio/surface-reconstruction.md @@ -0,0 +1,81 @@ +# Surface Reconstruction + +ShapeWorks reconstructs dense surface meshes from sparse correspondence particles. Since correspondence particles are placed sparsely across the surface, a dense mesh must be generated for visualization and downstream analysis. ShapeWorks does this by deforming a reference (template) mesh to match each shape's particle positions. + +## Reconstruction Methods + +Three surface reconstruction methods are available in the Analysis module, depending on your data: + +| Method | Requirements | Description | +|--------|-------------|-------------| +| **Mesh Warping** | Meshes or distance transforms | Deforms a template mesh using correspondence particles as control points. Recommended for most workflows. | +| **Distance Transform** | Distance transforms | Reconstructs surfaces from distance transforms. Requires groomed distance transforms in the project. | +| **Legacy** | Particles only | Fallback method for older XML projects that only contain particle files. | + +The reconstruction method can be selected from the radio buttons in the Surface Reconstruction panel. + +![ShapeWorks Studio Analysis Surface Reconstruction Panel](../img/studio/studio_analyze_surface_reconstruction.png) + +Below is an example of the difference between reconstruction methods. + +![ShapeWorks Studio Surface Reconstruction Result](../img/studio/studio_analyze_surface_reconstruction_result.png) + +## Mesh Warping + +Mesh warping is the default and recommended reconstruction method. It works by: + +1. Selecting a **template shape** (typically the median) from the population. +2. **Embedding** the template's correspondence particles as vertices in its mesh (splitting triangles as needed). +3. **Deforming** the enriched template mesh so that the embedded particle vertices move to each target shape's particle positions, dragging the rest of the mesh along smoothly. + +This produces a dense mesh for every shape in the population that shares the same connectivity (triangle topology), which is useful for computing point-to-point correspondences across the full surface. + +### Template Selection + +The template shape serves as the reference mesh that gets deformed to reconstruct all other shapes. You can choose the template using the options in the Mesh Warp panel: + +- **Median** -- Automatically selects the shape closest to the population median. This is the default and generally works well. +- **Template Sample** -- Manually select a specific shape index to use as the template. + +After changing the template or warp method, click **Run Mesh Warp** to regenerate the reconstructions. + +![Mesh Warp Options Panel](images/surface_reconstruction_panel.png){: width="400" } + +### Warp Methods + +Two warp methods are available for deforming the template mesh. The method can be selected from the **Warp Method** dropdown in the Mesh Warp options. The warp method is saved per project. The default for new projects can be changed in **Edit > Preferences** under the Mesh Optimization section. + +#### Biharmonic (Default) + +Biharmonic deformation computes a weight matrix that expresses every mesh vertex as a weighted combination of the particle positions. Once computed, each per-shape warp is a single matrix multiplication. This is the default and works well for most cases. + +#### Laplacian + +Laplacian surface deformation preserves the local curvature of the reference mesh during warping. Vertices far from any particle maintain their original surface shape, which can reduce artificial thinning in regions with sparse particle coverage. + +Both methods use [libigl](https://libigl.github.io/) for the underlying geometry processing (biharmonic weights and cotangent Laplacian respectively). + +**Laplacian** may produce better results when particles are unevenly distributed or on thin structures. **Biharmonic** is a good default for most datasets. + +| Biharmonic | Laplacian | +|:---:|:---:| +| ![Biharmonic](images/surface_reconstruction_biharmonic.jpg) | ![Laplacian](images/surface_reconstruction_laplacian.jpg) | + +### Python API + +The warp method can also be set programmatically using the Python API: + +```python +import shapeworks + +warper = shapeworks.MeshWarper() + +# Set warp method (default is Biharmonic) +warper.setWarpMethod(shapeworks.WarpMethod.Laplacian) +# or +warper.setWarpMethod(shapeworks.WarpMethod.Biharmonic) + +# Set up and run warping +warper.generateWarp(reference_mesh, reference_particles) +warped_mesh = warper.buildMesh(target_particles) +``` diff --git a/mkdocs.yml b/mkdocs.yml index 5e6a702c02..e6bd100b20 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -46,6 +46,7 @@ nav: - 'Groom Module': 'studio/studio-groom.md' - 'Optimize Module': 'studio/studio-optimize.md' - 'Analyze Module': 'studio/studio-analyze.md' + - 'Surface Reconstruction': 'studio/surface-reconstruction.md' - 'DeepSSM Module': 'studio/deepssm-in-studio.md' - 'Multiple Domains SSM': 'studio/multiple-domains.md' - 'Shared Boundaries': 'studio/studio-shared-boundary.md'