Fix shared boundary donut and contour correspondence issues

- MeshUtils: replace fragile is_clockwise atan2 test with summed signed-area
  vector (uses all vertices, immune to branch-cut flips); rotate boundary
  loops to start at canonical max-Y vertex for consistent ordering across
  subjects.
- ContourDomain: return length² from GetSurfaceArea so contours participate
  in sampling-scale auto-scaling at a magnitude comparable to meshes.
- SamplingFunction: apply auto-scale uniformly; contours no longer require
  special-case skipping.
- GroomParameters: skip phantom shared-boundary entry on migration when
  domain names are empty.
- Regenerate shared boundary test ground truth for canonical start vertex.

Author: akenmorris

Libs/Groom/GroomParameters.cpp

@@ -589,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/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();

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;
 

Testing/data/shared_boundary/00_out_c.vtp

@@ -1,3 +1,3 @@
 version https://git-lfs.github.com/spec/v1
-oid sha256:06fb0a7bfbf2fe03d29f855a7156adacd2bd46596f26a36409df0531995daff6
-size 2851
+oid sha256:95657cfe31442ba0bdeaeb7116752883c8dae291077d5e0c5cf95cee52c41a23
+size 4074