Add skin surface computation tools (#149)

[ebrahimebrahim]

Close #149

The goal of this PR is to add tooling to openlifu to create the skin surface interpolator needed for #147 (virtual fitting).

A useful bonus is that we will also get a skin surface model that can be used in transducer tracking.

Roughly the workflow from an MRI volume array and associated affine is this:

graph LR
A[array] --> B{"(1) fg mask"} --> C[mask] --> E{"(2) meshify"}
D[affine]-->E-->F[mesh]-->G{"(3)make<br>interpolator"}
G-->H[interpolator]
I[origin,<br>axes]-->H[interpolator]

Loading

The goal of this PR is to introduce algorithms (1), (2), and (3).

(1) compute_foreground_mask

compute_foreground_mask is a python implementation of the BRAINSTools foreground masking algorithm, with some simplifications and (hopefully) improvements

• Original algorithm documentation: https://slicer.readthedocs.io/en/latest/user_guide/modules/brainsroiauto.html
• Original algorithm code: https://github.com/BRAINSia/BRAINSTools/tree/7c37d9e8c238f66f8a83f997d9c9bb659c494c90/BRAINSROIAuto

(2) create_closed_surface_from_labelmap

create_closed_surface_from_labelmap is based on the closed surface representation functionality in Slicer's segmentation module.

• Original algorithm: https://github.com/Slicer/Slicer/blob/677932127c73a6c78654d4afd9458a655a4eef63/Libs/vtkSegmentationCore/vtkBinaryLabelmapToClosedSurfaceConversionRule.cxx#L246-L476

(3) spherical_interpolator_from_mesh

Here a "spherical interpolator" is a function that maps angles from a spherical coordinate system to r values (radial spherical coordinate values) by interpolating over a set of known values. It's essentially a "spherical plotter."

Summary of the algorithm:

• Transform the input mesh based on the desired origin and orientation of the spherical coordinate system.
• We will gather some points into a set $S$. For each point $P$ on the mesh consider the ray $\vec{OP}$ from the origin through $P$ and look at all the intersections of this ray $\vec{OP}$ with the mesh. If none of those intersections are further out from the origin than $P$ is, then we put $P$ into our set $S$.
• Using the spherical coordinates of the points in $S$, build a scipy.interpolate.LinearNDInterpolator that interpolates spherical $r$ values from the spherical $(\theta,\phi)$ values.
  • Problem: All the gathered $(\theta,\phi)$ values are likely strictly inside the square $[0,\pi]\times[-\pi,\pi]$, and LinearNDInterpolator does not extrapolate, and so angles close to the "seams" of the spherical coordinate system (the boundaries of that square) generate NaNs through the interpolator. The solution used here is to first clone the gathered points with appropriate angular shifts so as to cover those seams, and then give that larger set of points to the interpolator.
• Return the interpolator.

Example visualized

Here is a visualization of the whole pipeline on an example (this one is of a defaced MRI):

We start with a volume (shown as a slice), then we get a binary labelmap (shown as a slice), then we get a mesh (shown as a 3d render), and finally we get a spherical interpolator (shown applied to the vertices of a vtkSphereSource and then 3d rendered).