From 94d25fa4d6bf25b592691b4dc0b282349d947077 Mon Sep 17 00:00:00 2001
From: nawazishkhan1-nk <nawazishkhan1.nk@gmail.com>
Date: Wed, 25 Mar 2026 21:13:49 -0600
Subject: [PATCH] add coefficients extractions for morphodev

---
 .../EarlyStop/MorphologicalDeviationScore.cpp | 28 ++++++++++++++
 .../EarlyStop/MorphologicalDeviationScore.h   |  9 +++++
 Libs/Python/ShapeworksPython.cpp              | 38 +++++++++++++++++++
 3 files changed, 75 insertions(+)

diff --git a/Libs/Optimize/Function/EarlyStop/MorphologicalDeviationScore.cpp b/Libs/Optimize/Function/EarlyStop/MorphologicalDeviationScore.cpp
index 925a5ac287..b865068c90 100644
--- a/Libs/Optimize/Function/EarlyStop/MorphologicalDeviationScore.cpp
+++ b/Libs/Optimize/Function/EarlyStop/MorphologicalDeviationScore.cpp
@@ -110,4 +110,32 @@ Eigen::VectorXd MorphologicalDeviationScore::GetMorphoDevScore(const Eigen::Matr
   }
 }
 
+//---------------------------------------------------------------------------
+Eigen::MatrixXd MorphologicalDeviationScore::GetPCACoefficients(const Eigen::MatrixXd& X) {
+  try {
+    if (!is_fitted_) {
+      throw std::runtime_error("PPCA model is not fitted on control shapes.");
+    }
+
+    if (all_components_.cols() == 0) {
+      throw std::runtime_error("PPCA basis is empty.");
+    }
+
+    if (X.cols() != mean_.cols()) {
+      throw std::runtime_error("Input feature dimension does not match fitted PPCA model.");
+    }
+
+    if (X.rows() == 0) {
+      return Eigen::MatrixXd(0, all_components_.cols());
+    }
+
+    Eigen::MatrixXd X_bar = X.rowwise() - mean_;  // (n x d)
+    return X_bar * all_components_;               // (n x rank) scores in PCA basis
+
+  } catch (std::exception& e) {
+    SW_ERROR("Exception in computing PCA coefficients for early stopping {}", e.what());
+    return Eigen::MatrixXd();
+  }
+}
+
 }  // namespace shapeworks
diff --git a/Libs/Optimize/Function/EarlyStop/MorphologicalDeviationScore.h b/Libs/Optimize/Function/EarlyStop/MorphologicalDeviationScore.h
index a5316c07fe..537c5a1b75 100644
--- a/Libs/Optimize/Function/EarlyStop/MorphologicalDeviationScore.h
+++ b/Libs/Optimize/Function/EarlyStop/MorphologicalDeviationScore.h
@@ -1,6 +1,7 @@
 #pragma once
 
 #include <Eigen/Dense>
+#include <stdexcept>
 
 namespace shapeworks {
 class MorphologicalDeviationScore {
@@ -11,6 +12,14 @@ class MorphologicalDeviationScore {
   /// Get Mahalanobis-based deviation score for test samples (non-fixed
   /// shapes/domains)
   Eigen::VectorXd GetMorphoDevScore(const Eigen::MatrixXd& X);  // (n,)
+  Eigen::MatrixXd GetPCACoefficients(const Eigen::MatrixXd& X);  // (n, rank)
+  void SetRetainedVarianceRatio(double ratio) {
+    if (ratio <= 0.0 || ratio > 1.0) {
+      throw std::invalid_argument("retained_variance_ratio must be in the interval (0, 1].");
+    }
+    retained_variance_ratio_ = ratio;
+  }
+  double GetRetainedVarianceRatio() const { return retained_variance_ratio_; }
 
  private:
   /// Flag to ensure control shapes are set and PCA model is in place
diff --git a/Libs/Python/ShapeworksPython.cpp b/Libs/Python/ShapeworksPython.cpp
index a63f18aa39..2f392b51be 100644
--- a/Libs/Python/ShapeworksPython.cpp
+++ b/Libs/Python/ShapeworksPython.cpp
@@ -1744,6 +1744,28 @@ PYBIND11_MODULE(shapeworks_py, m) {
                      True if fitting was successful, False otherwise.
              )pbdoc")
 
+      .def("SetRetainedVarianceRatio",
+           &MorphologicalDeviationScore::SetRetainedVarianceRatio, py::arg("ratio"),
+           R"pbdoc(
+                 Set the retained variance ratio used to choose PPCA components.
+
+                 Parameters
+                 ----------
+                 ratio : float
+                     Target cumulative variance ratio in the interval (0, 1].
+             )pbdoc")
+
+      .def("GetRetainedVarianceRatio",
+           &MorphologicalDeviationScore::GetRetainedVarianceRatio,
+           R"pbdoc(
+                 Get the retained variance ratio used to choose PPCA components.
+
+                 Returns
+                 -------
+                 float
+                     The current retained variance ratio.
+             )pbdoc")
+
       .def("GetMorphoDevScore",
            &MorphologicalDeviationScore::GetMorphoDevScore, py::arg("X"),
            R"pbdoc(
@@ -1758,5 +1780,21 @@ PYBIND11_MODULE(shapeworks_py, m) {
                  -------
                  numpy.ndarray
                      Vector of Mahalanobis distances for each sample.
+             )pbdoc")
+
+      .def("GetPCACoefficients",
+           &MorphologicalDeviationScore::GetPCACoefficients, py::arg("X"),
+           R"pbdoc(
+                 Project samples onto the fitted PCA basis.
+
+                 Parameters
+                 ----------
+                 X : numpy.ndarray
+                     Matrix of samples (n_samples x n_features)
+
+                 Returns
+                 -------
+                 numpy.ndarray
+                     Matrix of PCA coefficients with shape (n_samples, rank).
              )pbdoc");
 }  // PYBIND11_MODULE(shapeworks_py)