Many issues with contour_labeled, consider deprecating and rewriting as contour_labels

[user27182]

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 #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:

pyvista/pyvista/core/filters/image_data.py

Lines 835 to 836 in a01b7bf

	n_labels : int, optional
	Number of labels to be extracted (all are extracted if None is given).

But the implementation calls GenerateLabels:

pyvista/pyvista/core/filters/image_data.py

Lines 917 to 918 in a01b7bf

	if n_labels is not None:
	alg.GenerateLabels(n_labels, 1, n_labels)

Which is odd because the vtk docs for GenerateLabels 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.:

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:

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:

The reason this isn't currently working is because SetLabel is required, e.g. something like this is needed:

[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.

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:

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 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

PyVista 0.44.dev0

Screenshots

No response