From 784c7b20546e0aca73af8d916be67ccc42932296 Mon Sep 17 00:00:00 2001 From: Alan Morris Date: Thu, 23 Jul 2026 12:36:49 -0600 Subject: [PATCH 1/4] Add correspondence-quality command New `shapeworks correspondence-quality` subcommand evaluates per-subject correspondence by reconstructing each shape from its local particles via biharmonic mesh warp from the cohort L1-medoid template (matching Studio's median selection) and measuring distance to the groomed mesh. - Per-subject mean/max point-to-cell (or point-to-point) distance, plus bbox-diagonal-normalized values for scale-invariant good/bad judgement. - Aggregates (mean/median/p95/max) exclude the template row so small cohorts aren't skewed by its near-identity reconstruction. - Handles both mesh and distance-transform groomed inputs. - Optional per-subject CSV and --output_meshes dump with embedded per-vertex distance field for visual inspection. --- Applications/shapeworks/Commands.cpp | 381 +++++++++++++++++++++++++ Applications/shapeworks/Commands.h | 1 + Applications/shapeworks/shapeworks.cpp | 1 + 3 files changed, 383 insertions(+) diff --git a/Applications/shapeworks/Commands.cpp b/Applications/shapeworks/Commands.cpp index e196b400a30..75e69c85363 100644 --- a/Applications/shapeworks/Commands.cpp +++ b/Applications/shapeworks/Commands.cpp @@ -5,13 +5,24 @@ #include #include #include +#include +#include +#include #include #include #include +#include #include +#include #include #include +#include +#include +#include +#include +#include + #include #include @@ -276,6 +287,376 @@ bool AnalyzeCommand::execute(const optparse::Values& options, SharedCommandData& } } +/////////////////////////////////////////////////////////////////////////////// +// CorrespondenceQuality +/////////////////////////////////////////////////////////////////////////////// +namespace { + +Mesh load_groomed_as_mesh(const std::string& path) { + const std::string ext = StringUtils::getLowerExtension(path); + if (ext == ".nrrd" || ext == ".mha" || ext == ".mhd" || ext == ".nii" || ext == ".gz" || ext == ".tif" || + ext == ".tiff") { + Image img(path); + return img.toMesh(0.0); + } + return Mesh(path); +} + +Eigen::MatrixXd load_particles_matrix(const std::string& filename) { + Eigen::VectorXd v = particles::read_particles(filename); + const int n = static_cast(v.size() / 3); + // read_particles returns interleaved [x0,y0,z0, x1,y1,z1, ...] + Eigen::MatrixXd m(n, 3); + for (int i = 0; i < n; ++i) { + m(i, 0) = v(3 * i + 0); + m(i, 1) = v(3 * i + 1); + m(i, 2) = v(3 * i + 2); + } + return m; +} + +struct DistanceStats { + double mean = 0.0; + double max = 0.0; + double median = 0.0; + double p95 = 0.0; +}; + +DistanceStats summarize(std::vector values) { + DistanceStats s; + if (values.empty()) return s; + std::sort(values.begin(), values.end()); + s.mean = std::accumulate(values.begin(), values.end(), 0.0) / values.size(); + s.max = values.back(); + s.median = values[values.size() / 2]; + const size_t p95_idx = std::min(values.size() - 1, static_cast(0.95 * values.size())); + s.p95 = values[p95_idx]; + return s; +} + +struct SubjectResult { + std::string subject; + int domain; + double mean_dist; + double max_dist; + double bbox_diag; + double norm_mean; // mean_dist / bbox_diag + double norm_max; // max_dist / bbox_diag + bool is_template; +}; + +} // namespace + +void CorrespondenceQualityCommand::buildParser() { + const std::string prog = "correspondence-quality"; + const std::string desc = + "Evaluate per-subject correspondence quality by reconstructing each shape from its local particles " + "(biharmonic mesh warp from the cohort L1-medoid template, matching Studio's median selection) and " + "measuring distance to the groomed mesh. Reports per-subject and aggregate statistics."; + parser.prog(prog).description(desc); + + parser.add_option("--name").action("store").type("string").set_default("").help( + "Path to project file (.swproj or .xlsx)."); + parser.add_option("--output").action("store").type("string").set_default("").help( + "Optional path to write per-subject CSV."); + parser.add_option("--output_meshes").action("store").type("string").set_default("").help( + "Optional directory; write reconstructed meshes with per-vertex distance field for visual inspection."); + std::list methods{"point-to-cell", "point-to-point"}; + parser.add_option("--method") + .action("store") + .type("choice") + .choices(methods.begin(), methods.end()) + .set_default("point-to-cell") + .help("Distance method [default: %default]."); + parser.add_option("--worst") + .action("store") + .type("int") + .set_default(5) + .help("Number of worst-ranked subjects to print [default: %default]."); + + Command::buildParser(); +} + +bool CorrespondenceQualityCommand::execute(const optparse::Values& options, SharedCommandData& sharedData) { + const std::string& projectFile(static_cast(options.get("name"))); + const std::string& outputFile(static_cast(options.get("output"))); + const std::string& outputMeshesDir(static_cast(options.get("output_meshes"))); + const std::string methodopt(options.get("method")); + const int worst_n = static_cast(options.get("worst")); + + if (projectFile.empty()) { + std::cerr << "Must specify project name with --name \n"; + return false; + } + + const Mesh::DistanceMethod distance_method = + (methodopt == "point-to-point") ? Mesh::DistanceMethod::PointToPoint : Mesh::DistanceMethod::PointToCell; + + try { + ProjectHandle project = std::make_shared(); + + const auto oldBasePath = boost::filesystem::current_path(); + auto base = StringUtils::getPath(projectFile); + auto filename = StringUtils::getFilename(projectFile); + if (base != projectFile) { + SW_LOG("Setting path to {}", base); + boost::filesystem::current_path(base.c_str()); + project->set_filename(filename); + } + + SW_LOG("Loading project: {}", filename); + project->load(filename); + + auto subjects = project->get_non_excluded_subjects(); + if (subjects.empty()) { + SW_ERROR("No (non-excluded) subjects in project"); + boost::filesystem::current_path(oldBasePath); + return false; + } + const int num_domains = project->get_number_of_domains_per_subject(); + if (num_domains <= 0) { + SW_ERROR("Project has no domains"); + boost::filesystem::current_path(oldBasePath); + return false; + } + + // Pass 1: load every (subject, domain) particles + groomed paths. Keep only subjects + // that have complete data across all domains, so the L1-medoid is computed over a + // consistent cohort and the same global template is used in every domain. + std::vector name_per_subject; + std::vector> particles_per_subject_domain; // [subject][domain] + std::vector> groomed_per_subject_domain; // [subject][domain] + + for (auto& subj : subjects) { + auto particle_files = subj->get_local_particle_filenames(); + auto groomed_files = subj->get_groomed_filenames(); + + bool complete = true; + std::vector subj_particles; + std::vector subj_groomed; + subj_particles.reserve(num_domains); + subj_groomed.reserve(num_domains); + + for (int d = 0; d < num_domains; ++d) { + if (d >= static_cast(particle_files.size()) || particle_files[d].empty() || + d >= static_cast(groomed_files.size()) || groomed_files[d].empty()) { + SW_LOG("Skipping subject '{}': missing particles or groomed for domain {}", subj->get_display_name(), d); + complete = false; + break; + } + try { + subj_particles.push_back(load_particles_matrix(particle_files[d])); + } catch (const std::exception& e) { + SW_LOG("Skipping subject '{}' domain {}: particle load failed: {}", subj->get_display_name(), d, e.what()); + complete = false; + break; + } + subj_groomed.push_back(groomed_files[d]); + } + + if (!complete) continue; + name_per_subject.push_back(subj->get_display_name()); + particles_per_subject_domain.push_back(std::move(subj_particles)); + groomed_per_subject_domain.push_back(std::move(subj_groomed)); + } + + const int num_subjects = static_cast(name_per_subject.size()); + if (num_subjects < 2) { + SW_ERROR("Need at least 2 subjects with complete data across all domains (have {})", num_subjects); + boost::filesystem::current_path(oldBasePath); + return false; + } + + // Verify particle counts match across subjects per domain + std::vector num_particles_per_domain(num_domains); + for (int d = 0; d < num_domains; ++d) { + num_particles_per_domain[d] = static_cast(particles_per_subject_domain[0][d].rows()); + for (int s = 1; s < num_subjects; ++s) { + if (particles_per_subject_domain[s][d].rows() != num_particles_per_domain[d]) { + SW_ERROR("Domain {}: subject '{}' has {} particles, expected {}", d, name_per_subject[s], + particles_per_subject_domain[s][d].rows(), num_particles_per_domain[d]); + boost::filesystem::current_path(oldBasePath); + return false; + } + } + } + + // L1-medoid template selection over concatenated per-domain local particles (matches + // Studio's compute_median_shape: argmin_i sum_j ||x_i - x_j||_1). + int template_idx = 0; + double best_l1_sum = std::numeric_limits::infinity(); + for (int i = 0; i < num_subjects; ++i) { + double sum_l1 = 0.0; + for (int j = 0; j < num_subjects; ++j) { + if (i == j) continue; + double pair_l1 = 0.0; + for (int d = 0; d < num_domains; ++d) { + pair_l1 += (particles_per_subject_domain[i][d] - particles_per_subject_domain[j][d]).cwiseAbs().sum(); + } + sum_l1 += pair_l1; + } + if (sum_l1 < best_l1_sum) { + best_l1_sum = sum_l1; + template_idx = i; + } + } + SW_LOG("Template subject (L1-medoid): '{}' (sum of L1 distances to others = {:.4f})", + name_per_subject[template_idx], best_l1_sum); + + // Resolve output meshes dir (create if missing) + boost::filesystem::path meshes_dir; + if (!outputMeshesDir.empty()) { + meshes_dir = boost::filesystem::path(outputMeshesDir); + if (!meshes_dir.is_absolute()) { + meshes_dir = boost::filesystem::path(oldBasePath) / meshes_dir; + } + boost::filesystem::create_directories(meshes_dir); + SW_LOG("Writing reconstructed meshes to: {}", meshes_dir.string()); + } + + std::vector all_results; + std::vector all_means; // pooled raw mean distances (excluding template) + std::vector all_norm_means; // pooled bbox-normalized mean distances (excluding template) + + // Pass 2: per-domain warp + distance using the single global template. + for (int domain = 0; domain < num_domains; ++domain) { + SW_LOG("=== Domain {} ===", domain); + + Mesh template_mesh = load_groomed_as_mesh(groomed_per_subject_domain[template_idx][domain]); + MeshWarper warper; + warper.set_warp_method(WarpMethod::Biharmonic); + warper.set_reference_mesh(template_mesh.getVTKMesh(), particles_per_subject_domain[template_idx][domain]); + if (!warper.generate_warp()) { + SW_ERROR("Domain {}: failed to generate warp from template '{}'", domain, name_per_subject[template_idx]); + boost::filesystem::current_path(oldBasePath); + return false; + } + + for (int i = 0; i < num_subjects; ++i) { + vtkSmartPointer reconstructed = warper.build_mesh(particles_per_subject_domain[i][domain]); + if (!reconstructed) { + SW_LOG("Domain {}: subject '{}' reconstruction returned null, skipping", domain, name_per_subject[i]); + 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; + + 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; + + const bool is_template = (i == template_idx); + all_results.push_back({name_per_subject[i], domain, mean_d, maxv, bbox_diag, norm_mean, norm_max, is_template}); + if (!is_template) { + all_means.push_back(mean_d); + all_norm_means.push_back(norm_mean); + } + + // Optional: write reconstructed mesh with distance field + if (!meshes_dir.empty()) { + field->SetName("distance"); + recon_mesh.setField("distance", field, Mesh::FieldType::Point); + std::string fname = name_per_subject[i] + "_domain" + std::to_string(domain) + "_reconstructed.vtk"; + if (is_template) fname = name_per_subject[i] + "_domain" + std::to_string(domain) + "_reconstructed_TEMPLATE.vtk"; + boost::filesystem::path out_mesh = meshes_dir / fname; + recon_mesh.write(out_mesh.string()); + } + } + } + + if (all_results.empty()) { + SW_ERROR("No subjects evaluated"); + boost::filesystem::current_path(oldBasePath); + return false; + } + + // Aggregate (excluding template — its reconstruction is near-identity and would skew + // small cohorts; e.g. with N=4 it dominates 25% of the average). + DistanceStats agg = summarize(all_means); + DistanceStats agg_norm = summarize(all_norm_means); + const size_t num_template_rows = all_results.size() - all_means.size(); + + std::cout << "\n=== Correspondence Quality Summary ===\n"; + std::cout << "Subjects evaluated: " << all_results.size() << " (" << all_means.size() << " in aggregate, " + << num_template_rows << " template row(s) excluded)\n"; + std::cout << "Template subject: " << name_per_subject[template_idx] << "\n"; + std::cout << "Distance method: " << methodopt << "\n"; + std::cout << std::fixed << std::setprecision(6); + std::cout << "Per-subject mean distance (reconstructed -> groomed):\n"; + std::cout << " mean = " << agg.mean << "\n"; + std::cout << " median = " << agg.median << "\n"; + std::cout << " p95 = " << agg.p95 << "\n"; + std::cout << " max = " << agg.max << "\n"; + std::cout << "Normalized by per-subject groomed-mesh bbox diagonal (fraction):\n"; + std::cout << " mean = " << agg_norm.mean << " (" << std::setprecision(3) << (agg_norm.mean * 100.0) << "%)\n"; + std::cout << std::setprecision(6); + std::cout << " median = " << agg_norm.median << "\n"; + std::cout << " p95 = " << agg_norm.p95 << "\n"; + std::cout << " max = " << agg_norm.max << "\n"; + + // Worst-N (sorted by normalized mean for scale-invariant ranking; template excluded) + if (worst_n > 0) { + std::vector sorted_results; + sorted_results.reserve(all_results.size()); + for (const auto& r : all_results) { + if (!r.is_template) sorted_results.push_back(r); + } + std::sort(sorted_results.begin(), sorted_results.end(), + [](const SubjectResult& a, const SubjectResult& b) { return a.norm_mean > b.norm_mean; }); + const int limit = std::min(worst_n, static_cast(sorted_results.size())); + std::cout << "\nWorst " << limit << " subjects (ranked by normalized mean; template excluded):\n"; + for (int i = 0; i < limit; ++i) { + 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"; + } + } + + // CSV + if (!outputFile.empty()) { + boost::filesystem::path out_path(outputFile); + if (!out_path.is_absolute()) { + out_path = boost::filesystem::path(oldBasePath) / out_path; + } + std::ofstream csv(out_path.string()); + if (!csv) { + SW_ERROR("Failed to open output file: {}", out_path.string()); + boost::filesystem::current_path(oldBasePath); + return false; + } + csv << "subject,domain,is_template,mean_dist,max_dist,bbox_diag,norm_mean,norm_max\n"; + csv << std::fixed << std::setprecision(8); + for (const auto& r : all_results) { + 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"; + } + SW_LOG("Wrote per-subject CSV: {}", out_path.string()); + } + + boost::filesystem::current_path(oldBasePath); + return true; + } catch (std::exception& e) { + std::cerr << "Error: " << e.what() << "\n"; + return false; + } +} + /////////////////////////////////////////////////////////////////////////////// // Convert Project /////////////////////////////////////////////////////////////////////////////// diff --git a/Applications/shapeworks/Commands.h b/Applications/shapeworks/Commands.h index 77caf86c285..764eab54dee 100644 --- a/Applications/shapeworks/Commands.h +++ b/Applications/shapeworks/Commands.h @@ -101,6 +101,7 @@ COMMAND_DECLARE(Seed, ShapeworksCommand); COMMAND_DECLARE(OptimizeCommand, OptimizeCommandGroup); COMMAND_DECLARE(GroomCommand, GroomCommandGroup); COMMAND_DECLARE(AnalyzeCommand, AnalyzeCommandGroup); +COMMAND_DECLARE(CorrespondenceQualityCommand, AnalyzeCommandGroup); COMMAND_DECLARE(ConvertProjectCommand, ProjectCommandGroup); COMMAND_DECLARE(DeepSSMCommand, DeepSSMCommandGroup); diff --git a/Applications/shapeworks/shapeworks.cpp b/Applications/shapeworks/shapeworks.cpp index a77516a39ec..a7300a99d15 100644 --- a/Applications/shapeworks/shapeworks.cpp +++ b/Applications/shapeworks/shapeworks.cpp @@ -110,6 +110,7 @@ int main(int argc, char *argv[]) shapeworks.addCommand(OptimizeCommand::getCommand()); shapeworks.addCommand(GroomCommand::getCommand()); shapeworks.addCommand(AnalyzeCommand::getCommand()); + shapeworks.addCommand(CorrespondenceQualityCommand::getCommand()); shapeworks.addCommand(ConvertProjectCommand::getCommand()); shapeworks.addCommand(DeepSSMCommand::getCommand()); From 1369b5525438355df36bc8efcb75c35fe52e1248 Mon Sep 17 00:00:00 2001 From: Alan Morris Date: Thu, 23 Jul 2026 12:31:31 -0600 Subject: [PATCH 2/4] Extract CorrespondenceEvaluation and expose via Python Move the correspondence-quality algorithm out of the shapeworks CLI handler into a reusable class in Libs/Particles, alongside ShapeEvaluation, and add pybind11 bindings so it can be called from Python without shelling out. - New Libs/Particles/CorrespondenceEvaluation.{h,cpp} with the CorrespondenceEvaluation class and CorrespondenceQualityReport / Row / Stats result structs. - shapeworks correspondence-quality command is now a thin CLI wrapper around the class; output is byte-identical to before. - Python: sw.CorrespondenceEvaluation.evaluate(project, method, output_meshes_dir) returns a CorrespondenceQualityReport with per- subject rows and aggregate stats (template row excluded). - Libs/Particles now depends on Project. --- Applications/shapeworks/Commands.cpp | 292 +++----------------- Libs/Particles/CMakeLists.txt | 3 + Libs/Particles/CorrespondenceEvaluation.cpp | 252 +++++++++++++++++ Libs/Particles/CorrespondenceEvaluation.h | 78 ++++++ Libs/Python/ShapeworksPython.cpp | 38 +++ 5 files changed, 404 insertions(+), 259 deletions(-) create mode 100644 Libs/Particles/CorrespondenceEvaluation.cpp create mode 100644 Libs/Particles/CorrespondenceEvaluation.h diff --git a/Applications/shapeworks/Commands.cpp b/Applications/shapeworks/Commands.cpp index 75e69c85363..8591862fc02 100644 --- a/Applications/shapeworks/Commands.cpp +++ b/Applications/shapeworks/Commands.cpp @@ -5,13 +5,10 @@ #include #include #include -#include -#include -#include #include #include #include -#include +#include #include #include #include @@ -20,8 +17,6 @@ #include #include #include -#include -#include #include #include @@ -290,63 +285,6 @@ bool AnalyzeCommand::execute(const optparse::Values& options, SharedCommandData& /////////////////////////////////////////////////////////////////////////////// // CorrespondenceQuality /////////////////////////////////////////////////////////////////////////////// -namespace { - -Mesh load_groomed_as_mesh(const std::string& path) { - const std::string ext = StringUtils::getLowerExtension(path); - if (ext == ".nrrd" || ext == ".mha" || ext == ".mhd" || ext == ".nii" || ext == ".gz" || ext == ".tif" || - ext == ".tiff") { - Image img(path); - return img.toMesh(0.0); - } - return Mesh(path); -} - -Eigen::MatrixXd load_particles_matrix(const std::string& filename) { - Eigen::VectorXd v = particles::read_particles(filename); - const int n = static_cast(v.size() / 3); - // read_particles returns interleaved [x0,y0,z0, x1,y1,z1, ...] - Eigen::MatrixXd m(n, 3); - for (int i = 0; i < n; ++i) { - m(i, 0) = v(3 * i + 0); - m(i, 1) = v(3 * i + 1); - m(i, 2) = v(3 * i + 2); - } - return m; -} - -struct DistanceStats { - double mean = 0.0; - double max = 0.0; - double median = 0.0; - double p95 = 0.0; -}; - -DistanceStats summarize(std::vector values) { - DistanceStats s; - if (values.empty()) return s; - std::sort(values.begin(), values.end()); - s.mean = std::accumulate(values.begin(), values.end(), 0.0) / values.size(); - s.max = values.back(); - s.median = values[values.size() / 2]; - const size_t p95_idx = std::min(values.size() - 1, static_cast(0.95 * values.size())); - s.p95 = values[p95_idx]; - return s; -} - -struct SubjectResult { - std::string subject; - int domain; - double mean_dist; - double max_dist; - double bbox_diag; - double norm_mean; // mean_dist / bbox_diag - double norm_max; // max_dist / bbox_diag - bool is_template; -}; - -} // namespace - void CorrespondenceQualityCommand::buildParser() { const std::string prog = "correspondence-quality"; const std::string desc = @@ -389,8 +327,8 @@ bool CorrespondenceQualityCommand::execute(const optparse::Values& options, Shar return false; } - const Mesh::DistanceMethod distance_method = - (methodopt == "point-to-point") ? Mesh::DistanceMethod::PointToPoint : Mesh::DistanceMethod::PointToCell; + const auto method = (methodopt == "point-to-point") ? CorrespondenceEvaluation::DistanceMethod::PointToPoint + : CorrespondenceEvaluation::DistanceMethod::PointToCell; try { ProjectHandle project = std::make_shared(); @@ -407,216 +345,53 @@ bool CorrespondenceQualityCommand::execute(const optparse::Values& options, Shar SW_LOG("Loading project: {}", filename); project->load(filename); - auto subjects = project->get_non_excluded_subjects(); - if (subjects.empty()) { - SW_ERROR("No (non-excluded) subjects in project"); - boost::filesystem::current_path(oldBasePath); - return false; - } - const int num_domains = project->get_number_of_domains_per_subject(); - if (num_domains <= 0) { - SW_ERROR("Project has no domains"); - boost::filesystem::current_path(oldBasePath); - return false; - } - - // Pass 1: load every (subject, domain) particles + groomed paths. Keep only subjects - // that have complete data across all domains, so the L1-medoid is computed over a - // consistent cohort and the same global template is used in every domain. - std::vector name_per_subject; - std::vector> particles_per_subject_domain; // [subject][domain] - std::vector> groomed_per_subject_domain; // [subject][domain] - - for (auto& subj : subjects) { - auto particle_files = subj->get_local_particle_filenames(); - auto groomed_files = subj->get_groomed_filenames(); - - bool complete = true; - std::vector subj_particles; - std::vector subj_groomed; - subj_particles.reserve(num_domains); - subj_groomed.reserve(num_domains); - - for (int d = 0; d < num_domains; ++d) { - if (d >= static_cast(particle_files.size()) || particle_files[d].empty() || - d >= static_cast(groomed_files.size()) || groomed_files[d].empty()) { - SW_LOG("Skipping subject '{}': missing particles or groomed for domain {}", subj->get_display_name(), d); - complete = false; - break; - } - try { - subj_particles.push_back(load_particles_matrix(particle_files[d])); - } catch (const std::exception& e) { - SW_LOG("Skipping subject '{}' domain {}: particle load failed: {}", subj->get_display_name(), d, e.what()); - complete = false; - break; - } - subj_groomed.push_back(groomed_files[d]); - } - - if (!complete) continue; - name_per_subject.push_back(subj->get_display_name()); - particles_per_subject_domain.push_back(std::move(subj_particles)); - groomed_per_subject_domain.push_back(std::move(subj_groomed)); - } - - const int num_subjects = static_cast(name_per_subject.size()); - if (num_subjects < 2) { - SW_ERROR("Need at least 2 subjects with complete data across all domains (have {})", num_subjects); - boost::filesystem::current_path(oldBasePath); - return false; - } - - // Verify particle counts match across subjects per domain - std::vector num_particles_per_domain(num_domains); - for (int d = 0; d < num_domains; ++d) { - num_particles_per_domain[d] = static_cast(particles_per_subject_domain[0][d].rows()); - for (int s = 1; s < num_subjects; ++s) { - if (particles_per_subject_domain[s][d].rows() != num_particles_per_domain[d]) { - SW_ERROR("Domain {}: subject '{}' has {} particles, expected {}", d, name_per_subject[s], - particles_per_subject_domain[s][d].rows(), num_particles_per_domain[d]); - boost::filesystem::current_path(oldBasePath); - return false; - } - } - } - - // L1-medoid template selection over concatenated per-domain local particles (matches - // Studio's compute_median_shape: argmin_i sum_j ||x_i - x_j||_1). - int template_idx = 0; - double best_l1_sum = std::numeric_limits::infinity(); - for (int i = 0; i < num_subjects; ++i) { - double sum_l1 = 0.0; - for (int j = 0; j < num_subjects; ++j) { - if (i == j) continue; - double pair_l1 = 0.0; - for (int d = 0; d < num_domains; ++d) { - pair_l1 += (particles_per_subject_domain[i][d] - particles_per_subject_domain[j][d]).cwiseAbs().sum(); - } - sum_l1 += pair_l1; - } - if (sum_l1 < best_l1_sum) { - best_l1_sum = sum_l1; - template_idx = i; - } - } - SW_LOG("Template subject (L1-medoid): '{}' (sum of L1 distances to others = {:.4f})", - name_per_subject[template_idx], best_l1_sum); - - // Resolve output meshes dir (create if missing) - boost::filesystem::path meshes_dir; + // Resolve output-meshes dir relative to the caller's original CWD. + std::string meshes_dir; if (!outputMeshesDir.empty()) { - meshes_dir = boost::filesystem::path(outputMeshesDir); - if (!meshes_dir.is_absolute()) { - meshes_dir = boost::filesystem::path(oldBasePath) / meshes_dir; - } - boost::filesystem::create_directories(meshes_dir); - SW_LOG("Writing reconstructed meshes to: {}", meshes_dir.string()); - } - - std::vector all_results; - std::vector all_means; // pooled raw mean distances (excluding template) - std::vector all_norm_means; // pooled bbox-normalized mean distances (excluding template) - - // Pass 2: per-domain warp + distance using the single global template. - for (int domain = 0; domain < num_domains; ++domain) { - SW_LOG("=== Domain {} ===", domain); - - Mesh template_mesh = load_groomed_as_mesh(groomed_per_subject_domain[template_idx][domain]); - MeshWarper warper; - warper.set_warp_method(WarpMethod::Biharmonic); - warper.set_reference_mesh(template_mesh.getVTKMesh(), particles_per_subject_domain[template_idx][domain]); - if (!warper.generate_warp()) { - SW_ERROR("Domain {}: failed to generate warp from template '{}'", domain, name_per_subject[template_idx]); - boost::filesystem::current_path(oldBasePath); - return false; - } - - for (int i = 0; i < num_subjects; ++i) { - vtkSmartPointer reconstructed = warper.build_mesh(particles_per_subject_domain[i][domain]); - if (!reconstructed) { - SW_LOG("Domain {}: subject '{}' reconstruction returned null, skipping", domain, name_per_subject[i]); - 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; - - 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; - - const bool is_template = (i == template_idx); - all_results.push_back({name_per_subject[i], domain, mean_d, maxv, bbox_diag, norm_mean, norm_max, is_template}); - if (!is_template) { - all_means.push_back(mean_d); - all_norm_means.push_back(norm_mean); - } - - // Optional: write reconstructed mesh with distance field - if (!meshes_dir.empty()) { - field->SetName("distance"); - recon_mesh.setField("distance", field, Mesh::FieldType::Point); - std::string fname = name_per_subject[i] + "_domain" + std::to_string(domain) + "_reconstructed.vtk"; - if (is_template) fname = name_per_subject[i] + "_domain" + std::to_string(domain) + "_reconstructed_TEMPLATE.vtk"; - boost::filesystem::path out_mesh = meshes_dir / fname; - recon_mesh.write(out_mesh.string()); - } - } + boost::filesystem::path mp(outputMeshesDir); + if (!mp.is_absolute()) mp = boost::filesystem::path(oldBasePath) / mp; + meshes_dir = mp.string(); } - if (all_results.empty()) { - SW_ERROR("No subjects evaluated"); + CorrespondenceQualityReport report; + try { + report = CorrespondenceEvaluation::evaluate(project, method, meshes_dir); + } catch (const std::exception& e) { + SW_ERROR("{}", e.what()); boost::filesystem::current_path(oldBasePath); return false; } - // Aggregate (excluding template — its reconstruction is near-identity and would skew - // small cohorts; e.g. with N=4 it dominates 25% of the average). - DistanceStats agg = summarize(all_means); - DistanceStats agg_norm = summarize(all_norm_means); - const size_t num_template_rows = all_results.size() - all_means.size(); - std::cout << "\n=== Correspondence Quality Summary ===\n"; - std::cout << "Subjects evaluated: " << all_results.size() << " (" << all_means.size() << " in aggregate, " - << num_template_rows << " template row(s) excluded)\n"; - std::cout << "Template subject: " << name_per_subject[template_idx] << "\n"; + std::cout << "Subjects evaluated: " << report.rows.size() << " (" << report.num_evaluated << " in aggregate, " + << report.num_template_rows << " template row(s) excluded)\n"; + std::cout << "Template subject: " << report.template_subject << "\n"; std::cout << "Distance method: " << methodopt << "\n"; std::cout << std::fixed << std::setprecision(6); std::cout << "Per-subject mean distance (reconstructed -> groomed):\n"; - std::cout << " mean = " << agg.mean << "\n"; - std::cout << " median = " << agg.median << "\n"; - std::cout << " p95 = " << agg.p95 << "\n"; - std::cout << " max = " << agg.max << "\n"; + std::cout << " mean = " << report.agg_raw.mean << "\n"; + std::cout << " median = " << report.agg_raw.median << "\n"; + std::cout << " p95 = " << report.agg_raw.p95 << "\n"; + std::cout << " max = " << report.agg_raw.max << "\n"; std::cout << "Normalized by per-subject groomed-mesh bbox diagonal (fraction):\n"; - std::cout << " mean = " << agg_norm.mean << " (" << std::setprecision(3) << (agg_norm.mean * 100.0) << "%)\n"; + std::cout << " mean = " << report.agg_norm.mean << " (" << std::setprecision(3) + << (report.agg_norm.mean * 100.0) << "%)\n"; std::cout << std::setprecision(6); - std::cout << " median = " << agg_norm.median << "\n"; - std::cout << " p95 = " << agg_norm.p95 << "\n"; - std::cout << " max = " << agg_norm.max << "\n"; + std::cout << " median = " << report.agg_norm.median << "\n"; + std::cout << " p95 = " << report.agg_norm.p95 << "\n"; + std::cout << " max = " << report.agg_norm.max << "\n"; - // Worst-N (sorted by normalized mean for scale-invariant ranking; template excluded) + // Worst-N (sorted by normalized mean; template excluded) if (worst_n > 0) { - std::vector sorted_results; - sorted_results.reserve(all_results.size()); - for (const auto& r : all_results) { + std::vector sorted_results; + sorted_results.reserve(report.rows.size()); + for (const auto& r : report.rows) { if (!r.is_template) sorted_results.push_back(r); } std::sort(sorted_results.begin(), sorted_results.end(), - [](const SubjectResult& a, const SubjectResult& b) { return a.norm_mean > b.norm_mean; }); + [](const CorrespondenceQualityRow& a, const CorrespondenceQualityRow& b) { + return a.norm_mean > b.norm_mean; + }); const int limit = std::min(worst_n, static_cast(sorted_results.size())); std::cout << "\nWorst " << limit << " subjects (ranked by normalized mean; template excluded):\n"; for (int i = 0; i < limit; ++i) { @@ -628,7 +403,6 @@ bool CorrespondenceQualityCommand::execute(const optparse::Values& options, Shar } } - // CSV if (!outputFile.empty()) { boost::filesystem::path out_path(outputFile); if (!out_path.is_absolute()) { @@ -642,7 +416,7 @@ bool CorrespondenceQualityCommand::execute(const optparse::Values& options, Shar } csv << "subject,domain,is_template,mean_dist,max_dist,bbox_diag,norm_mean,norm_max\n"; csv << std::fixed << std::setprecision(8); - for (const auto& r : all_results) { + 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"; } diff --git a/Libs/Particles/CMakeLists.txt b/Libs/Particles/CMakeLists.txt index dd45f3632c5..63095d5006c 100644 --- a/Libs/Particles/CMakeLists.txt +++ b/Libs/Particles/CMakeLists.txt @@ -4,6 +4,7 @@ set(Particles_sources ParticleSystemEvaluation.cpp ParticleShapeStatistics.cpp ShapeEvaluation.cpp + CorrespondenceEvaluation.cpp ReconstructSurface.cpp ParticleNormalEvaluation.cpp ParticleFile.cpp @@ -14,6 +15,7 @@ set(Particles_headers ParticleShapeStatistics.h EvaluationUtil.h ShapeEvaluation.h + CorrespondenceEvaluation.h ReconstructSurface.h ParticleNormalEvaluation.h ParticleFile.h @@ -32,6 +34,7 @@ target_include_directories(Particles PUBLIC target_link_libraries(Particles PUBLIC Mesh Optimize + Project tinyxml Eigen3::Eigen ) diff --git a/Libs/Particles/CorrespondenceEvaluation.cpp b/Libs/Particles/CorrespondenceEvaluation.cpp new file mode 100644 index 00000000000..39c59c729c0 --- /dev/null +++ b/Libs/Particles/CorrespondenceEvaluation.cpp @@ -0,0 +1,252 @@ +#include "CorrespondenceEvaluation.h" + +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include + +namespace shapeworks { + +namespace { + +Mesh load_groomed_as_mesh(const std::string& path) { + const std::string ext = StringUtils::getLowerExtension(path); + if (ext == ".nrrd" || ext == ".mha" || ext == ".mhd" || ext == ".nii" || ext == ".gz" || ext == ".tif" || + ext == ".tiff") { + Image img(path); + return img.toMesh(0.0); + } + return Mesh(path); +} + +Eigen::MatrixXd load_particles_matrix(const std::string& filename) { + Eigen::VectorXd v = particles::read_particles(filename); + const int n = static_cast(v.size() / 3); + Eigen::MatrixXd m(n, 3); + for (int i = 0; i < n; ++i) { + m(i, 0) = v(3 * i + 0); + m(i, 1) = v(3 * i + 1); + m(i, 2) = v(3 * i + 2); + } + return m; +} + +CorrespondenceQualityStats summarize(std::vector values) { + CorrespondenceQualityStats s; + if (values.empty()) return s; + std::sort(values.begin(), values.end()); + s.mean = std::accumulate(values.begin(), values.end(), 0.0) / values.size(); + s.max = values.back(); + s.median = values[values.size() / 2]; + const size_t p95_idx = std::min(values.size() - 1, static_cast(0.95 * values.size())); + s.p95 = values[p95_idx]; + return s; +} + +} // namespace + +CorrespondenceQualityReport CorrespondenceEvaluation::evaluate(ProjectHandle project, DistanceMethod method, + const std::string& output_meshes_dir) { + if (!project) { + throw std::runtime_error("null project"); + } + + auto subjects = project->get_non_excluded_subjects(); + if (subjects.empty()) { + throw std::runtime_error("no (non-excluded) subjects in project"); + } + const int num_domains = project->get_number_of_domains_per_subject(); + if (num_domains <= 0) { + 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. + std::vector name_per_subject; + std::vector> particles_per_subject_domain; + std::vector> groomed_per_subject_domain; + + for (auto& subj : subjects) { + auto particle_files = subj->get_local_particle_filenames(); + auto groomed_files = subj->get_groomed_filenames(); + + bool complete = true; + std::vector subj_particles; + std::vector subj_groomed; + subj_particles.reserve(num_domains); + subj_groomed.reserve(num_domains); + + for (int d = 0; d < num_domains; ++d) { + if (d >= static_cast(particle_files.size()) || particle_files[d].empty() || + d >= static_cast(groomed_files.size()) || groomed_files[d].empty()) { + SW_LOG("Skipping subject '{}': missing particles or groomed for domain {}", subj->get_display_name(), d); + complete = false; + break; + } + try { + subj_particles.push_back(load_particles_matrix(particle_files[d])); + } catch (const std::exception& e) { + SW_LOG("Skipping subject '{}' domain {}: particle load failed: {}", subj->get_display_name(), d, e.what()); + complete = false; + break; + } + subj_groomed.push_back(groomed_files[d]); + } + + if (!complete) continue; + name_per_subject.push_back(subj->get_display_name()); + particles_per_subject_domain.push_back(std::move(subj_particles)); + groomed_per_subject_domain.push_back(std::move(subj_groomed)); + } + + const int num_subjects = static_cast(name_per_subject.size()); + if (num_subjects < 2) { + throw std::runtime_error("need at least 2 subjects with complete data across all domains (have " + + std::to_string(num_subjects) + ")"); + } + + // Verify particle counts match across subjects per domain. + std::vector num_particles_per_domain(num_domains); + for (int d = 0; d < num_domains; ++d) { + num_particles_per_domain[d] = static_cast(particles_per_subject_domain[0][d].rows()); + for (int s = 1; s < num_subjects; ++s) { + if (particles_per_subject_domain[s][d].rows() != num_particles_per_domain[d]) { + throw std::runtime_error("domain " + std::to_string(d) + ": subject '" + name_per_subject[s] + "' has " + + std::to_string(particles_per_subject_domain[s][d].rows()) + " particles, expected " + + std::to_string(num_particles_per_domain[d])); + } + } + } + + // L1-medoid template selection over concatenated per-domain local particles + // (matches ParticleShapeStatistics::compute_median_shape). + int template_idx = 0; + double best_l1_sum = std::numeric_limits::infinity(); + for (int i = 0; i < num_subjects; ++i) { + double sum_l1 = 0.0; + for (int j = 0; j < num_subjects; ++j) { + if (i == j) continue; + double pair_l1 = 0.0; + for (int d = 0; d < num_domains; ++d) { + pair_l1 += (particles_per_subject_domain[i][d] - particles_per_subject_domain[j][d]).cwiseAbs().sum(); + } + sum_l1 += pair_l1; + } + if (sum_l1 < best_l1_sum) { + best_l1_sum = sum_l1; + template_idx = i; + } + } + SW_LOG("Template subject (L1-medoid): '{}' (sum of L1 distances to others = {:.4f})", name_per_subject[template_idx], + best_l1_sum); + + // Resolve output meshes dir. + boost::filesystem::path meshes_dir; + if (!output_meshes_dir.empty()) { + meshes_dir = boost::filesystem::path(output_meshes_dir); + boost::filesystem::create_directories(meshes_dir); + SW_LOG("Writing reconstructed meshes to: {}", meshes_dir.string()); + } + + CorrespondenceQualityReport report; + report.template_subject = name_per_subject[template_idx]; + std::vector all_means; // pooled raw mean distances (template excluded) + std::vector 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); + + Mesh template_mesh = load_groomed_as_mesh(groomed_per_subject_domain[template_idx][domain]); + MeshWarper warper; + warper.set_warp_method(WarpMethod::Biharmonic); + warper.set_reference_mesh(template_mesh.getVTKMesh(), particles_per_subject_domain[template_idx][domain]); + if (!warper.generate_warp()) { + throw std::runtime_error("domain " + std::to_string(domain) + ": failed to generate warp from template '" + + name_per_subject[template_idx] + "'"); + } + + for (int i = 0; i < num_subjects; ++i) { + vtkSmartPointer reconstructed = warper.build_mesh(particles_per_subject_domain[i][domain]); + if (!reconstructed) { + SW_LOG("Domain {}: subject '{}' reconstruction returned null, skipping", domain, name_per_subject[i]); + 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; + + 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"); + 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) { + fname = name_per_subject[i] + "_domain" + std::to_string(domain) + "_reconstructed_TEMPLATE.vtk"; + } + boost::filesystem::path out_mesh = meshes_dir / fname; + recon_mesh.write(out_mesh.string()); + } + } + } + + if (report.rows.empty()) { + throw std::runtime_error("no subjects evaluated"); + } + + report.num_template_rows = static_cast(report.rows.size() - all_means.size()); + report.num_evaluated = static_cast(all_means.size()); + report.agg_raw = summarize(all_means); + report.agg_norm = summarize(all_norm_means); + return report; +} + +} // namespace shapeworks diff --git a/Libs/Particles/CorrespondenceEvaluation.h b/Libs/Particles/CorrespondenceEvaluation.h new file mode 100644 index 00000000000..3de1a61a9d5 --- /dev/null +++ b/Libs/Particles/CorrespondenceEvaluation.h @@ -0,0 +1,78 @@ +#pragma once + +#include +#include +#include +#include + +namespace shapeworks { + +class Project; +using ProjectHandle = std::shared_ptr; + +//! Per-subject-per-domain correspondence quality result row. +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 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_max = 0.0; //!< max_dist / bbox_diag + bool is_template = false; //!< true for the L1-medoid template row (excluded from aggregates) +}; + +//! Aggregate summary statistics. +struct CorrespondenceQualityStats { + double mean = 0.0; + double median = 0.0; + double p95 = 0.0; + double max = 0.0; +}; + +//! Full evaluation report. +struct CorrespondenceQualityReport { + std::vector rows; + std::string template_subject; + int num_evaluated = 0; //!< rows.size() - template rows + int num_template_rows = 0; + CorrespondenceQualityStats agg_raw; //!< aggregates over raw mean_dist (template excluded) + CorrespondenceQualityStats agg_norm; //!< aggregates over bbox-normalized values (template excluded) +}; + +/** + * \class CorrespondenceEvaluation + * \ingroup Group-Particles + * + * Per-subject correspondence-quality metric: reconstruct each subject's shape + * from its local particles via biharmonic mesh warp from the cohort L1-medoid + * template (matches Studio's median-subject selection), then measure distance + * from the reconstruction to that subject's groomed mesh. Distances are also + * normalized by each subject's bounding-box diagonal so the metric is + * scale-invariant. + * + * 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. + */ +class CorrespondenceEvaluation { + public: + enum class DistanceMethod { PointToCell, PointToPoint }; + + //! Evaluate. Project is expected to be already loaded, and the current + //! working directory must be one from which the project's relative paths + //! (groomed, local particles) resolve. + //! + //! If \p output_meshes_dir is non-empty, per-subject reconstructed meshes + //! are written there as .vtk with an embedded per-vertex "distance" field. + //! The path is used verbatim (interpreted relative to the current CWD if + //! not absolute). + //! + //! Throws std::runtime_error on setup failures (no subjects, warp failure, + //! mismatched particle counts across subjects). + static CorrespondenceQualityReport evaluate(ProjectHandle project, + DistanceMethod method = DistanceMethod::PointToCell, + const std::string& output_meshes_dir = ""); +}; + +} // namespace shapeworks diff --git a/Libs/Python/ShapeworksPython.cpp b/Libs/Python/ShapeworksPython.cpp index 4d6eda78096..0a5ce1d30ce 100644 --- a/Libs/Python/ShapeworksPython.cpp +++ b/Libs/Python/ShapeworksPython.cpp @@ -44,6 +44,7 @@ using namespace pybind11::literals; #include "PythonAnalyze.h" #include "PythonGroom.h" #include "ReconstructSurface.h" +#include "CorrespondenceEvaluation.h" #include "ShapeEvaluation.h" #include "Shapeworks.h" #include "ShapeworksUtils.h" @@ -1313,6 +1314,43 @@ PYBIND11_MODULE(shapeworks_py, m) { "Computes the specificity measure for a particle system, all modes", "particleSystem"_a, "progress_callback"_a = nullptr, "check_abort"_a = nullptr, "surface_distance_mode"_a = false); + // CorrespondenceEvaluation + py::class_(m, "CorrespondenceQualityRow") + .def_readonly("subject", &CorrespondenceQualityRow::subject) + .def_readonly("domain", &CorrespondenceQualityRow::domain) + .def_readonly("mean_dist", &CorrespondenceQualityRow::mean_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_max", &CorrespondenceQualityRow::norm_max) + .def_readonly("is_template", &CorrespondenceQualityRow::is_template); + + py::class_(m, "CorrespondenceQualityStats") + .def_readonly("mean", &CorrespondenceQualityStats::mean) + .def_readonly("median", &CorrespondenceQualityStats::median) + .def_readonly("p95", &CorrespondenceQualityStats::p95) + .def_readonly("max", &CorrespondenceQualityStats::max); + + py::class_(m, "CorrespondenceQualityReport") + .def_readonly("rows", &CorrespondenceQualityReport::rows) + .def_readonly("template_subject", &CorrespondenceQualityReport::template_subject) + .def_readonly("num_evaluated", &CorrespondenceQualityReport::num_evaluated) + .def_readonly("num_template_rows", &CorrespondenceQualityReport::num_template_rows) + .def_readonly("agg_raw", &CorrespondenceQualityReport::agg_raw) + .def_readonly("agg_norm", &CorrespondenceQualityReport::agg_norm); + + py::class_ corresp_eval(m, "CorrespondenceEvaluation"); + + py::enum_(corresp_eval, "DistanceMethod") + .value("PointToCell", CorrespondenceEvaluation::DistanceMethod::PointToCell) + .value("PointToPoint", CorrespondenceEvaluation::DistanceMethod::PointToPoint) + .export_values(); + + corresp_eval.def_static( + "evaluate", &CorrespondenceEvaluation::evaluate, + "Evaluate per-subject correspondence quality (biharmonic mesh warp from L1-medoid template).", + "project"_a, "method"_a = CorrespondenceEvaluation::DistanceMethod::PointToCell, "output_meshes_dir"_a = ""); + py::class_(m, "ParticleShapeStatistics") .def(py::init<>()) From 6e9eaa2019b6bf56fe7378d889710cd02fbaf8e8 Mon Sep 17 00:00:00 2001 From: Alan Morris Date: Thu, 23 Jul 2026 17:00:32 -0600 Subject: [PATCH 3/4] Link Project as PRIVATE to Particles Project is only referenced inside CorrespondenceEvaluation.cpp, not in any header. Linking it PUBLIC to Particles was propagating it as a transitive dep of every consumer (Optimize, shapeworks_py, tests) and perturbing VTK static-initialization / factory-registration order in the shared-library graph, which broke .vtp reads and mesh-domain tests on Linux and Mac Arm64 in CI while working locally. --- Libs/Particles/CMakeLists.txt | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/Libs/Particles/CMakeLists.txt b/Libs/Particles/CMakeLists.txt index 63095d5006c..93c33344e9b 100644 --- a/Libs/Particles/CMakeLists.txt +++ b/Libs/Particles/CMakeLists.txt @@ -34,11 +34,17 @@ target_include_directories(Particles PUBLIC target_link_libraries(Particles PUBLIC Mesh Optimize - Project tinyxml Eigen3::Eigen ) +# Project is only used inside CorrespondenceEvaluation.cpp — keep PRIVATE so it +# does not propagate to consumers of Particles (which was found to perturb VTK +# static-init/link ordering in shapeworks_py and break .vtp reads). +target_link_libraries(Particles PRIVATE + Project + ) + # set set_target_properties(Particles PROPERTIES PUBLIC_HEADER "${Particles_headers}" From 1467e8417587b8e65871003d0c24685492d828bf Mon Sep 17 00:00:00 2001 From: Alan Morris Date: Fri, 24 Jul 2026 11:29:42 -0600 Subject: [PATCH 4/4] Fix .vtp-read segfault from xlnt/expat symbol collision Linking Project into Particles hoisted libxlnt (which bundles expat and exports the public XML_* symbols) ahead of VTK's/ITK's static expat archives, so xlnt's expat interposed VTK's and crashed on the first .vtp read (OptimizeTests, shapeworks_py). Drop the Particles->Project link edge: CorrespondenceEvaluation.o is only pulled from libParticles.a by the CLI and Python module, which already link Project. Order Particles before Project in shapeworks_py so its symbols still resolve on Linux. --- Libs/Particles/CMakeLists.txt | 18 ++++++++++++------ Libs/Python/CMakeLists.txt | 2 +- 2 files changed, 13 insertions(+), 7 deletions(-) diff --git a/Libs/Particles/CMakeLists.txt b/Libs/Particles/CMakeLists.txt index 93c33344e9b..d6374f4c060 100644 --- a/Libs/Particles/CMakeLists.txt +++ b/Libs/Particles/CMakeLists.txt @@ -38,12 +38,18 @@ target_link_libraries(Particles PUBLIC Eigen3::Eigen ) -# Project is only used inside CorrespondenceEvaluation.cpp — keep PRIVATE so it -# does not propagate to consumers of Particles (which was found to perturb VTK -# static-init/link ordering in shapeworks_py and break .vtp reads). -target_link_libraries(Particles PRIVATE - Project - ) +# NOTE: CorrespondenceEvaluation.cpp uses Project, but we deliberately do NOT +# link Particles against Project here. Project pulls in xlnt, which bundles its +# own copy of expat and exports the public XML_* symbols. Adding a Project -> +# xlnt edge to this low-level library hoists xlnt ahead of VTK's/ITK's static +# expat archives in every consumer's link line, corrupting expat's encoding +# tables and segfaulting on the first .vtp read (broke OptimizeTests and +# shapeworks_py's .vtp reads). Because CorrespondenceEvaluation is only +# referenced by the shapeworks CLI and the Python module -- both of which +# already link Project -- its object is only pulled from libParticles.a by +# those targets, and its Project symbols resolve there. Targets that link +# Particles but never reference CorrespondenceEvaluation (OptimizeTests, etc.) +# never pull the object and keep their original, safe link order. # set set_target_properties(Particles PROPERTIES PUBLIC_HEADER diff --git a/Libs/Python/CMakeLists.txt b/Libs/Python/CMakeLists.txt index 44120688e06..2d33cff798e 100644 --- a/Libs/Python/CMakeLists.txt +++ b/Libs/Python/CMakeLists.txt @@ -16,8 +16,8 @@ target_link_libraries(shapeworks_py PUBLIC Mesh Optimize Groom - Project Particles + Project ${ITK_LIBRARIES} Eigen3::Eigen ${DTAG_EXTRA_FLAGS}