Processing data

Here, we provide functionalities designed specifically for 3D image analysis and processing. From filter pipelines to structure tensor computation and blob detection, qim3d equips you with the tools you need to extract meaningful insights from your data.

qim3d.processing.structure_tensor

structure_tensor(
    volume,
    sigma=1.0,
    rho=6.0,
    base_noise=True,
    smallest=True,
    visualize=False,
    **viz_kwargs,
)

Computes the Structure Tensor for 3D local orientation analysis.

The structure tensor is a matrix-based representation of partial derivatives. It is used to analyze the local orientation of features (like fibers, veins, or layers) in a 3D volume. By examining the eigenvalues and eigenvectors of this tensor, one can determine the dominant direction of structures at every voxel.

This implementation wraps the structure-tensor package, which offers fast, GPU-accelerated computation suitable for large datasets (e.g., micro-CT).

Parameters:

Name	Type	Description	Default
volume	ndarray	The input 3D volume.	required
sigma	float	The noise scale. This defines the standard deviation of the Gaussian kernel used to smooth the image before calculating gradients. Features smaller than this scale are suppressed. Defaults to 1.0.	1.0
rho	float	The integration scale. This defines the size of the neighborhood over which orientation information is averaged (integrated). Larger values yield smoother orientation fields but may blur distinct features. Defaults to 6.0.	6.0
base_noise	bool	If True, adds infinitesimal noise to the volume before processing. This prevents numerical instability (NaNs) in uniform regions where gradients are zero. Defaults to True.	True
smallest	bool	If True, returns only the eigenvector corresponding to the smallest eigenvalue. In gradient-based analysis, the smallest eigenvalue typically corresponds to the direction of least change (i.e., along the fiber/structure axis). If False, returns all eigenvectors. Defaults to True.	True
visualize	bool	If True, immediately displays a visualization of the vector field using qim3d.viz.vectors. Defaults to False.	False
**viz_kwargs	Any	Additional keyword arguments passed to the visualization function (e.g., n_vectors, opacity).	{}

Returns:

Type	Description
(val, vec)	A tuple containing: * val (np.ndarray): An array of shape (3, Z, Y, X) containing the eigenvalues sorted in ascending order. * vec (np.ndarray): * If smallest is True: An array of shape (3, Z, Y, X) representing the vector components (z, y, x) of the primary orientation. * If smallest is False: An array of shape (3, 3, Z, Y, X) containing all three eigenvectors sorted by eigenvalue.

Raises:

Type	Description
ValueError	If the input volume is not 3D.

Example

import qim3d

vol = qim3d.examples.fibers_150x150x150
val, vec = qim3d.processing.structure_tensor(vol, visualize = True, axis = 1)

Runtime and memory usage

Performance computed on Intel(R) Xeon(R) Gold 6226 CPU @ 2.70GHz.

Reference

Jeppesen, N., et al. "Quantifying effects of manufacturing methods on fiber orientation in unidirectional composites using structure tensor analysis." Composites Part A: Applied Science and Manufacturing 149 (2021): 106541. https://doi.org/10.1016/j.compositesa.2021.106541

@article{JEPPESEN2021106541,
title = {Quantifying effects of manufacturing methods on fiber orientation in unidirectional composites using structure tensor analysis},
journal = {Composites Part A: Applied Science and Manufacturing},
volume = {149},
pages = {106541},
year = {2021},
issn = {1359-835X},
doi = {[https://doi.org/10.1016/j.compositesa.2021.106541](https://doi.org/10.1016/j.compositesa.2021.106541)},
url = {[https://www.sciencedirect.com/science/article/pii/S1359835X21002633](https://www.sciencedirect.com/science/article/pii/S1359835X21002633)},
author = {N. Jeppesen and L.P. Mikkelsen and A.B. Dahl and A.N. Christensen and V.A. Dahl}
}

Source code in qim3d/processing/_structure_tensor.py

def structure_tensor(
    volume: np.ndarray,
    sigma: float = 1.0,
    rho: float = 6.0,
    base_noise: bool = True,
    smallest: bool = True,
    visualize: bool = False,
    **viz_kwargs,
) -> tuple[np.ndarray, np.ndarray]:
    """
    Computes the Structure Tensor for 3D local orientation analysis.

    The structure tensor is a matrix-based representation of partial derivatives. It is used to analyze the local orientation of features (like fibers, veins, or layers) in a 3D volume. By examining the eigenvalues and eigenvectors of this tensor, one can determine the dominant direction of structures at every voxel.

    This implementation wraps the [structure-tensor package](https://github.com/Skielex/structure-tensor/), which offers fast, GPU-accelerated computation suitable for large datasets (e.g., micro-CT).

    Args:
        volume (np.ndarray): The input 3D volume.
        sigma (float, optional): The noise scale. This defines the standard deviation of the Gaussian kernel used to smooth the image *before* calculating gradients. Features smaller than this scale are suppressed. Defaults to 1.0.
        rho (float, optional): The integration scale. This defines the size of the neighborhood over which orientation information is averaged (integrated). Larger values yield smoother orientation fields but may blur distinct features. Defaults to 6.0.
        base_noise (bool, optional): If `True`, adds infinitesimal noise to the volume before processing. This prevents numerical instability (NaNs) in uniform regions where gradients are zero. Defaults to `True`.
        smallest (bool, optional): If `True`, returns only the eigenvector corresponding to the **smallest** eigenvalue. In gradient-based analysis, the smallest eigenvalue typically corresponds to the direction of least change (i.e., **along** the fiber/structure axis). If `False`, returns all eigenvectors. Defaults to `True`.
        visualize (bool, optional): If `True`, immediately displays a visualization of the vector field using `qim3d.viz.vectors`. Defaults to `False`.
        **viz_kwargs (Any): Additional keyword arguments passed to the visualization function (e.g., `n_vectors`, `opacity`).

    Returns:
        (val, vec): A tuple containing:
            * **val** (np.ndarray): An array of shape `(3, Z, Y, X)` containing the eigenvalues sorted in ascending order.
            * **vec** (np.ndarray):
                * If `smallest` is `True`: An array of shape `(3, Z, Y, X)` representing the vector components (z, y, x) of the primary orientation.
                * If `smallest` is `False`: An array of shape `(3, 3, Z, Y, X)` containing all three eigenvectors sorted by eigenvalue.

    Raises:
        ValueError: If the input `volume` is not 3D.

    Example:
        ```python
        import qim3d

        vol = qim3d.examples.fibers_150x150x150
        val, vec = qim3d.processing.structure_tensor(vol, visualize = True, axis = 1)
        ```
        ![structure tensor](../../assets/screenshots/structure_tensor_visualization_fibers.gif)

    !!! info "Runtime and memory usage"
        ![structure tensor estimate time and mem](../../assets/screenshots/Structure_tensor_time_mem_estimation.png)
        Performance computed on Intel(R) Xeon(R) Gold 6226 CPU @ 2.70GHz.

    !!! quote "Reference"
        Jeppesen, N., et al. "Quantifying effects of manufacturing methods on fiber orientation in unidirectional composites using structure tensor analysis." Composites Part A: Applied Science and Manufacturing 149 (2021): 106541.
        <https://doi.org/10.1016/j.compositesa.2021.106541>

        ```bibtex
        @article{JEPPESEN2021106541,
        title = {Quantifying effects of manufacturing methods on fiber orientation in unidirectional composites using structure tensor analysis},
        journal = {Composites Part A: Applied Science and Manufacturing},
        volume = {149},
        pages = {106541},
        year = {2021},
        issn = {1359-835X},
        doi = {[https://doi.org/10.1016/j.compositesa.2021.106541](https://doi.org/10.1016/j.compositesa.2021.106541)},
        url = {[https://www.sciencedirect.com/science/article/pii/S1359835X21002633](https://www.sciencedirect.com/science/article/pii/S1359835X21002633)},
        author = {N. Jeppesen and L.P. Mikkelsen and A.B. Dahl and A.N. Christensen and V.A. Dahl}
        }
        ```
    """
    previous_logging_level = logging.getLogger().getEffectiveLevel()
    logging.getLogger().setLevel(logging.CRITICAL)
    import structure_tensor as st

    logging.getLogger().setLevel(previous_logging_level)

    if volume.ndim != 3:
        msg = 'The input volume must be 3D'
        raise ValueError(msg)

    # Ensure volume is a float
    if volume.dtype != np.float32 and volume.dtype != np.float64:
        volume = volume.astype(np.float32)

    if base_noise:
        # Add small noise to the volume
        # FIXME: This is a temporary solution to avoid uniform regions with constant values
        # in the volume, which lead to numerical issues in the structure tensor computation
        vol_noisy = volume + np.random.default_rng(seed=0).uniform(
            0, 1e-10, size=volume.shape
        )

        # Compute the structure tensor (of volume with noise)
        s_vol = st.structure_tensor_3d(vol_noisy, sigma, rho)

    else:
        # Compute the structure tensor (of volume without noise)
        s_vol = st.structure_tensor_3d(volume, sigma, rho)

    # Compute the eigenvalues and eigenvectors of the structure tensor
    full = not smallest
    print(
        f'Computing eigenvalues and eigenvectors of the structure tensor, full = {full}'
    )
    val, vec = st.eig_special_3d(s_vol, full=full, eigenvalue_order='asc')

    if visualize:
        from qim3d.viz import vectors

        display(vectors(volume, vec, **viz_kwargs))

    return val, vec

qim3d.processing.local_thickness

local_thickness(
    image, scale=1, mask=None, visualize=False, **viz_kwargs
)

Computes the local thickness map for a 2D or 3D image.

Local Thickness is a morphometric measure defined as the diameter of the largest sphere that fits entirely within the object boundary and contains the point. Intuitively, it represents the "width" of the structure at any given voxel. It is widely used to analyze pore sizes, trabecular bone thickness, or vessel diameters.

This function wraps the localthickness package, which implements the "Fast Local Thickness" algorithm. Unlike traditional methods that use computationally expensive large structuring elements, this algorithm uses iterative dilation with small kernels, making it significantly faster and feasible for high-resolution 3D microscopy data.

Note: This function requires a binary image. If a grayscale image is provided, it will be automatically binarized using Otsu's thresholding method.

Parameters:

Name	Type	Description	Default
image	ndarray	The input 2D or 3D array. * Binary: Processed directly. * Grayscale: Automatically binarized using Otsu's method (a warning will be logged).	required
scale	float	A downscaling factor to speed up computation. For example, 0.5 downsamples the image by half in each dimension before processing. Defaults to 1 (no scaling).	1
mask	ndarray	A binary mask of the same shape as image. Local thickness will only be computed for regions where the mask is True. Defaults to None.	None
visualize	bool	If True, immediately displays a visualization of the thickness map using qim3d.viz.local_thickness. Defaults to False.	False
**viz_kwargs	Any	Additional keyword arguments passed to the visualization function.	{}

Returns:

Name	Type	Description
local_thickness	ndarray	A NumPy array of the same shape as the input, where each pixel/voxel value represents the local thickness at that point.

Example

import qim3d

vol = qim3d.examples.fly_150x256x256
lt_vol = qim3d.processing.local_thickness(vol, visualize=True, axis=0)

import qim3d

# Generate synthetic collection of blobs
vol, labels = qim3d.generate.volume_collection(num_volumes=15)

# Extract one slice to show that local thickness works on 2D slices too
slice = vol[:,:,50]
lt_blobs = qim3d.processing.local_thickness(slice, visualize=True)

Runtime and memory usage

Performance computed on Intel(R) Xeon(R) Gold 6226 CPU @ 2.70GHz.

Reference

Dahl, V. A., & Dahl, A. B. (2023, June). Fast Local Thickness. 2023 IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops (CVPRW). https://doi.org/10.1109/cvprw59228.2023.00456

@inproceedings{Dahl_2023, title={Fast Local Thickness},
url={[http://dx.doi.org/10.1109/CVPRW59228.2023.00456](http://dx.doi.org/10.1109/CVPRW59228.2023.00456)},
DOI={10.1109/cvprw59228.2023.00456},
booktitle={2023 IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops (CVPRW)},
publisher={IEEE},
author={Dahl, Vedrana Andersen and Dahl, Anders Bjorholm},
year={2023},
month=jun }

Source code in qim3d/processing/_local_thickness.py

def local_thickness(
    image: np.ndarray,
    scale: float = 1,
    mask: Optional[np.ndarray] = None,
    visualize: bool = False,
    **viz_kwargs,
) -> np.ndarray:
    """
    Computes the local thickness map for a 2D or 3D image.

    Local Thickness is a morphometric measure defined as the diameter of the largest sphere that fits entirely within the object boundary and contains the point. Intuitively, it represents the "width" of the structure at any given voxel. It is widely used to analyze pore sizes, trabecular bone thickness, or vessel diameters.

    This function wraps the [localthickness](https://github.com/vedranaa/local-thickness) package, which implements the "Fast Local Thickness" algorithm. Unlike traditional methods that use computationally expensive large structuring elements, this algorithm uses iterative dilation with small kernels, making it significantly faster and feasible for high-resolution 3D microscopy data.

    **Note:** This function requires a **binary** image. If a grayscale image is provided, it will be automatically binarized using Otsu's thresholding method.

    Args:
        image (np.ndarray): The input 2D or 3D array.
            * **Binary:** Processed directly.
            * **Grayscale:** Automatically binarized using Otsu's method (a warning will be logged).
        scale (float, optional): A downscaling factor to speed up computation. For example, `0.5` downsamples the image by half in each dimension before processing. Defaults to 1 (no scaling).
        mask (np.ndarray, optional): A binary mask of the same shape as `image`. Local thickness will only be computed for regions where the mask is `True`. Defaults to `None`.
        visualize (bool, optional): If `True`, immediately displays a visualization of the thickness map using `qim3d.viz.local_thickness`. Defaults to `False`.
        **viz_kwargs (Any): Additional keyword arguments passed to the visualization function.

    Returns:
        local_thickness (np.ndarray):
            A NumPy array of the same shape as the input, where each pixel/voxel value represents the local thickness at that point.

    Example:
        ```python
        import qim3d

        vol = qim3d.examples.fly_150x256x256
        lt_vol = qim3d.processing.local_thickness(vol, visualize=True, axis=0)
        ```
        ![local thickness 3d](../../assets/screenshots/local_thickness_3d.gif)

        ```python
        import qim3d

        # Generate synthetic collection of blobs
        vol, labels = qim3d.generate.volume_collection(num_volumes=15)

        # Extract one slice to show that local thickness works on 2D slices too
        slice = vol[:,:,50]
        lt_blobs = qim3d.processing.local_thickness(slice, visualize=True)

        ```
        ![local thickness 2d](../../assets/screenshots/local_thickness_2d.png)

    !!! info "Runtime and memory usage"
        ![local thickness estimate time and mem](../../assets/screenshots/Local_thickness_time_mem_estimation.png)
        Performance computed on Intel(R) Xeon(R) Gold 6226 CPU @ 2.70GHz.

    !!! quote "Reference"
        Dahl, V. A., & Dahl, A. B. (2023, June). Fast Local Thickness. 2023 IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops (CVPRW).
        <https://doi.org/10.1109/cvprw59228.2023.00456>

        ```bibtex
        @inproceedings{Dahl_2023, title={Fast Local Thickness},
        url={[http://dx.doi.org/10.1109/CVPRW59228.2023.00456](http://dx.doi.org/10.1109/CVPRW59228.2023.00456)},
        DOI={10.1109/cvprw59228.2023.00456},
        booktitle={2023 IEEE/CVF Conference on Computer Vision and Pattern Recognition Workshops (CVPRW)},
        publisher={IEEE},
        author={Dahl, Vedrana Andersen and Dahl, Anders Bjorholm},
        year={2023},
        month=jun }
        ```
    """
    import localthickness as lt
    from skimage.filters import threshold_otsu

    # Check if input is binary
    if np.unique(image).size > 2:
        # If not, binarize it using Otsu's method, log the threshold and compute the local thickness
        threshold = threshold_otsu(image=image)
        log.warning(
            f"Input image is not binary. It will be binarized using Otsu's method with threshold: {threshold}"
        )
        local_thickness = lt.local_thickness(image > threshold, scale=scale, mask=mask)
    else:
        # If it is binary, compute the local thickness directly
        local_thickness = lt.local_thickness(image, scale=scale, mask=mask)

    # Visualize the local thickness if requested
    if visualize:
        display(qim3d.viz.local_thickness(image, local_thickness, **viz_kwargs))

    return local_thickness

qim3d.processing.get_lines

get_lines(segmentations)

Extracts 1D line coordinates from 2D binary segmentation masks.

This utility function is designed to work with the output of qim3d.processing.segment_layers. It converts the binary masks (which split the image into "above" and "below" a layer) into a 1D array of height indices. This allows for easy plotting of the layer boundary using Matplotlib.

Parameters:

Name	Type	Description	Default
segmentations	list[ndarray]	A list of 2D binary arrays, typically returned by segment_layers. Each array should contain two classes (0 and 1) separated by a boundary.	required

Returns:

Name	Type	Description
segmentation_lines	list[ndarray]	A list of 1D arrays. Each array contains the vertical index (y-coordinate) of the layer boundary for every horizontal position (x-coordinate).

Example

# Assuming 'layers' is the output from segment_layers
lines = qim3d.processing.get_lines(layers)

# Plotting the first layer
plt.plot(lines[0], color='red')

Source code in qim3d/processing/_layers.py

def get_lines(segmentations: list[np.ndarray]) -> list:
    """
    Extracts 1D line coordinates from 2D binary segmentation masks.

    This utility function is designed to work with the output of `qim3d.processing.segment_layers`. It converts the binary masks (which split the image into "above" and "below" a layer) into a 1D array of height indices. This allows for easy plotting of the layer boundary using Matplotlib.

    Args:
        segmentations (list[np.ndarray]): A list of 2D binary arrays, typically returned by `segment_layers`. Each array should contain two classes (0 and 1) separated by a boundary.

    Returns:
        segmentation_lines (list[np.ndarray]):
            A list of 1D arrays. Each array contains the vertical index (y-coordinate) of the layer boundary for every horizontal position (x-coordinate).

    Example:
        ```python
        # Assuming 'layers' is the output from segment_layers
        lines = qim3d.processing.get_lines(layers)

        # Plotting the first layer
        plt.plot(lines[0], color='red')
        ```
    """
    segmentation_lines = [np.argmin(s, axis=0) - 0.5 for s in segmentations]
    return segmentation_lines

qim3d.processing.segment_layers

segment_layers(
    data,
    inverted=False,
    n_layers=1,
    delta=1,
    min_margin=10,
    max_margin=None,
    wrap=False,
)

Segments distinct layers in 2D or 3D data using graph-based optimal surface detection.

This function identifies continuous surfaces (layers) within the volume that minimize a specific cost function. It uses a Graph Cut (Max-Flow/Min-Cut) algorithm to find the global optimum. This is particularly useful for detecting boundaries in geological strata, biological tissues (e.g., retinal layers), or material interfaces.

It acts as a wrapper around the slgbuilder library.

Parameters:

Name	Type	Description	Default
data	ndarray	The 2D or 3D input array. The algorithm assumes layers are stacked along the first dimension (axis 0).	required
inverted	bool	If True, inverts the intensity of the image before processing. Use this if the boundaries you are looking for are dark instead of bright (or vice-versa, depending on the cost function logic). Defaults to False.	False
n_layers	int	The number of surfaces/boundaries to detect. Defaults to 1.	1
delta	float	The smoothness penalty. Higher values enforce smoother (stiffer) layer boundaries, resisting sudden changes in height. Defaults to 1.	1
min_margin	int or None	The minimum vertical distance (in pixels/voxels) allowed between consecutive layers. Used to prevent surfaces from crossing or collapsing onto each other. Defaults to 10.	10
max_margin	int or None	The maximum vertical distance allowed between consecutive layers. Defaults to None.	None
wrap	bool	If True, enforces a periodic boundary condition where the start and end of the layer (along the width) must connect. Useful for cylindrical or unwrapped data. Defaults to False.	False

Returns:

Name	Type	Description
segmentations	list[ndarray]	A list of binary masks (0s and 1s), where each mask represents the region defined by a detected layer. The list length equals n_layers.

Raises:

Type	Description
TypeError	If data is not a numpy array or n_layers is not an integer.
ValueError	If n_layers is less than 1 or delta is non-positive.

Example

import qim3d
import matplotlib.pyplot as plt

# Load data (using a 2D slice for this example)
# In this image, we want to find 2 distinct boundaries
layers_image = qim3d.io.load('layers3d.tif')[:,:,0]

# Segment the layers
layers = qim3d.processing.segment_layers(layers_image, n_layers=2)

# Extract the line coordinates for plotting
layer_lines = qim3d.processing.get_lines(layers)

# Visualize
plt.imshow(layers_image, cmap='gray')
plt.axis('off')
for layer_line in layer_lines:
    plt.plot(layer_line, linewidth=3, label='Detected Layer')
plt.legend()
plt.show()

Source code in qim3d/processing/_layers.py

def segment_layers(
    data: np.ndarray,
    inverted: bool = False,
    n_layers: int = 1,
    delta: float = 1,
    min_margin: int = 10,
    max_margin: int = None,
    wrap: bool = False,
) -> list:
    """
    Segments distinct layers in 2D or 3D data using graph-based optimal surface detection.

    This function identifies continuous surfaces (layers) within the volume that minimize a specific cost function. It uses a **Graph Cut (Max-Flow/Min-Cut)** algorithm to find the global optimum. This is particularly useful for detecting boundaries in geological strata, biological tissues (e.g., retinal layers), or material interfaces.

    It acts as a wrapper around the [slgbuilder](https://github.com/Skielex/slgbuilder) library.

    Args:
        data (np.ndarray): The 2D or 3D input array. The algorithm assumes layers are stacked along the first dimension (axis 0).
        inverted (bool, optional): If `True`, inverts the intensity of the image before processing. Use this if the boundaries you are looking for are dark instead of bright (or vice-versa, depending on the cost function logic). Defaults to `False`.
        n_layers (int, optional): The number of surfaces/boundaries to detect. Defaults to 1.
        delta (float, optional): The smoothness penalty. Higher values enforce smoother (stiffer) layer boundaries, resisting sudden changes in height. Defaults to 1.
        min_margin (int or None, optional): The minimum vertical distance (in pixels/voxels) allowed between consecutive layers. Used to prevent surfaces from crossing or collapsing onto each other. Defaults to 10.
        max_margin (int or None, optional): The maximum vertical distance allowed between consecutive layers. Defaults to `None`.
        wrap (bool, optional): If `True`, enforces a periodic boundary condition where the start and end of the layer (along the width) must connect. Useful for cylindrical or unwrapped data. Defaults to `False`.

    Returns:
        segmentations (list[np.ndarray]):
            A list of binary masks (0s and 1s), where each mask represents the region defined by a detected layer. The list length equals `n_layers`.

    Raises:
        TypeError: If `data` is not a numpy array or `n_layers` is not an integer.
        ValueError: If `n_layers` is less than 1 or `delta` is non-positive.

    Example:
        ```python
        import qim3d
        import matplotlib.pyplot as plt

        # Load data (using a 2D slice for this example)
        # In this image, we want to find 2 distinct boundaries
        layers_image = qim3d.io.load('layers3d.tif')[:,:,0]

        # Segment the layers
        layers = qim3d.processing.segment_layers(layers_image, n_layers=2)

        # Extract the line coordinates for plotting
        layer_lines = qim3d.processing.get_lines(layers)

        # Visualize
        plt.imshow(layers_image, cmap='gray')
        plt.axis('off')
        for layer_line in layer_lines:
            plt.plot(layer_line, linewidth=3, label='Detected Layer')
        plt.legend()
        plt.show()
        ```
        ![layer_segmentation](../../assets/screenshots/layers.png)
        ![layer_segmentation](../../assets/screenshots/segmented_layers.png)
    """
    if isinstance(data, np.ndarray):
        data = data.astype(np.int32)
        if inverted:
            data = ~data
    else:
        raise TypeError(
            f'Data has to be type np.ndarray. Your data is of type {type(data)}'
        )

    helper = MaxflowBuilder()
    if not isinstance(n_layers, int):
        raise TypeError(
            f'Number of layers has to be positive integer. You passed {type(n_layers)}'
        )

    if n_layers == 1:
        layer = GraphObject(data)
        helper.add_object(layer)
    elif n_layers > 1:
        layers = [GraphObject(data) for _ in range(n_layers)]
        helper.add_objects(layers)
        for i in range(len(layers) - 1):
            helper.add_layered_containment(
                layers[i], layers[i + 1], min_margin=min_margin, max_margin=max_margin
            )

    else:
        raise ValueError(
            f'Number of layers has to be positive integer. You passed {n_layers}'
        )

    helper.add_layered_boundary_cost()

    if delta > 1:
        delta = int(delta)
    elif delta <= 0:
        raise ValueError(f'Delta has to be positive number. You passed {delta}')
    helper.add_layered_smoothness(delta=delta, wrap=bool(wrap))
    helper.solve()
    if n_layers == 1:
        segmentations = [helper.what_segments(layer)]
    else:
        segmentations = [helper.what_segments(l).astype(np.int32) for l in layers]

    return segmentations