Add conversion of 3D label maps to labeled meshes with SurfaceNets.

[jmargeta]

Overview

This PR adds support for conversion of 3D label maps to labeled surfaces for colored visualization using vtkSurfaceNets3D algorithm from VTK. This is particularly useful for visualization of anatomical structures in medical images from their labels or the outputs of multi-label segmentation maps.

Resolves #5170

Details

• added contour_labeled to ImageData - contours can be extracted with label_map.contour_labeled()
• added example into the gallery with the frog tissue labels and a simple test for the number of extracted labels in the output mesh
  
• Note, vtk 9.3+ must be used which is currently in the pre-release stage on PyPI. Otherwise, the support with throw an exception. The example will, however, not be generated, and the test will be blank.

As mentioned in the issue, another appealing example could be made using the great Neuroparc project. Is it acceptable to make an example with an external dataset?

import pooch
import pyvista as pv
file_name = pooch.retrieve(
    "https://github.com/neurodata/neuroparc/raw/master/atlases/label/Human/Schaefer1000_space-MNI152NLin6_res-1x1x1.nii.gz",
    known_hash="c3efe797aab3b3d9e705645bf29fac4e932c88dbe54ccaeb03982f11e66b3249",
)
label_map = pv.read(file_name)
label_map.contour_labeled(smoothing=True).plot(cmap="glasbey_bw")

Does the name contour_labeled fit well with PyVista, or would you suggest a better one?

Also, I added types to the function, and had to silence these Mypy errors (as in #3837):

pyvista/core/filters/image_data.py:853: error: Argument 1 to "set_default_active_scalars" has incompatible type "ImageDataFilters"; expected "DataSet"  [arg-type]
pyvista/core/filters/image_data.py:854: error: "ImageDataFilters" has no attribute "active_scalars_info"  [attr-defined]
pyvista/core/filters/image_data.py:858: error: "ImageDataFilters" has no attribute "get_array_association"  [attr-defined]

More info about the algorithm: Sarah F. Frisken, SurfaceNets for Multi-Label Segmentations with Preservation of Sharp Boundaries, Journal of Computer Graphics Techniques (JCGT), vol. 11, no. 1, 34-54, 2022. Available online http://jcgt.org/published/0011/01/03/

[codecov]

Codecov Report

Merging #5176 (b1bb05e) into main (4e804c4) will increase coverage by 0.00%.
The diff coverage is 100.00%.

Additional details and impacted files

@@           Coverage Diff           @@
##             main    #5176   +/-   ##
=======================================
  Coverage   96.11%   96.11%           
=======================================
  Files         131      131           
  Lines       21877    21914   +37     
=======================================
+ Hits        21026    21063   +37     
  Misses        851      851           

[user27182]

Awesome feature! Thanks for your contribution. This should be linked with #4802 in terms of PyVista VTK 9.3 support. PackLabels is used internally by SurfaceNets.

[jmargeta]

Thanks, @user27182, for the pointer and for linking the two as it should be! Great to see the support for PackLabels.
Once your PR gets merged, I will be happy to adopt your vtk 9.3 detection logic.
Cannot wait to see the latest capabilities of vtk in PyVista.

[MatthewFlamm]

vtk 9.3 is now in the CI, so it can be tested directly

[jmargeta]

@MatthewFlamm, thanks a lot for the update on the CI and your review!
I believe to have addressed all your comments and consider the wrapper complete functionally-wise.

Any feedback on what else I can do here to make the API better, please? In particular, does the contour_labeled name align well with the spirit of PyVista? I chose it primarily for easier discoverability under the contour«TAB»-completion.

def contour_labeled(
    self,
    n_labels: Optional[int] = None,
    smoothing: bool = False,
    smoothing_num_iterations: int = 50,
    smoothing_relaxation_factor: float = 0.5,
    smoothing_constraint_distance: float = 1,
    output_mesh_type: str = 'quads',
    output_style: str = 'default',
    scalars: Optional[str] = None,
    progress_bar: bool = False,
) -> 'pyvista.PolyData'

More explicit alternatives could be something like segmentaton, label_contours, label_map, surface_nets, label, label_surfaces, extract_label_map, segmentation_contours...

@tkoyama010, any chance for this to get into the November Release (#5138)?

[tkoyama010]

@jmargeta I added this PR to the waiting list in the release. Let's add this to the next release!

[jmargeta]

Awesome, thanks @tkoyama010 !

[banesullivan]

Really awesome work! Happy to approve this, but I'm wondering if you can bump up the code coverage in the tests since some of the conditional options are not currently covered.

Otherwise small comments

[jmargeta]

How about the name of the function itself? labels.contour_labeled() works well.
However, after using this internally for a while, it feels that labels.surfaces_labeled() or just labels.surfaces() would be slightly more ergonomic alternatives.

'contour' suggests computation of an iso-contour. Here, 'labels' are a label map as an output from an image segmentation algorithm (3D image with discrete labels per voxel) we want to extract a surface "per label" rather than an iso-contour.
What do you think? Or is the original name fine?

[banesullivan]

How about the name of the function itself? labels.contour_labeled() works well.

+1 for contour_labeled as that is consistent with the few other contouring methods and makes it a bit more discovered with tab-complete. see also: contour_banded()

In my mind "contouring" is simply a geometry reducing operation: volume -> surface, surface -> line, line -> point. Since this is taking a volume (ImageData) and producing surfaces (PolyData), this naming feels appropriate to me.

Otherwise, since this is a direct wrapping of vtkSurfaceNets3D I don't think anyone could argue against surface_nets()

[jmargeta]

@banesullivan, thank you very much for your thoughts. contour_labeled it is, and I am fairly happy with the PR now.
Please, do not hesitate to ping me if you see anything else to add/change.

[banesullivan]

LGTM! Thanks for your contributions, @jmargeta!

[eigenvivek]

Thank you for this PR, it is very very helpful!

When I use contour_labeled, it prints some internal messages from VTK. Is there a way to suppress this by default? I don't want this showing up on a docs page I'm building. Also, messages like this don't show up when you use marching cubes, for example.

Thank you very much for your work, have been looking forward to the possibility of using SurfaceNets in PyVista for a few months now!

[banesullivan]

@eigenvivek, can you share a pv.Report()? I mainly want to know what your VTK version is

[jmargeta]

Hi @eigenvivek,
Great to hear you are using this functionality!

You can try to change the vtk logging level before the contour_labeled execution. For example, this mutes the messages for me:

import vtk
vtk.vtkLogger.SetStderrVerbosity(vtk.vtkLogger.ConvertToVerbosity(-1))

Disclaimer, I am not sure, this is the official way to do it :)

[jmargeta]

All this is happening within the compiled part of VTK.
vtkLog(INFO, "Executing Surface Nets 3D");

[eigenvivek]

@jmargeta thank you for the suggestion, that worked perfectly :)

@banesullivan using vtk==9.3.0

[banesullivan]

Glad you got it fixed! FYI, this discussion also has some tips for suppressing these logs: #3495

@jmargeta and @eigenvivek, if either of you think we should perform this suppression in this filter directly, don't hesitate to add that

[jmargeta]

Hi @banesullivan, thanks for the pointer, vtk.vtkLogger.SetStderrVerbosity(vtk.vtkLogger.VERBOSITY_OFF) from the discussion you mentioned looks good. PR incoming.