Repository: pyvista/pyvista

Issue:
Many issues with `contour_labeled`, consider deprecating and rewriting as `contour_labels`
### Describe the bug, what's wrong, and what you expected.

### Overview
There are multiple issues with the current implemention of `contour_labeled`:

- Parameter `n_labels` does not work as expected
- Parameter `output_style` is broken/has no effect (the output is currently always `'boundary'`)
- The output `'BoundaryLabels'` array is a confusing (and undocumented) 2-component array by default, whereas a single component output array is more appropriate IMO as this also matches the shape of the input array
- Parameter `smoothing_constraint_distance` is broken/has no effect

Fixing some or part of this is relatively straightforward (see details), but some fixes, like changing the default size of the output array, are breaking changes or require a lot of effort to support a graceful deprecation.

In addition, adding new parameters to support new features, e.g. as proposed in #5968 , would technically require placing parameters at the end of the kwargs dict due the use of positional kwargs, e.g. see https://github.com/pyvista/pyvista/issues/5506. So this would mean the docs would be a bit out of order.

To address all of these issues, I propose deprecating  `contour_labeled` altogether and implementing the changes / fixes detailed  below as part of a new `contour_labels` filter.

### Details
### Issue with `n_labels`
The docs state:
https://github.com/pyvista/pyvista/blob/a01b7bfeb9b4b348ed5bba17c7985798bd531995/pyvista/core/filters/image_data.py#L835-L836

But the implementation calls `GenerateLabels`:
https://github.com/pyvista/pyvista/blob/a01b7bfeb9b4b348ed5bba17c7985798bd531995/pyvista/core/filters/image_data.py#L917-L918

Which is odd because the [vtk docs for  GenerateLabels](https://vtk.org/doc/nightly/html/classvtkSurfaceNets3D.html#a58cddd531effd33eb3f2da16ca5f594b) state that this is for generating `n` _equally spaced_ labels. 

As such, the  way the current implementation for `n_labels` is written, if I have a label map with two label values, e.g. `[2, 5]`, and I want two labels and set `n_labels=2`, this won't work as expected, because it sets `alg.GenerateLabels(2, 1, 2)`, which will actually sample `2` labels in the range `[1, 2]` and will therefore only return one label, not two (label `5` will be missing).

So, it seems that `alg.GenerateLabels` is more for sub-sampling labels rather than specifying the number of labels to generate surfaces for.

At the least I think the docs would need to be updated to clarify this, but maybe even renaming the parameter to `generate_labels` would be best. This would require deprecating `n_labels`

### Issue with `output_style`

The current implementation is such that setting `output_style` has no effect, and all outputs are actually `boundary`. E.g.:

``` python
import pyvista as pv
import numpy as np
# Create 4x3x3 image with two adjacent labels
dim = (4, 3, 3)
labels = np.zeros(np.prod(dim))
labels[17] = 2
labels[[18,19]] = 5
image = pv.ImageData(dimensions=dim)
image.point_data['labels'] = labels

mesh = image.contour_labeled().plot()
```
Will generate:
<img width="624" alt="image" src="https://github.com/pyvista/pyvista/assets/89109579/f82854da-eedf-43d4-9fda-380a903d3913">

This is _supposed_ to be the `'default'` output, but the default mode is meant to include faces at the boundaries between regions. In this case, there should be a quad cell separating the yellow from the purple regions. But, you can see the far end of the inner purple cell because this boundary cell is missing. So, what you're seeing is actually what `boundary` is supposed to do, where internal boundary cells between regions are excluded.

The correct, expected 'default' output is this:
<img width="624" alt="image" src="https://github.com/pyvista/pyvista/assets/89109579/84fdd914-a3c3-48c0-b3d9-2d24519e3c01">

The reason this isn't currently working is because `SetLabel` is required, e.g. something like this is needed:
``` python
[alg.SetLabel(i, val) for i, val in enumerate(np.unique(self.active_scalars))]
```
To be fair, this is really unclear in the vtk docs, I had to go searching in the vtk tests to find this, see the test code for [`TestSurfaceNets3D`](https://gitlab.kitware.com/vtk/vtk/-/blob/master/Filters/Core/Testing/Python/TestSurfaceNets3D3.py#L77).

In the second image, notice the internal quad separating the two regions. Also, notice the scalar values now are extending above `5`. This is a bit surprising because I would have expected the output scalars to only include `[2, 5]`, nothing else. This brings me to the next issue...

### Issue with output array

By default, the output `BoundaryLables` array is actually a 2-component array. From the [vtkSurfaceNets3D docs](https://vtk.org/doc/nightly/html/classvtkSurfaceNets3D.html#a58cddd531effd33eb3f2da16ca5f594b):
>the filter outputs a two-component, cell data array indicating the labels/regions on either side of the polygons composing the output

This explains the odd scalars, because we have `[2, 0]` or `[5,0]` for most cells that share a boundary with the background, but we have `[2, 5]` at the boundary between these two regions. This gets computed as the scalar value `norm([2, 5])=5.38`. 

In my view, a two-component output like this is a bit confusing (the vtks docs even suggests it's for advanced use). I believe the output should be a single-component by default, with either using two arrays, one for each component, or adding a kwarg to toggle this option on/off.

But, I think changing this default would be a bit complicated in terms of warnings/deprecations.

### Issue with `smoothing_constraint_distance`

`smoothing_constraint_distance` is broken/has no effect since by default the [AutomaticSmoothingConstraints](https://vtk.org/doc/nightly/html/classvtkSurfaceNets3D.html#af7e082c10eaba014ded16f791e0629aa) are on, which mean the manually specified `smoothing_constraint_distance` is ignored. This is a simple bugfix by setting `alg.AutomaticSmoothingConstraintsOff()`.




### Steps to reproduce the bug.

See sample code above

### System Information

```shell
PyVista 0.44.dev0
```


### Screenshots

_No response_

Discussion:

Patch:
diff --git a/examples/01-filter/contouring.py b/examples/01-filter/contouring.py
index 9fa7079a..258a3f18 100644
--- a/examples/01-filter/contouring.py
+++ b/examples/01-filter/contouring.py
@@ -55,7 +55,7 @@ pl.show()
 # %%
 # Banded Contours
 # +++++++++++++++
-# Create banded contours for surface meshes using :func:`contour_banded() <pyvista.PolyDataFilters.contour_banded>`.
+# Create banded contours for surface meshes using :func:`~pyvista.PolyDataFilters.contour_banded`.
 mesh = examples.load_random_hills()
 
 # %%
@@ -83,11 +83,11 @@ pl.show()
 # +++++++++++++++++++++++++
 #
 # Create labeled surfaces from 3D label maps (e.f. multi-label image segmentation)
-# using :func:`contour_labeled() <pyvista.ImageDataFilters.contour_labeled>`.
+# using :func:`~pyvista.ImageDataFilters.contour_labels`.
 # Requires VTK version 9.3
 if pv.vtk_version_info >= (9, 3):
     label_map = pv.examples.load_frog_tissues()
-    mesh = label_map.contour_labeled(smoothing=True)
-    mesh.plot(cmap='glasbey', cpos='yx', show_scalar_bar=False, clim=mesh.get_data_range())
+    mesh = label_map.contour_labels()
+    mesh.plot(cmap='glasbey', cpos='yx', categories=True)
 # %%
 # .. tags:: filter
diff --git a/pyvista/core/filters/image_data.py b/pyvista/core/filters/image_data.py
index f9bd34c2..9edc7689 100644
--- a/pyvista/core/filters/image_data.py
+++ b/pyvista/core/filters/image_data.py
@@ -868,6 +868,24 @@ class ImageDataFilters(DataSetFilters):
         .. note::
            Requires ``vtk>=9.3.0``.
 
+        .. deprecated:: 0.45
+            This filter produces unexpected results and is deprecated.
+            Use :meth:`~pyvista.ImageDataFilters.contour_labels` instead.
+            See https://github.com/pyvista/pyvista/issues/5981 for details.
+
+            To replicate the default behavior from this filter, call `contour_labels`
+            with the following arguments:
+
+            .. code::
+
+                image.contour_labels(
+                    smoothing=False,  # old filter does not apply smoothing
+                    output_mesh_type='quads',  # old filter generates quads
+                    pad_background=False,  # old filter generates open surfaces at input edges
+                    compute_normals=False,  # old filter does not compute normals
+                    multi_component_output=True,  # old filter returns multi-component scalars
+                )
+
         Parameters
         ----------
         n_labels : int, optional
@@ -927,6 +945,13 @@ class ImageDataFilters(DataSetFilters):
             Function used internally by SurfaceNets to generate contiguous label data.
 
         """
+        warnings.warn(
+            'This filter produces unexpected results and is deprecated. Use `contour_labels` instead.'
+            '\nRefer to the documentation for `contour_labeled` for details on how to transition to the new filter.'
+            '\nSee https://github.com/pyvista/pyvista/issues/5981 for details.',
+            PyVistaDeprecationWarning,
+        )
+
         if not hasattr(_vtk, 'vtkSurfaceNets3D'):  # pragma: no cover
             from pyvista.core.errors import VTKVersionError
 
diff --git a/pyvista/core/filters/poly_data.py b/pyvista/core/filters/poly_data.py
index 69c3859e..79b33683 100644
--- a/pyvista/core/filters/poly_data.py
+++ b/pyvista/core/filters/poly_data.py
@@ -4173,7 +4173,7 @@ class PolyDataFilters(DataSetFilters):
         pyvista.voxelize_volume
             Similar function that returns a :class:`~pyvista.RectilinearGrid` with cell data.
 
-        pyvista.ImageDataFilters.contour_labeled
+        pyvista.ImageDataFilters.contour_labels
             Filter that generates surface contours from labeled image data. Can be
             loosely considered as an inverse of this filter.
 
@@ -4272,7 +4272,7 @@ class PolyDataFilters(DataSetFilters):
         an existing mask.
 
         >>> volume = examples.load_frog_tissues()
-        >>> poly = volume.contour_labeled(smoothing=True)
+        >>> poly = volume.contour_labels()
 
         Now create the mask from the polydata using the volume as a reference.
 
diff --git a/pyvista/examples/downloads.py b/pyvista/examples/downloads.py
index 514587eb..5deadf12 100644
--- a/pyvista/examples/downloads.py
+++ b/pyvista/examples/downloads.py
@@ -7545,7 +7545,7 @@ def download_whole_body_ct_male(load=True):  # pragma: no cover
 
     Create a surface mesh of the segmentation labels
 
-    >>> labels_mesh = label_map.contour_labeled(smoothing=True)
+    >>> labels_mesh = label_map.contour_labels(smoothing=True)
 
     Plot the CT image and segmentation labels together.
 
@@ -7689,7 +7689,7 @@ def download_whole_body_ct_female(load=True):  # pragma: no cover
 
     Create a surface mesh of the segmentation labels
 
-    >>> labels_mesh = label_map.contour_labeled(smoothing=True)
+    >>> labels_mesh = label_map.contour_labels(smoothing=True)
 
     Plot the CT image and segmentation labels together.
 


Test patch:
diff --git a/tests/core/test_imagedata_filters.py b/tests/core/test_imagedata_filters.py
index 9fd17971..506e39bd 100644
--- a/tests/core/test_imagedata_filters.py
+++ b/tests/core/test_imagedata_filters.py
@@ -11,7 +11,6 @@ from pyvista import examples
 from pyvista.core._validation._cast_array import _cast_to_tuple
 from pyvista.core.errors import PyVistaDeprecationWarning
 
-VTK93 = pv.vtk_version_info >= (9, 3)
 BOUNDARY_LABELS = 'boundary_labels'
 
 
@@ -48,124 +47,121 @@ def three_dimensionality_image():
     return variable_dimensionality_image(dimensions=(2, 2, 2))
 
 
-@pytest.mark.skipif(not VTK93, reason='At least VTK 9.3 is required')
-def test_contour_labeled():
-    # Load a 3D label map (segmentation of a frog's tissue)
-    label_map = examples.load_frog_tissues()
+@pytest.fixture
+def frog_tissues():
+    return examples.load_frog_tissues()
+
+
+def test_contour_labeled_deprecated():
+    match = 'This filter produces unexpected results and is deprecated.'
+    with pytest.raises(PyVistaDeprecationWarning, match=match):
+        pv.ImageData().contour_labeled()
 
+
+@pytest.mark.needs_vtk_version(9, 3, 0)
+def test_contour_labeled(frog_tissues):
     # Extract surface for each label
-    mesh = label_map.contour_labeled()
+    with pytest.warns(PyVistaDeprecationWarning):
+        mesh = frog_tissues.contour_labeled()
 
-    assert label_map.point_data.active_scalars.max() == 29
+    assert frog_tissues.point_data.active_scalars.max() == 29
     assert 'BoundaryLabels' in mesh.cell_data
     assert np.max(mesh['BoundaryLabels'][:, 0]) == 29
 
 
-@pytest.mark.skipif(not VTK93, reason='At least VTK 9.3 is required')
-def test_contour_labeled_with_smoothing():
-    # Load a 3D label map (segmentation of a frog's tissue)
-    label_map = examples.load_frog_tissues()
-
+@pytest.mark.needs_vtk_version(9, 3, 0)
+def test_contour_labeled_with_smoothing(frog_tissues):
     # Extract smooth surface for each label
-    mesh = label_map.contour_labeled(smoothing=True)
+    with pytest.warns(PyVistaDeprecationWarning):
+        mesh = frog_tissues.contour_labeled(smoothing=True)
     # this somehow mutates the object... also the n_labels is likely not correct
 
     assert 'BoundaryLabels' in mesh.cell_data
     assert np.max(mesh['BoundaryLabels'][:, 0]) == 29
 
 
-@pytest.mark.skipif(not VTK93, reason='At least VTK 9.3 is required')
-def test_contour_labeled_with_reduced_labels_count():
-    # Load a 3D label map (segmentation of a frog's tissue)
-    label_map = examples.load_frog_tissues()
-
+@pytest.mark.needs_vtk_version(9, 3, 0)
+def test_contour_labeled_with_reduced_labels_count(frog_tissues):
     # Extract surface for each label
-    mesh = label_map.contour_labeled(n_labels=2)
+    with pytest.warns(PyVistaDeprecationWarning):
+        mesh = frog_tissues.contour_labeled(n_labels=2)
     # this somehow mutates the object... also the n_labels is likely not correct
 
     assert 'BoundaryLabels' in mesh.cell_data
     assert np.max(mesh['BoundaryLabels'][:, 0]) == 2
 
 
-@pytest.mark.skipif(not VTK93, reason='At least VTK 9.3 is required')
-def test_contour_labeled_with_triangle_output_mesh():
-    # Load a 3D label map (segmentation of a frog's tissue)
-    label_map = examples.load_frog_tissues()
-
+@pytest.mark.needs_vtk_version(9, 3, 0)
+def test_contour_labeled_with_triangle_output_mesh(frog_tissues):
     # Extract surface for each label
-    mesh = label_map.contour_labeled(scalars='MetaImage', output_mesh_type='triangles')
+    with pytest.warns(PyVistaDeprecationWarning):
+        mesh = frog_tissues.contour_labeled(scalars='MetaImage', output_mesh_type='triangles')
 
     assert 'BoundaryLabels' in mesh.cell_data
     assert np.max(mesh['BoundaryLabels'][:, 0]) == 29
 
 
-@pytest.mark.skipif(not VTK93, reason='At least VTK 9.3 is required')
-def test_contour_labeled_with_boundary_output_style():
-    # Load a 3D label map (segmentation of a frog's tissue)
-    label_map = examples.load_frog_tissues()
-
+@pytest.mark.needs_vtk_version(9, 3, 0)
+def test_contour_labeled_with_boundary_output_style(frog_tissues):
     # Extract surface for each label
-    mesh = label_map.contour_labeled(output_style='boundary')
+    with pytest.warns(PyVistaDeprecationWarning):
+        mesh = frog_tissues.contour_labeled(output_style='boundary')
 
     assert 'BoundaryLabels' in mesh.cell_data
     assert np.max(mesh['BoundaryLabels'][:, 0]) == 29
 
 
-@pytest.mark.skipif(not VTK93, reason='At least VTK 9.3 is required')
-def test_contour_labeled_with_invalid_output_mesh_type():
-    # Load a 3D label map (segmentation of a frog's tissue)
-    label_map = examples.load_frog_tissues()
-
+@pytest.mark.needs_vtk_version(9, 3, 0)
+def test_contour_labeled_with_invalid_output_mesh_type(frog_tissues):
     # Extract surface for each label
-    with pytest.raises(ValueError):  # noqa: PT011
-        label_map.contour_labeled(output_mesh_type='invalid')
-
+    with pytest.warns(PyVistaDeprecationWarning):
+        with pytest.raises(ValueError):  # noqa: PT011
+            frog_tissues.contour_labeled(output_mesh_type='invalid')
 
-@pytest.mark.skipif(not VTK93, reason='At least VTK 9.3 is required')
-def test_contour_labeled_with_invalid_output_style():
-    # Load a 3D label map (segmentation of a frog's tissue)
-    label_map = examples.load_frog_tissues()
 
+@pytest.mark.needs_vtk_version(9, 3, 0)
+def test_contour_labeled_with_invalid_output_style(frog_tissues):
     # Extract surface for each label
-    with pytest.raises(NotImplementedError):
-        label_map.contour_labeled(output_style='selected')
+    with pytest.warns(PyVistaDeprecationWarning):
+        with pytest.raises(NotImplementedError):
+            frog_tissues.contour_labeled(output_style='selected')
 
-    with pytest.raises(ValueError):  # noqa: PT011
-        label_map.contour_labeled(output_style='invalid')
+    with pytest.warns(PyVistaDeprecationWarning):
+        with pytest.raises(ValueError):  # noqa: PT011
+            frog_tissues.contour_labeled(output_style='invalid')
 
 
-@pytest.mark.skipif(not VTK93, reason='At least VTK 9.3 is required')
-def test_contour_labeled_with_scalars():
-    # Load a 3D label map (segmentation of a frog's tissue)
-    # and create a new array with reduced number of labels
-    label_map = examples.load_frog_tissues()
-    label_map['labels'] = label_map['MetaImage'] // 2
+@pytest.mark.needs_vtk_version(9, 3, 0)
+def test_contour_labeled_with_scalars(frog_tissues):
+    # Create a new array with reduced number of labels
+    frog_tissues['labels'] = frog_tissues['MetaImage'] // 2
 
     # Extract surface for each label
-    mesh = label_map.contour_labeled(scalars='labels')
+    with pytest.warns(PyVistaDeprecationWarning):
+        mesh = frog_tissues.contour_labeled(scalars='labels')
 
     assert 'BoundaryLabels' in mesh.cell_data
     assert np.max(mesh['BoundaryLabels'][:, 0]) == 14
 
 
-@pytest.mark.skipif(not VTK93, reason='At least VTK 9.3 is required')
-def test_contour_labeled_with_invalid_scalars():
-    # Load a 3D label map (segmentation of a frog's tissue)
-    label_map = examples.load_frog_tissues()
-
+@pytest.mark.needs_vtk_version(9, 3, 0)
+def test_contour_labeled_with_invalid_scalars(frog_tissues):
     # Nonexistent scalar key
-    with pytest.raises(KeyError):
-        label_map.contour_labeled(scalars='nonexistent_key')
+    with pytest.warns(PyVistaDeprecationWarning):
+        with pytest.raises(KeyError):
+            frog_tissues.contour_labeled(scalars='nonexistent_key')
 
     # Using cell data
-    label_map.cell_data['cell_data'] = np.zeros(label_map.n_cells)
-    with pytest.raises(ValueError, match='Can only process point data'):
-        label_map.contour_labeled(scalars='cell_data')
+    frog_tissues.cell_data['cell_data'] = np.zeros(frog_tissues.n_cells)
+    with pytest.warns(PyVistaDeprecationWarning):
+        with pytest.raises(ValueError, match='Can only process point data'):
+            frog_tissues.contour_labeled(scalars='cell_data')
 
     # When no scalas are given and active scalars are not point data
-    label_map.set_active_scalars('cell_data', preference='cell')
-    with pytest.raises(ValueError, match='active scalars must be point array'):
-        label_map.contour_labeled()
+    frog_tissues.set_active_scalars('cell_data', preference='cell')
+    with pytest.warns(PyVistaDeprecationWarning):
+        with pytest.raises(ValueError, match='active scalars must be point array'):
+            frog_tissues.contour_labeled()
 
 
 @pytest.fixture
diff --git a/tests/core/test_polydata_filters.py b/tests/core/test_polydata_filters.py
index 617f33f6..e67b514a 100644
--- a/tests/core/test_polydata_filters.py
+++ b/tests/core/test_polydata_filters.py
@@ -114,7 +114,7 @@ def frog_tissues_image():
 
 @pytest.fixture
 def frog_tissues_contour(frog_tissues_image):
-    return frog_tissues_image.contour_labeled(smoothing=False)
+    return frog_tissues_image.contour_labels(smoothing=False)
 
 
 @pytest.mark.needs_vtk_version(9, 3, 0)
@@ -244,7 +244,7 @@ def oriented_image():
 
 @pytest.fixture
 def oriented_polydata(oriented_image):
-    oriented_poly = oriented_image.pad_image().contour_labeled(smoothing=False)
+    oriented_poly = oriented_image.pad_image().contour_labels(smoothing=False)
     assert np.allclose(oriented_poly.bounds, oriented_image.points_to_cells().bounds, atol=0.1)
     return oriented_poly
 
@@ -253,7 +253,7 @@ def oriented_polydata(oriented_image):
 def test_voxelize_binary_mask_orientation(oriented_image, oriented_polydata):
     mask = oriented_polydata.voxelize_binary_mask(reference_volume=oriented_image)
     assert mask.bounds == oriented_image.bounds
-    mask_as_surface = mask.pad_image().contour_labeled(smoothing=False)
+    mask_as_surface = mask.pad_image().contour_labels(smoothing=False)
     assert mask_as_surface.bounds == oriented_polydata.bounds