Merge pull request #2535 from SCIInstitute/amorris/fix-shared-boundary

Fix #2534 - shared boundary issue: skip extract-largest-component for shared boundary domains

Author: akenmorris

Libs/Groom/Groom.cpp

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

Libs/Groom/GroomParameters.cpp

@@ -544,6 +544,19 @@ bool GroomParameters::get_shared_boundaries_enabled() {
 //---------------------------------------------------------------------------
 void GroomParameters::set_shared_boundaries_enabled(bool enabled) { params_.set(Keys::SHARED_BOUNDARY, enabled); }
 
+//---------------------------------------------------------------------------
+bool GroomParameters::is_shared_boundary_domain() {
+  if (!get_shared_boundaries_enabled()) {
+    return false;
+  }
+  for (const auto& boundary : get_shared_boundaries()) {
+    if (boundary.first_domain == domain_name_ || boundary.second_domain == domain_name_) {
+      return true;
+    }
+  }
+  return false;
+}
+
 //---------------------------------------------------------------------------
 std::vector<GroomParameters::SharedBoundary> GroomParameters::get_shared_boundaries() {
   std::vector<SharedBoundary> boundaries;
@@ -576,7 +589,9 @@ std::vector<GroomParameters::SharedBoundary> GroomParameters::get_shared_boundar
     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);
+    if (!boundary.first_domain.empty() && !boundary.second_domain.empty()) {
+      boundaries.push_back(boundary);
+    }
 
     // Migrate to new format
     set_shared_boundaries(boundaries);

Libs/Groom/GroomParameters.h

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

Libs/Mesh/MeshUtils.cpp

@@ -266,20 +266,31 @@ int MeshUtils::findReferenceMesh(std::vector<Mesh>& meshes, int random_subset_si
  */
 
 //---------------------------------------------------------------------------
+// Robust orientation test: compute the loop's signed-area vector by summing edge
+// cross products. The dominant axis of this vector gives the loop's normal; its sign
+// determines orientation. Uses ALL vertices, so it is insensitive to where loop[0]
+// happens to start and does not suffer from atan2 branch-cut flips.
 static bool is_clockwise(const Eigen::MatrixXd& V, const Eigen::MatrixXi& F, const std::vector<int>& loop) {
   Eigen::RowVector3d centroid{0.0, 0.0, 0.0};
   for (const auto& i : loop) {
     centroid += V.row(i);
   }
   centroid /= loop.size();
 
-  // todo this is arbitrary and works for the peanut data and initial tests on LA+Septum data
-  // it enforces a consistent ordering in the boundary loop
-  const auto v0 = V.row(loop[0]) - centroid;
-  const auto v1 = V.row(loop[1]) - centroid;
-  const double angle0 = atan2(v0.z(), v0.y());
-  const double angle1 = atan2(v1.z(), v1.y());
-  return angle0 > angle1;
+  Eigen::RowVector3d area_vec{0.0, 0.0, 0.0};
+  for (size_t i = 0; i < loop.size(); i++) {
+    const Eigen::RowVector3d v0 = V.row(loop[i]) - centroid;
+    const Eigen::RowVector3d v1 = V.row(loop[(i + 1) % loop.size()]) - centroid;
+    area_vec += v0.cross(v1);
+  }
+
+  // Use the dominant axis of the area vector as the reference normal. Calling "clockwise"
+  // the case where the area vector points in the negative dominant-axis direction yields
+  // consistent results across samples that share the same boundary plane.
+  int max_axis = 0;
+  if (std::abs(area_vec[1]) > std::abs(area_vec[max_axis])) max_axis = 1;
+  if (std::abs(area_vec[2]) > std::abs(area_vec[max_axis])) max_axis = 2;
+  return area_vec[max_axis] < 0;
 }
 
 //---------------------------------------------------------------------------
@@ -293,7 +304,32 @@ Mesh MeshUtils::extract_boundary_loop(Mesh mesh) {
     throw std::runtime_error("Expected at least one boundary loop in the mesh");
   }
 
-  const auto& loop = loops[0];
+  auto loop = loops[0];  // copy so we can rotate it to a canonical start vertex
+
+  // Rotate the loop so it always starts at a canonical vertex (the one with the
+  // largest Y coordinate; lexicographic tiebreakers on Z then X). igl::boundary_loop
+  // returns loops starting at an arbitrary vertex, which produces per-subject
+  // rotational offsets that destroy inter-subject correspondence on contour domains.
+  {
+    size_t canonical = 0;
+    for (size_t i = 1; i < loop.size(); i++) {
+      const double cur_y = V(loop[i], 1);
+      const double cur_z = V(loop[i], 2);
+      const double cur_x = V(loop[i], 0);
+      const double best_y = V(loop[canonical], 1);
+      const double best_z = V(loop[canonical], 2);
+      const double best_x = V(loop[canonical], 0);
+      if (cur_y > best_y ||
+          (cur_y == best_y && cur_z > best_z) ||
+          (cur_y == best_y && cur_z == best_z && cur_x > best_x)) {
+        canonical = i;
+      }
+    }
+    if (canonical != 0) {
+      std::rotate(loop.begin(), loop.begin() + canonical, loop.end());
+    }
+  }
+
   const auto is_cw = is_clockwise(V, F, loop);
 
   auto pts = vtkSmartPointer<vtkPoints>::New();
@@ -389,6 +425,12 @@ static std::tuple<Eigen::MatrixXd, Eigen::MatrixXi, Eigen::MatrixXd, Eigen::Matr
   Eigen::MatrixXi out_F;
   Eigen::MatrixXd rem_V;
   Eigen::MatrixXi rem_F;
+
+  // If either mesh is empty, there can be no shared surface
+  if (is_empty(src_V, src_F) || is_empty(other_V, other_F)) {
+    return std::make_tuple(out_V, out_F, src_V, src_F);
+  }
+
   igl::AABB<Eigen::MatrixXd, 3> tree;
   tree.init(other_V, other_F);
 
@@ -482,6 +524,10 @@ std::array<Mesh, 3> MeshUtils::shared_boundary_extractor(const Mesh& mesh_l, con
   V_r = mesh_r.points();
   F_r = mesh_r.faces();
 
+  if (is_empty(V_l, F_l) || is_empty(V_r, F_r)) {
+    throw std::runtime_error("Input mesh is empty. Cannot extract shared boundary from empty meshes");
+  }
+
   Eigen::MatrixXd shared_V_l, shared_V_r, rem_V_l, rem_V_r;
   Eigen::MatrixXi shared_F_l, shared_F_r, rem_F_l, rem_F_r;
   std::tie(shared_V_l, shared_F_l, rem_V_l, rem_F_l) = find_shared_surface(V_l, F_l, V_r, F_r, tol);
@@ -1073,18 +1119,23 @@ vtkSmartPointer<vtkPolyData> MeshUtils::recreate_mesh(vtkSmartPointer<vtkPolyDat
 }
 
 //---------------------------------------------------------------------------
-vtkSmartPointer<vtkPolyData> MeshUtils::repair_mesh(vtkSmartPointer<vtkPolyData> mesh) {
+vtkSmartPointer<vtkPolyData> MeshUtils::repair_mesh(vtkSmartPointer<vtkPolyData> mesh, bool extract_largest) {
   auto triangle_filter = vtkSmartPointer<vtkTriangleFilter>::New();
   triangle_filter->SetInputData(mesh);
   triangle_filter->PassLinesOff();
   triangle_filter->Update();
 
-  auto connectivity = vtkSmartPointer<vtkPolyDataConnectivityFilter>::New();
-  connectivity->SetInputConnection(triangle_filter->GetOutputPort());
-  connectivity->SetExtractionModeToLargestRegion();
-  connectivity->Update();
+  vtkSmartPointer<vtkPolyData> triangulated = triangle_filter->GetOutput();
+
+  if (extract_largest) {
+    auto connectivity = vtkSmartPointer<vtkPolyDataConnectivityFilter>::New();
+    connectivity->SetInputData(triangulated);
+    connectivity->SetExtractionModeToLargestRegion();
+    connectivity->Update();
+    triangulated = connectivity->GetOutput();
+  }
 
-  auto cleaned = MeshUtils::clean_mesh(connectivity->GetOutput());
+  auto cleaned = MeshUtils::clean_mesh(triangulated);
 
   auto fixed = Mesh(cleaned).fixNonManifold();
 

Libs/Mesh/MeshUtils.h

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

Libs/Optimize/Domain/ContourDomain.cpp

@@ -344,13 +344,13 @@ int ContourDomain::NumberOfLinesIncidentOnPoint(int i) const {
 }
 
 void ContourDomain::ComputeAvgEdgeLength() {
-  const double total_length = std::accumulate(lines_.begin(), lines_.end(),
+  total_length_ = std::accumulate(lines_.begin(), lines_.end(),
                                               0.0, [&](double s, const vtkSmartPointer<vtkLine> &line) {
             const auto pt_a = GetPoint(line->GetPointId(0));
             const auto pt_b = GetPoint(line->GetPointId(1));
             return s + (pt_a - pt_b).norm();
           });
-  avg_edge_length_ = total_length / lines_.size();
+  avg_edge_length_ = total_length_ / lines_.size();
 }
 
 }

Libs/Optimize/Domain/ContourDomain.h

@@ -86,8 +86,11 @@ class ContourDomain : public ParticleDomain {
   }
 
   double GetSurfaceArea() const override {
-    // TODO: Implement something analogous for scaling purposes
-    return 1.0;
+    // Return length² as an area-equivalent for a contour so it participates in the
+    // sampling-scale auto-scaling. For a circle of radius r this is (2πr)² = 4π²r²,
+    // which is π times the matching sphere's surface area — comparable magnitude so
+    // the scale factor doesn't crush the contour gradient.
+    return total_length_ * total_length_;
   }
 
   void DeleteImages() override {
@@ -132,6 +135,7 @@ class ContourDomain : public ParticleDomain {
   mutable double geo_lq_dist_ = -1;
 
   double avg_edge_length_{0.0};
+  double total_length_{0.0};
 
   void ComputeBounds();
   void ComputeGeodesics(vtkSmartPointer<vtkPolyData> poly_data);

Libs/Optimize/Function/SamplingFunction.cpp

@@ -250,7 +250,10 @@ SamplingFunction::VectorType SamplingFunction::evaluate(unsigned int idx, unsign
 
   gradE = gradE / m_avgKappa;
 
-  // Apply sampling scale if enabled
+  // Apply sampling scale if enabled. Contour domains return length² from GetSurfaceArea,
+  // giving a scale factor comparable in magnitude to the equivalent mesh's — without this,
+  // contour particles would move at native gradient magnitudes (~1000× stronger than scaled
+  // mesh particles), causing them to spin rapidly around the loop instead of settling.
   if (m_SamplingScale) {
     double scale_factor = 1.0;
 

Libs/Optimize/OptimizeParameters.cpp

@@ -328,6 +328,9 @@ std::vector<std::vector<itk::Point<double>>> OptimizeParameters::get_initial_poi
     for (auto s : subjects) {
       if (s->is_fixed()) {
         count++;
+        if (d >= s->get_world_particle_filenames().size()) {
+          throw std::runtime_error("Subject " + s->get_display_name() + " does not have enough world particle files");
+        }
         // read the world points that are in the shared coordinate space
         auto filename = s->get_world_particle_filenames()[d];
         auto particles = read_particles_as_vector(filename);

Studio/Optimize/OptimizeTool.cpp

@@ -2,8 +2,11 @@
 
 // qt
 #include <QFileDialog>
+#include <QFontMetrics>
 #include <QIntValidator>
+#include <QLabel>
 #include <QMessageBox>
+#include <QResizeEvent>
 #include <QThread>
 #include <QTimer>
 
@@ -22,6 +25,47 @@
 
 using namespace shapeworks;
 
+namespace {
+// QLabel that elides its text with "..." when it's narrower than the text.
+// The full text is kept as a tooltip so the user can see it on hover.
+// sizeHint() reports the full-text width so the layout gives the label its natural
+// space when available; elision only happens when the layout must shrink it.
+class ElidedLabel : public QLabel {
+ public:
+  explicit ElidedLabel(const QString& text, QWidget* parent = nullptr) : QLabel(parent), full_text_(text) {
+    setToolTip(text);
+    updateElidedText();
+  }
+
+  QSize sizeHint() const override {
+    QFontMetrics metrics(font());
+    return QSize(metrics.horizontalAdvance(full_text_) + 4, QLabel::sizeHint().height());
+  }
+
+  QSize minimumSizeHint() const override {
+    // Minimum: enough for ellipsis plus a couple characters so the label can shrink
+    // further if the panel is very narrow, but not to zero width.
+    QFontMetrics metrics(font());
+    return QSize(metrics.horizontalAdvance(QStringLiteral("X...")), QLabel::minimumSizeHint().height());
+  }
+
+ protected:
+  void resizeEvent(QResizeEvent* event) override {
+    QLabel::resizeEvent(event);
+    updateElidedText();
+  }
+
+ private:
+  void updateElidedText() {
+    QFontMetrics metrics(font());
+    QString elided = metrics.elidedText(full_text_, Qt::ElideRight, width());
+    QLabel::setText(elided);
+  }
+
+  QString full_text_;
+};
+}  // namespace
+
 //---------------------------------------------------------------------------
 OptimizeTool::OptimizeTool(Preferences& prefs, Telemetry& telemetry) : preferences_(prefs), telemetry_(telemetry) {
   ui_ = new Ui_OptimizeTool;
@@ -473,6 +517,8 @@ void OptimizeTool::setup_domain_boxes() {
     QLineEdit* box = new QLineEdit(this);
     last_box = box;
     box->setAlignment(Qt::AlignHCenter);
+    box->setMinimumWidth(50);
+    box->setMaximumWidth(100);
     box->setValidator(above_zero);
     box->setText(ui_->number_of_particles->text());
     connect(box, &QLineEdit::textChanged, this, &OptimizeTool::update_run_button);
@@ -482,7 +528,7 @@ void OptimizeTool::setup_domain_boxes() {
   } else {
     auto domain_names = session_->get_project()->get_domain_names();
     for (int i = 0; i < domain_names.size(); i++) {
-      auto label = new QLabel(QString::fromStdString(domain_names[i]), this);
+      auto label = new ElidedLabel(QString::fromStdString(domain_names[i]), this);
       label->setAlignment(Qt::AlignRight | Qt::AlignVCenter);
       grid->addWidget(label, i, 1);
       domain_grid_widgets_.push_back(label);