Add correspondence quality panel into Studio analysis

Applications/shapeworks/Commands.cpp

@@ -398,8 +398,8 @@ bool CorrespondenceQualityCommand::execute(const optparse::Values& options, Shar
         const auto& r = sorted_results[i];
         std::cout << "  " << r.subject << "  (domain " << r.domain << ")"
                   << "  norm_mean=" << r.norm_mean << "  (" << std::setprecision(3) << (r.norm_mean * 100.0) << "%)"
-                  << std::setprecision(6) << "  mean=" << r.mean_dist << "  max=" << r.max_dist
-                  << "  bbox_diag=" << r.bbox_diag << "\n";
+                  << std::setprecision(6) << "  mean=" << r.mean_dist << "  median=" << r.median_dist
+                  << "  max=" << r.max_dist << "  bbox_diag=" << r.bbox_diag << "\n";
       }
     }
 
@@ -414,11 +414,13 @@ bool CorrespondenceQualityCommand::execute(const optparse::Values& options, Shar
         boost::filesystem::current_path(oldBasePath);
         return false;
       }
-      csv << "subject,domain,is_template,mean_dist,max_dist,bbox_diag,norm_mean,norm_max\n";
+      csv << "subject,domain,is_template,mean_dist,median_dist,p99_dist,max_dist,bbox_diag,norm_mean,norm_median,"
+             "norm_p99,norm_max\n";
       csv << std::fixed << std::setprecision(8);
       for (const auto& r : report.rows) {
         csv << r.subject << "," << r.domain << "," << (r.is_template ? 1 : 0) << "," << r.mean_dist << ","
-            << r.max_dist << "," << r.bbox_diag << "," << r.norm_mean << "," << r.norm_max << "\n";
+            << r.median_dist << "," << r.p99_dist << "," << r.max_dist << "," << r.bbox_diag << "," << r.norm_mean
+            << "," << r.norm_median << "," << r.norm_p99 << "," << r.norm_max << "\n";
       }
       SW_LOG("Wrote per-subject CSV: {}", out_path.string());
     }

Libs/Particles/CorrespondenceEvaluation.cpp

@@ -43,7 +43,10 @@ Eigen::MatrixXd load_particles_matrix(const std::string& filename) {
   return m;
 }
 
-CorrespondenceQualityStats summarize(std::vector<double> values) {
+}  // namespace
+
+//---------------------------------------------------------------------------
+CorrespondenceQualityStats CorrespondenceEvaluation::summarize(std::vector<double> values) {
   CorrespondenceQualityStats s;
   if (values.empty()) return s;
   std::sort(values.begin(), values.end());
@@ -55,8 +58,85 @@ CorrespondenceQualityStats summarize(std::vector<double> values) {
   return s;
 }
 
-}  // namespace
+//---------------------------------------------------------------------------
+CorrespondenceQualityRow CorrespondenceEvaluation::evaluate_reconstruction(
+    vtkSmartPointer<vtkPolyData> reconstructed, const Mesh& groomed, DistanceMethod method,
+    vtkSmartPointer<vtkDataArray>* out_distance) {
+  CorrespondenceQualityRow row;
+  if (!reconstructed || reconstructed->GetNumberOfPoints() == 0) {
+    return row;
+  }
+
+  const Mesh::DistanceMethod distance_method =
+      (method == DistanceMethod::PointToPoint) ? Mesh::DistanceMethod::PointToPoint : Mesh::DistanceMethod::PointToCell;
+
+  Mesh recon_mesh(reconstructed);
+  auto field = recon_mesh.distance(groomed, distance_method)[0];
+
+  const int n = field->GetNumberOfTuples();
+  if (n == 0) {
+    return row;
+  }
+
+  std::vector<double> values(n);
+  double sum = 0.0;
+  double maxv = 0.0;
+  for (int k = 0; k < n; ++k) {
+    const double v = std::fabs(field->GetTuple1(k));
+    values[k] = v;
+    sum += v;
+    if (v > maxv) maxv = v;
+  }
+
+  const size_t mid = values.size() / 2;
+  std::nth_element(values.begin(), values.begin() + mid, values.end());
+  const double median = values[mid];
+
+  const size_t p99_idx = std::min(values.size() - 1, static_cast<size_t>(0.99 * values.size()));
+  std::nth_element(values.begin(), values.begin() + p99_idx, values.end());
+
+  row.mean_dist = sum / n;
+  row.median_dist = median;
+  row.p99_dist = values[p99_idx];
+  row.max_dist = maxv;
+
+  const auto bbox = groomed.boundingBox();
+  row.bbox_diag = (bbox.max - bbox.min).GetNorm();
+  if (row.bbox_diag > 0.0) {
+    row.norm_mean = row.mean_dist / row.bbox_diag;
+    row.norm_median = row.median_dist / row.bbox_diag;
+    row.norm_p99 = row.p99_dist / row.bbox_diag;
+    row.norm_max = row.max_dist / row.bbox_diag;
+  }
+
+  if (out_distance) {
+    field->SetName("distance");
+    *out_distance = field;
+  }
+
+  return row;
+}
+
+//---------------------------------------------------------------------------
+void CorrespondenceEvaluation::compute_aggregates(CorrespondenceQualityReport& report) {
+  std::vector<double> means;
+  std::vector<double> norm_means;
+  int num_template_rows = 0;
+  for (const auto& r : report.rows) {
+    if (r.is_template) {
+      num_template_rows++;
+      continue;
+    }
+    means.push_back(r.mean_dist);
+    norm_means.push_back(r.norm_mean);
+  }
+  report.num_template_rows = num_template_rows;
+  report.num_evaluated = static_cast<int>(means.size());
+  report.agg_raw = summarize(means);
+  report.agg_norm = summarize(norm_means);
+}
 
+//---------------------------------------------------------------------------
 CorrespondenceQualityReport CorrespondenceEvaluation::evaluate(ProjectHandle project, DistanceMethod method,
                                                                const std::string& output_meshes_dir) {
   if (!project) {
@@ -72,9 +152,6 @@ CorrespondenceQualityReport CorrespondenceEvaluation::evaluate(ProjectHandle pro
     throw std::runtime_error("project has no domains");
   }
 
-  const Mesh::DistanceMethod distance_method =
-      (method == DistanceMethod::PointToPoint) ? Mesh::DistanceMethod::PointToPoint : Mesh::DistanceMethod::PointToCell;
-
   // Pass 1: load particles + groomed paths per (subject, domain). Only keep subjects
   // with complete data across all domains, so the L1-medoid is computed over a
   // consistent cohort and the same global template applies to every domain.
@@ -166,9 +243,6 @@ CorrespondenceQualityReport CorrespondenceEvaluation::evaluate(ProjectHandle pro
 
   CorrespondenceQualityReport report;
   report.template_subject = name_per_subject[template_idx];
-  std::vector<double> all_means;       // pooled raw mean distances (template excluded)
-  std::vector<double> all_norm_means;  // pooled bbox-normalized values (template excluded)
-
   // Pass 2: per-domain warp + distance using the single global template.
   for (int domain = 0; domain < num_domains; ++domain) {
     SW_LOG("=== Domain {} ===", domain);
@@ -189,44 +263,18 @@ CorrespondenceQualityReport CorrespondenceEvaluation::evaluate(ProjectHandle pro
         continue;
       }
 
-      Mesh recon_mesh(reconstructed);
       Mesh groomed_mesh = load_groomed_as_mesh(groomed_per_subject_domain[i][domain]);
-      auto field = recon_mesh.distance(groomed_mesh, distance_method)[0];
-
-      const int n = field->GetNumberOfTuples();
-      if (n == 0) continue;
-      double sum = 0.0;
-      double maxv = 0.0;
-      for (int k = 0; k < n; ++k) {
-        const double v = std::fabs(field->GetTuple1(k));
-        sum += v;
-        if (v > maxv) maxv = v;
-      }
-      const double mean_d = sum / n;
+      vtkSmartPointer<vtkDataArray> field;
+      CorrespondenceQualityRow row = evaluate_reconstruction(reconstructed, groomed_mesh, method, &field);
+      if (!field) continue;
 
-      const auto bbox = groomed_mesh.boundingBox();
-      const double bbox_diag = (bbox.max - bbox.min).GetNorm();
-      const double norm_mean = (bbox_diag > 0.0) ? mean_d / bbox_diag : 0.0;
-      const double norm_max = (bbox_diag > 0.0) ? maxv / bbox_diag : 0.0;
-
-      CorrespondenceQualityRow row;
       row.subject = name_per_subject[i];
       row.domain = domain;
-      row.mean_dist = mean_d;
-      row.max_dist = maxv;
-      row.bbox_diag = bbox_diag;
-      row.norm_mean = norm_mean;
-      row.norm_max = norm_max;
       row.is_template = (i == template_idx);
       report.rows.push_back(row);
 
-      if (!row.is_template) {
-        all_means.push_back(mean_d);
-        all_norm_means.push_back(norm_mean);
-      }
-
       if (!meshes_dir.empty()) {
-        field->SetName("distance");
+        Mesh recon_mesh(reconstructed);
         recon_mesh.setField("distance", field, Mesh::FieldType::Point);
         std::string fname = name_per_subject[i] + "_domain" + std::to_string(domain) + "_reconstructed.vtk";
         if (row.is_template) {
@@ -242,10 +290,7 @@ CorrespondenceQualityReport CorrespondenceEvaluation::evaluate(ProjectHandle pro
     throw std::runtime_error("no subjects evaluated");
   }
 
-  report.num_template_rows = static_cast<int>(report.rows.size() - all_means.size());
-  report.num_evaluated = static_cast<int>(all_means.size());
-  report.agg_raw = summarize(all_means);
-  report.agg_norm = summarize(all_norm_means);
+  compute_aggregates(report);
   return report;
 }
 

Libs/Particles/CorrespondenceEvaluation.h

@@ -1,12 +1,18 @@
 #pragma once
 
+#include <vtkSmartPointer.h>
+
 #include <Eigen/Core>
 #include <memory>
 #include <string>
 #include <vector>
 
+class vtkDataArray;
+class vtkPolyData;
+
 namespace shapeworks {
 
+class Mesh;
 class Project;
 using ProjectHandle = std::shared_ptr<Project>;
 
@@ -15,9 +21,13 @@ struct CorrespondenceQualityRow {
   std::string subject;
   int domain = 0;
   double mean_dist = 0.0;    //!< mean point-to-cell (or point-to-point) distance, reconstructed -> groomed
+  double median_dist = 0.0;  //!< median per-vertex distance
+  double p99_dist = 0.0;     //!< 99th percentile per-vertex distance, a max that ignores single outlier vertices
   double max_dist = 0.0;     //!< max per-vertex distance
   double bbox_diag = 0.0;    //!< diagonal of the subject's groomed-mesh bounding box
   double norm_mean = 0.0;    //!< mean_dist / bbox_diag (scale-invariant)
+  double norm_median = 0.0;  //!< median_dist / bbox_diag
+  double norm_p99 = 0.0;     //!< p99_dist / bbox_diag
   double norm_max = 0.0;     //!< max_dist / bbox_diag
   bool is_template = false;  //!< true for the L1-medoid template row (excluded from aggregates)
 };
@@ -54,6 +64,11 @@ struct CorrespondenceQualityReport {
  * The template row itself is included in `rows` (with is_template=true) but
  * excluded from aggregate statistics — its reconstruction is near-identity
  * and would skew small cohorts.
+ *
+ * `evaluate()` drives the whole thing from a project file. Callers that already
+ * have reconstructions in memory (Studio, which reconstructs through its own
+ * configured mesh warper) should use `evaluate_reconstruction()` and
+ * `compute_aggregates()` instead so the metric definition stays in one place.
  */
 class CorrespondenceEvaluation {
  public:
@@ -73,6 +88,24 @@ class CorrespondenceEvaluation {
   static CorrespondenceQualityReport evaluate(ProjectHandle project,
                                               DistanceMethod method = DistanceMethod::PointToCell,
                                               const std::string& output_meshes_dir = "");
+
+  //! Score a single already-reconstructed mesh against its groomed target.
+  //!
+  //! Fills everything on the row except `subject`, `domain` and `is_template`,
+  //! which the caller owns. If \p out_distance is non-null it receives the
+  //! per-vertex distance field (named "distance"), for surface display or
+  //! writing alongside the mesh.
+  //!
+  //! Returns a default-constructed row if \p reconstructed is null or empty.
+  static CorrespondenceQualityRow evaluate_reconstruction(vtkSmartPointer<vtkPolyData> reconstructed,
+                                                          const Mesh& groomed, DistanceMethod method,
+                                                          vtkSmartPointer<vtkDataArray>* out_distance = nullptr);
+
+  //! Summary statistics (mean/median/p95/max) over a set of values.
+  static CorrespondenceQualityStats summarize(std::vector<double> values);
+
+  //! Fill num_evaluated, num_template_rows, agg_raw and agg_norm from report.rows.
+  static void compute_aggregates(CorrespondenceQualityReport& report);
 };
 
 }  // namespace shapeworks

Libs/Python/ShapeworksPython.cpp

@@ -1401,9 +1401,13 @@ PYBIND11_MODULE(shapeworks_py, m) {
       .def_readonly("subject", &CorrespondenceQualityRow::subject)
       .def_readonly("domain", &CorrespondenceQualityRow::domain)
       .def_readonly("mean_dist", &CorrespondenceQualityRow::mean_dist)
+      .def_readonly("median_dist", &CorrespondenceQualityRow::median_dist)
+      .def_readonly("p99_dist", &CorrespondenceQualityRow::p99_dist)
       .def_readonly("max_dist", &CorrespondenceQualityRow::max_dist)
       .def_readonly("bbox_diag", &CorrespondenceQualityRow::bbox_diag)
       .def_readonly("norm_mean", &CorrespondenceQualityRow::norm_mean)
+      .def_readonly("norm_median", &CorrespondenceQualityRow::norm_median)
+      .def_readonly("norm_p99", &CorrespondenceQualityRow::norm_p99)
       .def_readonly("norm_max", &CorrespondenceQualityRow::norm_max)
       .def_readonly("is_template", &CorrespondenceQualityRow::is_template);