ENH: Rewrite PCA model registration around an analytic-gradient objec… - #117
ENH: Rewrite PCA model registration around an analytic-gradient objec…#117aylward wants to merge 1 commit into
Conversation
…tive
RegisterModelsPCA previously maximized mean intensity sampled from the fixed
distance map through ITK's LinearInterpolateImageFunction, with the optimizer
estimating gradients by finite differences. The objective is now stated and
minimized directly:
mean distance(model -> target)
+ w * mean distance(target -> model) # symmetric term
+ lambda * sum(b_i^2) # Mahalanobis shape prior
Because the PCA deformation is linear in the coefficients b, the gradient is
analytic and handed to the optimizer instead of being estimated, which removes
one objective evaluation per coefficient per step. The post-PCA transform is
folded into the mode directions so the gradient stays exact; when that
transform is not affine its Jacobian is not constant, so the analytic gradient
is disabled and finite differences are used with a logged warning.
- Add symmetric_weight (default 0.5) so partial target coverage is penalized,
and pca_prior_weight (default 0.0, disabled) for the Mahalanobis prior
- Sample the distance map and its gradient with scipy.ndimage.map_coordinates,
and build the target-to-model term with a scipy.spatial.cKDTree
- Replace the cached ITK interpolator and _create_itk_points with
_prepare_sampling, which builds the arrays the objective is made of once
- Add ContourTools.sample_mesh_faces for face-density-aware point sampling and
a negative_inside option on the signed distance map
- Log transform fidelity after computing the PCA transforms
Tests cover the pieces that were previously unverified: the analytic gradient
against finite differences, recovery of known coefficients, the symmetric term
penalizing partial coverage, the prior shrinking coefficients, eigenvector
scaling by standard deviation, deformation happening in the template frame,
and a transform round trip.
|
Warning Review limit reached
Next review available in: 59 minutes You've used all free OSS reviews for now. Wait for the free limit to reset to keep reviewing this public repository. How can I continue?After more reviews become available, a review can be triggered using the To avoid repeated limits, reduce automatic review volume by pausing incremental auto-reviews earlier, using label-based review opt-in, excluding WIP or generated PR titles, or requesting reviews manually when the PR is ready. If your team needs uninterrupted high-volume reviews, an organization admin can enable usage-based reviews. How do review limits work?CodeRabbit enforces per-developer PR review limits for each organization. Most developers receive the normal plan review availability. For paid Pro and Pro+ PR reviews, CodeRabbit uses adaptive limits for sustained high-volume activity. When a developer's recent PR review activity reaches the 95th percentile or higher among CodeRabbit users, additional reviews become available more gradually as earlier reviews age out of the rolling window. Please refer docs for additional details. Review details⚙️ Run configurationConfiguration used: Organization UI Review profile: CHILL Plan: Pro Plus Run ID: 📒 Files selected for processing (3)
Thanks for using CodeRabbit! It's free for OSS, and your support helps us grow. If you like it, consider giving us a shout-out. Comment |
|
Combined with PR #118 |
There was a problem hiding this comment.
Pull request overview
Refactors RegisterModelsPCA to directly minimize a symmetric distance-based objective with an (optionally) analytic gradient, improving optimizer efficiency and enabling additional regularization terms (symmetric coverage + Mahalanobis prior). This aligns the PCA registration implementation with a clearer “distance-to-target in mm” formulation and adds targeted tests for previously unverified behaviors.
Changes:
- Replaces ITK interpolator-based sampling with cached NumPy/SciPy sampling (
map_coordinates) plus analytic objective gradient when post-PCA transform is affine. - Adds a symmetric target-to-model term (via
cKDTree) and an optional PCA prior term (pca_prior_weight). - Improves distance-map construction robustness by optionally sampling triangle faces (
ContourTools.sample_mesh_faces) and adds extensive synthetic test coverage.
Reviewed changes
Copilot reviewed 3 out of 3 changed files in this pull request and generated 6 comments.
| File | Description |
|---|---|
src/physiotwin4d/register_models_pca.py |
New objective definition, sampling caches, analytic gradient path, symmetric/prior terms, and transform fidelity logging. |
src/physiotwin4d/contour_tools.py |
Adds face-aware mesh sampling and updates distance-map rasterization to use those samples (plus a new sample_faces flag). |
tests/test_register_models_pca.py |
Adds synthetic tests for mode scaling, frame correctness, analytic gradient, symmetric penalty, prior shrinkage, and transform round trips. |
Suppressed comments (6)
src/physiotwin4d/register_models_pca.py:38
- Non-ASCII characters in .py files are disallowed in this repo; this docstring line uses Σ and a superscript ². Please replace with an ASCII-only expression (e.g.,
sum(b_i**2)).
+ lambda * Σ b_i² # Mahalanobis prior
src/physiotwin4d/register_models_pca.py:53
- Non-ASCII glyphs are disallowed in .py files in this repo; this docstring line uses the multiplication sign (×). Please switch to ASCII (e.g.,
xor(modes, n_points*3)).
pca_eigenvectors (np.ndarray): PCA eigenvectors/components (modes × n_points*3)
src/physiotwin4d/register_models_pca.py:175
- Non-ASCII glyphs are disallowed in .py files in this repo; this error message uses the multiplication sign (×). Please switch to ASCII for robustness on Windows encoding.
f"Component dimension mismatch: expected {expected_size} "
f"(3 × {pca_template_model.n_points} points), got "
f"{self.pca_eigenvectors.shape[1]}"
src/physiotwin4d/register_models_pca.py:821
- Non-ASCII glyphs are disallowed in .py files in this repo; this log message uses the ± symbol. Please switch to ASCII (e.g., '+/-') for robustness.
self.log_info(
f"PCA coefficient bounds: ±{pca_coefficient_bounds} std deviations"
)
src/physiotwin4d/register_models_pca.py:792
- Non-ASCII glyphs are disallowed in .py files in this repo; this docstring line uses the ± symbol. Please switch to ASCII (e.g., '+/-3').
pca_coefficient_bounds: Bound on PCA coefficients in units of std deviations.
Default: 3.0 (±3 std deviations per mode)
src/physiotwin4d/register_models_pca.py:1048
- Non-ASCII glyphs are disallowed in .py files in this repo; this docstring line uses the ± symbol. Please switch to ASCII (e.g., '+/- std devs').
pca_number_of_modes: Number of PCA modes to use. Default: 0 (use all available modes)
pca_coefficient_bounds: PCA coefficient bounds (±std devs). Default: 3.5
method: Optimization method for scipy.optimize.minimize.
💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.
| - Optimizes PCA coefficients | ||
| - Model equation: P = mean + Σ(b_i * std_i * pca_eigenvector_i) | ||
| - Maximizes mean distance at deformed model points P | ||
| - Model equation: P = template + Σ(b_i * std_i * pca_eigenvector_i) |
| # Target points for the symmetric term. | ||
| self._target_points: Optional[np.ndarray] = None | ||
| if self.symmetric_weight > 0.0: | ||
| if self.fixed_model is None: | ||
| self.log_warning( | ||
| "symmetric_weight is %.3g but no fixed_model is available; " | ||
| "the target-to-model term is disabled.", | ||
| self.symmetric_weight, | ||
| ) | ||
| else: | ||
| self._target_points = np.asarray( | ||
| self.fixed_model.points, dtype=np.float64 | ||
| )[self._sample_slice] | ||
|
|
| return None | ||
| if self._post_pca_affine_key is not self.post_pca_transform: | ||
| self._post_pca_affine_key = self.post_pca_transform | ||
| self._post_pca_affine = self._affine_of_transform(self.post_pca_transform) | ||
| return self._post_pca_affine |
| options: dict = {"maxiter": max_iterations, "disp": disp, "gtol": 1e-6} | ||
| if not self._analytic_gradient: | ||
| options["eps"] = 1e-2 |
| points = np.asarray(mesh.points, dtype=np.float64) | ||
| surface = mesh.extract_surface() if not isinstance(mesh, pv.PolyData) else mesh | ||
| surface = surface.triangulate() |
| u, v = np.meshgrid(steps, steps, indexing="ij") | ||
| mask = (u + v) <= 1.0 | ||
| weights = np.column_stack([1.0 - u[mask] - v[mask], u[mask], v[mask]]) | ||
| samples.append(np.einsum("fca,kc->fka", selected, weights).reshape(-1, 3)) |
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #117 +/- ##
==========================================
+ Coverage 36.60% 39.50% +2.89%
==========================================
Files 72 72
Lines 8510 8627 +117
==========================================
+ Hits 3115 3408 +293
+ Misses 5395 5219 -176
Flags with carried forward coverage won't be shown. Click here to find out more. ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
…tive
RegisterModelsPCA previously maximized mean intensity sampled from the fixed distance map through ITK's LinearInterpolateImageFunction, with the optimizer estimating gradients by finite differences. The objective is now stated and minimized directly:
Because the PCA deformation is linear in the coefficients b, the gradient is analytic and handed to the optimizer instead of being estimated, which removes one objective evaluation per coefficient per step. The post-PCA transform is folded into the mode directions so the gradient stays exact; when that transform is not affine its Jacobian is not constant, so the analytic gradient is disabled and finite differences are used with a logged warning.
Tests cover the pieces that were previously unverified: the analytic gradient against finite differences, recovery of known coefficients, the symmetric term penalizing partial coverage, the prior shrinking coefficients, eigenvector scaling by standard deviation, deformation happening in the template frame, and a transform round trip.