Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add updated SPHARM-PDM functionality to SlicerSALT build #306

Open
bpaniagua opened this issue Sep 10, 2024 · 2 comments
Open

Add updated SPHARM-PDM functionality to SlicerSALT build #306

bpaniagua opened this issue Sep 10, 2024 · 2 comments
Assignees

Comments

@bpaniagua
Copy link
Contributor

bpaniagua commented Sep 10, 2024

https://github.com/Kitware/SlicerSALT/blob/master/CMakeLists.txt#L399

@allemangD could you please link here the SPHARM-PDM fork that includes your WIP on the python implementation of SPHARM-PDM?


Plan is to include the rework version as an "experimental" alternative for SPHARM-PDM.

@bpaniagua bpaniagua added this to the SlicerSALT 5.0 milestone Sep 10, 2024
@allemangD
Copy link
Collaborator

allemangD commented Sep 10, 2024

https://github.com/slicersalt/py-spharm-pdm/

Two attempts at optimization with branches torch and torch-optim - neither work very well for one reason or another.

Currently the two core challenges:

  1. scipy.optimize uses dense matrices for storing Jacobian and Hessian information, so its constrained optimizers use a prohibitive amount of memory.
    • The branch torch attempts to use pytorch backward autograd to compute these instead. Runtime improves but memory usage still prohibitive.
    • The branch torch-optim attempts to use torch.optim library, however it doesn't support constrained optimizers. I attempted a constraint penalty term, but it doesn't perform well: it optimizes for the constraint or the metric, not both. There are too many degrees of freedom and the constraints are too sparse.
  2. There is a bug somewhere in core.initial_parameterization that causes self-intersecting initial parameterization; since the initial conditions are wrong the optimizer fails to produce any meaningful result.

The problem with the python optimizers is we need sparse and constrained optimization. I'm curious to investigate the package pyOptSparse; it claims to support this. It doesn't provide an autograd interface but I can attempt with the same technique as in the torch branch. I'm not super confident since it's not a terribly popular library but I think it's worth looking at before diving into a custom implementation (that's the current situation with EqualAreaParametricMeshNewtonIterator!!)

The first issue to resolve is the discrepency in initial parameterization. My current Python implementation is based on the paper, not the code. I've isolated the relevant sections of EqualAreaParametricMeshNewtonIterator and am comparing them now. I will follow up here with details, or schedule a breakout with you and Jared if the bug is too subtle.

Note for posterity: the relevant methods in EqualAreaParametricMeshNewtonIterator are:

  • start_values invokes general workflow.
  • generate_matrix populates laplacian matrix for latitude problem.
  • set_lati_rhs populates rhs vector for latitude problem.
  • Solve Ax = b; extract latitudes from x.
  • modify_matrix updates the matrix for the longitude problem.
  • set_longi_rhs populates rhs vector for longitude problem.
  • Solve Ax = b; extract longitudes from x.
  • Compute cartesian coordinates.

@allemangD
Copy link
Collaborator

I believe I located and resolved the bug: slicersalt/py-spharm-pdm@e5a47f5 Now working on more detailed tests and automated validation to be sure, but I wanted to report this progress.

While investigating this I realize that the legacy TSurfaceNet and new vtkSurfaceNets3D output have different vertex ordering; then the spherical parameterization and resulting harmonics will differ. Therefore I do not think it will be possible to maintain correspondence between models generated with this rework and the legacy implementation.

It may be possible to resolve it in certain cases with ellipsoid alignment but I'm not confident. Assuming not, I plan to use model to model distance as the metric to validate.

@bpaniagua bpaniagua removed this from the SlicerSALT 5.0 milestone Oct 7, 2024
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants