Skip to content

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

qim3d.processing.blob_detection 

blob_detection(vol, background='dark', min_sigma=1, max_sigma=50, sigma_ratio=1.6, threshold=0.5, overlap=0.5, threshold_rel=None, exclude_border=False)

Extract blobs from a volume using Difference of Gaussian (DoG) method, and retrieve a binary volume with the blobs marked as True

Parameters:

Name Type Description Default
vol ndarray

The volume to detect blobs in

 required
background str

'dark' if background is darker than the blobs, 'bright' if background is lighter than the blobs

 'dark'
min_sigma float

The minimum standard deviation for Gaussian kernel

 1
max_sigma float

The maximum standard deviation for Gaussian kernel

 50
sigma_ratio float

The ratio between the standard deviation of Gaussian Kernels

 1.6
threshold float

The absolute lower bound for scale space maxima. Reduce this to detect blobs with lower intensities

 0.5
overlap float

The fraction of area of two blobs that overlap

 0.5
threshold_rel float

The relative lower bound for scale space maxima

 None
exclude_border bool

If True, exclude blobs that are too close to the border of the image

 False

Returns:

Name Type Description
blobs ndarray

The blobs found in the volume as (p, r, c, radius)
binary_volume ndarray

A binary volume with the blobs marked as True
Example 

import qim3d

# Get data
vol = qim3d.examples.cement_128x128x128
vol_blurred = qim3d.processing.gaussian(vol, sigma=2)

# Detect blobs, and get binary mask
blobs, mask = qim3d.processing.blob_detection(
    vol_blurred,
    min_sigma=1,
    max_sigma=8,
    threshold=0.001,
    overlap=0.1,
    background="bright"
    )

# Visualize detected blobs
qim3d.viz.circles(blobs, vol, alpha=0.8, color='blue')
blob detection

# Visualize binary mask
qim3d.viz.slicer(mask)
blob detection

Source code in qim3d/processing/detection.py 
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
def blob_detection(
    vol: np.ndarray,
    background: str = "dark",
    min_sigma: float = 1,
    max_sigma: float = 50,
    sigma_ratio: float = 1.6,
    threshold: float = 0.5,
    overlap: float = 0.5,
    threshold_rel: float = None,
    exclude_border: bool = False,
) -> np.ndarray:
    """
    Extract blobs from a volume using Difference of Gaussian (DoG) method, and retrieve a binary volume with the blobs marked as True

    Args:
        vol: The volume to detect blobs in
        background: 'dark' if background is darker than the blobs, 'bright' if background is lighter than the blobs
        min_sigma: The minimum standard deviation for Gaussian kernel
        max_sigma: The maximum standard deviation for Gaussian kernel
        sigma_ratio: The ratio between the standard deviation of Gaussian Kernels
        threshold: The absolute lower bound for scale space maxima. Reduce this to detect blobs with lower intensities
        overlap: The fraction of area of two blobs that overlap
        threshold_rel: The relative lower bound for scale space maxima
        exclude_border: If True, exclude blobs that are too close to the border of the image

    Returns:
        blobs: The blobs found in the volume as (p, r, c, radius)
        binary_volume: A binary volume with the blobs marked as True

    Example:
            ```python
            import qim3d

            # Get data
            vol = qim3d.examples.cement_128x128x128
            vol_blurred = qim3d.processing.gaussian(vol, sigma=2)

            # Detect blobs, and get binary mask
            blobs, mask = qim3d.processing.blob_detection(
                vol_blurred,
                min_sigma=1,
                max_sigma=8,
                threshold=0.001,
                overlap=0.1,
                background="bright"
                )

            # Visualize detected blobs
            qim3d.viz.circles(blobs, vol, alpha=0.8, color='blue')
            ```
            ![blob detection](assets/screenshots/blob_detection.gif)

            ```python
            # Visualize binary mask
            qim3d.viz.slicer(mask)
            ```
            ![blob detection](assets/screenshots/blob_get_mask.gif)
    """
    from skimage.feature import blob_dog

    if background == "bright":
        log.info("Bright background selected, volume will be inverted.")
        vol = np.invert(vol)

    blobs = blob_dog(
        vol,
        min_sigma=min_sigma,
        max_sigma=max_sigma,
        sigma_ratio=sigma_ratio,
        threshold=threshold,
        overlap=overlap,
        threshold_rel=threshold_rel,
        exclude_border=exclude_border,
    )

    # Change sigma to radius
    blobs[:, 3] = blobs[:, 3] * np.sqrt(3)

    # Create binary mask of detected blobs
    vol_shape = vol.shape
    binary_volume = np.zeros(vol_shape, dtype=bool)

    for z, y, x, radius in blobs:
        # Calculate the bounding box around the blob
        z_start = max(0, int(z - radius))
        z_end = min(vol_shape[0], int(z + radius) + 1)
        y_start = max(0, int(y - radius))
        y_end = min(vol_shape[1], int(y + radius) + 1)
        x_start = max(0, int(x - radius))
        x_end = min(vol_shape[2], int(x + radius) + 1)

        z_indices, y_indices, x_indices = np.indices(
            (z_end - z_start, y_end - y_start, x_end - x_start)
        )
        z_indices += z_start
        y_indices += y_start
        x_indices += x_start

        # Calculate distances from the center of the blob to voxels within the bounding box
        dist = np.sqrt(
            (x_indices - x) ** 2 + (y_indices - y) ** 2 + (z_indices - z) ** 2
        )

        binary_volume[z_start:z_end, y_start:y_end, x_start:x_end][
            dist <= radius
        ] = True

    return blobs, binary_volume

qim3d.processing.structure_tensor 

structure_tensor(vol, sigma=1.0, rho=6.0, base_noise=True, full=False, visualize=False, **viz_kwargs)

Wrapper for the 3D structure tensor implementation from the structure_tensor package.

The structure tensor algorithm is a method for analyzing the orientation of fiber-like structures in 3D images. The core of the algorithm involves computing a 3-by-3 matrix at each point in a volume, capturing the local orientation. This matrix, known as the structure tensor, is derived from the gradients of the image intensities and integrates neighborhood information using Gaussian kernels.

The implementation here used allows for fast and efficient computation using GPU-based processing, making it suitable for large datasets. This efficiency is particularly advantageous for high-resolution imaging techniques like X-ray computed microtomography (μCT).

Parameters:

Name Type Description Default
vol ndarray

3D NumPy array representing the volume.

 required
sigma float

A noise scale, structures smaller than sigma will be removed by smoothing.

 1.0
rho float

An integration scale giving the size over the neighborhood in which the orientation is to be analysed.

 6.0
base_noise bool

A flag indicating whether to add a small noise to the volume. Default is True.

 True
full bool

A flag indicating that all three eigenvalues should be returned. Default is False.

 False
visualize bool

Whether to visualize the structure tensor. Default is False.

 False
**viz_kwargs

Additional keyword arguments for passed to qim3d.viz.vectors. Only used if visualize=True.

 {}

Raises:

Type Description
ValueError

If the input volume is not 3D.

Returns:

Name Type Description
val ndarray

An array with shape (3, *vol.shape) containing the eigenvalues of the structure tensor.
vec ndarray

An array with shape (3, *vol.shape) if full is False, otherwise (3, 3, *vol.shape) containing eigenvectors.
Example 

import qim3d

vol = qim3d.examples.NT_128x128x128
val, vec = qim3d.processing.structure_tensor(vol, visualize = True, axis = 2)
structure tensor

Runtime and memory usage of the structure tensor method for different volume sizes

structure tensor estimate time and mem

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},
url = {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 
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
def structure_tensor(
    vol: np.ndarray,
    sigma: float = 1.0,
    rho: float = 6.0,
    base_noise: bool = True,
    full: bool = False,
    visualize=False,
    **viz_kwargs
) -> Tuple[np.ndarray, np.ndarray]:
    """Wrapper for the 3D structure tensor implementation from the [structure_tensor package](https://github.com/Skielex/structure-tensor/).


    The structure tensor algorithm is a method for analyzing the orientation of fiber-like structures in 3D images.
    The core of the algorithm involves computing a 3-by-3 matrix at each point in a volume, capturing the local orientation. This matrix, known as the structure tensor, is derived from the gradients of the image intensities and integrates neighborhood information using Gaussian kernels.

    The implementation here used allows for fast and efficient computation using GPU-based processing, making it suitable for large datasets.
    This efficiency is particularly advantageous for high-resolution imaging techniques like X-ray computed microtomography (μCT).

    Args:
        vol (np.ndarray): 3D NumPy array representing the volume.
        sigma (float, optional): A noise scale, structures smaller than sigma will be removed by smoothing.
        rho (float, optional): An integration scale giving the size over the neighborhood in which the orientation is to be analysed.
        base_noise (bool, optional): A flag indicating whether to add a small noise to the volume. Default is True.
        full (bool, optional): A flag indicating that all three eigenvalues should be returned. Default is False.
        visualize (bool, optional): Whether to visualize the structure tensor. Default is False.
        **viz_kwargs: Additional keyword arguments for passed to `qim3d.viz.vectors`. Only used if `visualize=True`.

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

    Returns:
        val: An array with shape `(3, *vol.shape)` containing the eigenvalues of the structure tensor.
        vec: An array with shape `(3, *vol.shape)` if `full` is `False`, otherwise `(3, 3, *vol.shape)` containing eigenvectors.

    Example:
        ```python
        import qim3d

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


    !!! info "Runtime and memory usage of the structure tensor method for different volume sizes"
        ![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},
        url = {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 vol.ndim != 3:
        raise ValueError("The input volume must be 3D")

    # Ensure volume is a float
    if vol.dtype != np.float32 and vol.dtype != np.float64:
        vol = vol.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 = vol + np.random.default_rng(seed=0).uniform(
            0, 1e-10, size=vol.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(vol, sigma, rho)

    # Compute the eigenvalues and eigenvectors of the structure tensor
    val, vec = st.eig_special_3d(s_vol, full=full)

    if visualize:
        from qim3d.viz.structure_tensor import vectors

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

    return val, vec

qim3d.processing.local_thickness 

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

Wrapper for the local thickness function from the local thickness package

The "Fast Local Thickness" by Vedrana Andersen Dahl and Anders Bjorholm Dahl from the Technical University of Denmark is a efficient algorithm for computing local thickness in 2D and 3D images. Their method significantly reduces computation time compared to traditional algorithms by utilizing iterative dilation with small structuring elements, rather than the large ones typically used. This approach allows the local thickness to be determined much faster, making it feasible for high-resolution volumetric data that are common in contemporary 3D microscopy.

Testing against conventional methods and other Python-based tools like PoreSpy shows that the new algorithm is both accurate and faster, offering significant improvements in processing time for large datasets.

Parameters:

Name Type Description Default
image ndarray

2D or 3D NumPy array representing the image/volume. If binary, it will be passed directly to the local thickness function. If grayscale, it will be binarized using Otsu's method.

 required
scale float

Downscaling factor, e.g. 0.5 for halving each dim of the image. Default is 1.

 1
mask ndarray

binary mask of the same size of the image defining parts of the image to be included in the computation of the local thickness. Default is None.

 None
visualize bool

Whether to visualize the local thickness. Default is False.

 False
**viz_kwargs

Additional keyword arguments passed to qim3d.viz.local_thickness. Only used if visualize=True.

 {}

Returns:

Name Type Description
local_thickness ndarray

2D or 3D NumPy array representing the local thickness of the input image/volume.
Example 

import qim3d

fly = qim3d.examples.fly_150x256x256 # 3D volume
lt_fly = qim3d.processing.local_thickness(fly, visualize=True, axis=0)
local thickness 3d

import qim3d

# Generate synthetic collection of blobs
num_objects = 15
synthetic_collection, labels = qim3d.generate.collection(num_objects = num_objects)

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

Runtime and memory usage of the local thickness method for different volume sizes

local thickness estimate time and mem

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},
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 
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
def local_thickness(
    image: np.ndarray,
    scale: float = 1,
    mask: Optional[np.ndarray] = None,
    visualize=False,
    **viz_kwargs
) -> np.ndarray:
    """Wrapper for the local thickness function from the [local thickness package](https://github.com/vedranaa/local-thickness)

    The "Fast Local Thickness" by Vedrana Andersen Dahl and Anders Bjorholm Dahl from the Technical University of Denmark is a efficient algorithm for computing local thickness in 2D and 3D images.
    Their method significantly reduces computation time compared to traditional algorithms by utilizing iterative dilation with small structuring elements, rather than the large ones typically used.
    This approach allows the local thickness to be determined much faster, making it feasible for high-resolution volumetric data that are common in contemporary 3D microscopy.

    Testing against conventional methods and other Python-based tools like PoreSpy shows that the new algorithm is both accurate and faster, offering significant improvements in processing time for large datasets.


    Args:
        image (np.ndarray): 2D or 3D NumPy array representing the image/volume.
            If binary, it will be passed directly to the local thickness function.
            If grayscale, it will be binarized using Otsu's method.
        scale (float, optional): Downscaling factor, e.g. 0.5 for halving each dim of the image.
            Default is 1.
        mask (np.ndarray, optional): binary mask of the same size of the image defining parts of the
            image to be included in the computation of the local thickness. Default is None.
        visualize (bool, optional): Whether to visualize the local thickness. Default is False.
        **viz_kwargs: Additional keyword arguments passed to `qim3d.viz.local_thickness`. Only used if `visualize=True`.

    Returns:
        local_thickness (np.ndarray): 2D or 3D NumPy array representing the local thickness of the input image/volume.

    Example:
        ```python
        import qim3d

        fly = qim3d.examples.fly_150x256x256 # 3D volume
        lt_fly = qim3d.processing.local_thickness(fly, visualize=True, axis=0)
        ```
        ![local thickness 3d](assets/screenshots/local_thickness_3d.gif)

        ```python
        import qim3d

        # Generate synthetic collection of blobs
        num_objects = 15
        synthetic_collection, labels = qim3d.generate.collection(num_objects = num_objects)

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

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

    !!! info "Runtime and memory usage of the local thickness method for different volume sizes"
        ![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},
        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(
            "Input image is not binary. It will be binarized using Otsu's method with threshold: {}".format(
                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_3d_cc 

get_3d_cc(image)

Returns an object (CC) containing the connected components of the input volume. Use plot_cc to visualize the connected components.

Parameters:

Name Type Description Default
image ndarray

An array-like object to be labeled. Any non-zero values in input are counted as features and zero values are considered the background.

 required

Returns:

Name Type Description
CC CC

A ConnectedComponents object containing the connected components and the number of connected components.
Example 
import qim3d
vol = qim3d.examples.cement_128x128x128[50:150]<60
cc = qim3d.processing.get_3d_cc(vol)
Source code in qim3d/processing/cc.py 
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
def get_3d_cc(image: np.ndarray) -> CC:
    """ Returns an object (CC) containing the connected components of the input volume. Use plot_cc to visualize the connected components.

    Args:
        image (np.ndarray): An array-like object to be labeled. Any non-zero values in `input` are
            counted as features and zero values are considered the background.

    Returns:
        CC: A ConnectedComponents object containing the connected components and the number of connected components.

    Example:
        ```python
        import qim3d
        vol = qim3d.examples.cement_128x128x128[50:150]<60
        cc = qim3d.processing.get_3d_cc(vol)
        ```
    """
    connected_components, num_connected_components = label(image)
    log.info(f"Total number of connected components found: {num_connected_components}")
    return CC(connected_components, num_connected_components)

qim3d.processing.gaussian 

gaussian(vol, dask=False, chunks='auto', *args, **kwargs)

Applies a Gaussian filter to the input volume using scipy.ndimage.gaussian_filter or dask_image.ndfilters.gaussian_filter.

Parameters:

Name Type Description Default
vol

The input image or volume.

 required
dask

Whether to use Dask for the Gaussian filter.

 False
chunks

Defines how to divide the array into blocks when using Dask. Can be an integer, tuple, size in bytes, or "auto" for automatic sizing.

 'auto'
*args

Additional positional arguments for the Gaussian filter.

 ()
**kwargs

Additional keyword arguments for the Gaussian filter.

 {}

Returns:

Type Description

The filtered image or volume.
Source code in qim3d/processing/filters.py 
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
def gaussian(vol, dask=False, chunks='auto', *args, **kwargs):
    """
    Applies a Gaussian filter to the input volume using scipy.ndimage.gaussian_filter or dask_image.ndfilters.gaussian_filter.

    Args:
        vol: The input image or volume.
        dask: Whether to use Dask for the Gaussian filter.
        chunks: Defines how to divide the array into blocks when using Dask. Can be an integer, tuple, size in bytes, or "auto" for automatic sizing.
        *args: Additional positional arguments for the Gaussian filter.
        **kwargs: Additional keyword arguments for the Gaussian filter.

    Returns:
        The filtered image or volume.
    """

    if dask:
        if not isinstance(vol, da.Array):
            vol = da.from_array(vol, chunks=chunks)
        dask_vol = dask_ndfilters.gaussian_filter(vol, *args, **kwargs)
        res = dask_vol.compute()
        return res
    else:
        res = ndimage.gaussian_filter(vol, *args, **kwargs)
        return res

qim3d.processing.median 

median(vol, dask=False, chunks='auto', **kwargs)

Applies a median filter to the input volume using scipy.ndimage.median_filter or dask_image.ndfilters.median_filter.

Parameters:

Name Type Description Default
vol

The input image or volume.

 required
dask

Whether to use Dask for the median filter.

 False
chunks

Defines how to divide the array into blocks when using Dask. Can be an integer, tuple, size in bytes, or "auto" for automatic sizing.

 'auto'
**kwargs

Additional keyword arguments for the median filter.

 {}

Returns:

Type Description

The filtered image or volume.
Source code in qim3d/processing/filters.py 
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
def median(vol, dask=False, chunks='auto', **kwargs):
    """
    Applies a median filter to the input volume using scipy.ndimage.median_filter or dask_image.ndfilters.median_filter.

    Args:
        vol: The input image or volume.
        dask: Whether to use Dask for the median filter.
        chunks: Defines how to divide the array into blocks when using Dask. Can be an integer, tuple, size in bytes, or "auto" for automatic sizing.
        **kwargs: Additional keyword arguments for the median filter.

    Returns:
        The filtered image or volume.
    """
    if dask:
        if not isinstance(vol, da.Array):
            vol = da.from_array(vol, chunks=chunks)
        dask_vol = dask_ndfilters.median_filter(vol, **kwargs)
        res = dask_vol.compute()
        return res
    else:
        res = ndimage.median_filter(vol, **kwargs)
        return res

qim3d.processing.maximum 

maximum(vol, dask=False, chunks='auto', **kwargs)

Applies a maximum filter to the input volume using scipy.ndimage.maximum_filter or dask_image.ndfilters.maximum_filter.

Parameters:

Name Type Description Default
vol

The input image or volume.

 required
dask

Whether to use Dask for the maximum filter.

 False
chunks

Defines how to divide the array into blocks when using Dask. Can be an integer, tuple, size in bytes, or "auto" for automatic sizing.

 'auto'
**kwargs

Additional keyword arguments for the maximum filter.

 {}

Returns:

Type Description

The filtered image or volume.
Source code in qim3d/processing/filters.py 
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
def maximum(vol, dask=False, chunks='auto', **kwargs):
    """
    Applies a maximum filter to the input volume using scipy.ndimage.maximum_filter or dask_image.ndfilters.maximum_filter.

    Args:
        vol: The input image or volume.
        dask: Whether to use Dask for the maximum filter.
        chunks: Defines how to divide the array into blocks when using Dask. Can be an integer, tuple, size in bytes, or "auto" for automatic sizing.
        **kwargs: Additional keyword arguments for the maximum filter.

    Returns:
        The filtered image or volume.
    """
    if dask:
        if not isinstance(vol, da.Array):
            vol = da.from_array(vol, chunks=chunks)
        dask_vol = dask_ndfilters.maximum_filter(vol, **kwargs)
        res = dask_vol.compute()
        return res
    else:
        res = ndimage.maximum_filter(vol, **kwargs)
        return res

qim3d.processing.minimum 

minimum(vol, dask=False, chunks='auto', **kwargs)

Applies a minimum filter to the input volume using scipy.ndimage.minimum_filter or dask_image.ndfilters.minimum_filter.

Parameters:

Name Type Description Default
vol

The input image or volume.

 required
dask

Whether to use Dask for the minimum filter.

 False
chunks

Defines how to divide the array into blocks when using Dask. Can be an integer, tuple, size in bytes, or "auto" for automatic sizing.

 'auto'
**kwargs

Additional keyword arguments for the minimum filter.

 {}

Returns:

Type Description

The filtered image or volume.
Source code in qim3d/processing/filters.py 
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
def minimum(vol, dask=False, chunks='auto', **kwargs):
    """
    Applies a minimum filter to the input volume using scipy.ndimage.minimum_filter or dask_image.ndfilters.minimum_filter.

    Args:
        vol: The input image or volume.
        dask: Whether to use Dask for the minimum filter.
        chunks: Defines how to divide the array into blocks when using Dask. Can be an integer, tuple, size in bytes, or "auto" for automatic sizing.
        **kwargs: Additional keyword arguments for the minimum filter.

    Returns:
        The filtered image or volume.
    """
    if dask:
        if not isinstance(vol, da.Array):
            vol = da.from_array(vol, chunks=chunks)
        dask_vol = dask_ndfilters.minimum_filter(vol, **kwargs)
        res = dask_vol.compute()
        return res
    else:
        res = ndimage.minimum_filter(vol, **kwargs)
        return res

qim3d.processing.tophat 

tophat(vol, dask=False, chunks='auto', **kwargs)

Remove background from the volume.

Parameters:

Name Type Description Default
vol

The volume to remove background from.

 required
radius

The radius of the structuring element (default: 3).

 required
background

Color of the background, 'dark' or 'bright' (default: 'dark'). If 'bright', volume will be inverted.

 required
dask

Whether to use Dask for the tophat filter (not supported, will default to SciPy).

 False
chunks

Defines how to divide the array into blocks when using Dask. Can be an integer, tuple, size in bytes, or "auto" for automatic sizing.

 'auto'
**kwargs

Additional keyword arguments.

 {}

Returns:

Name Type Description
vol

The volume with background removed.
Source code in qim3d/processing/filters.py 
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
def tophat(vol, dask=False, chunks='auto', **kwargs):
    """
    Remove background from the volume.

    Args:
        vol: The volume to remove background from.
        radius: The radius of the structuring element (default: 3).
        background: Color of the background, 'dark' or 'bright' (default: 'dark'). If 'bright', volume will be inverted.
        dask: Whether to use Dask for the tophat filter (not supported, will default to SciPy).
        chunks: Defines how to divide the array into blocks when using Dask. Can be an integer, tuple, size in bytes, or "auto" for automatic sizing.
        **kwargs: Additional keyword arguments.

    Returns:
        vol: The volume with background removed.
    """

    radius = kwargs["radius"] if "radius" in kwargs else 3
    background = kwargs["background"] if "background" in kwargs else "dark"

    if dask:
        log.info("Dask not supported for tophat filter, switching to scipy.")

    if background == "bright":
        log.info("Bright background selected, volume will be temporarily inverted when applying white_tophat")
        vol = np.invert(vol)

    selem = morphology.ball(radius)
    vol = vol - morphology.white_tophat(vol, selem)

    if background == "bright":
        vol = np.invert(vol)

    return vol

qim3d.processing.get_lines 

get_lines(segmentations)

Expects list of arrays where each array is 2D segmentation with only 2 classes. This function gets the border between those two so it could be plotted. Used with qim3d.processing.segment_layers

Parameters:

Name Type Description Default
segmentations list | ndarray

list of arrays where each array is 2D segmentation with only 2 classes

 required

Returns:

Name Type Description
segmentation_lines list

list of 1D numpy arrays
Source code in qim3d/processing/layers2d.py 
85
86
87
88
89
90
91
92
93
94
95
96
97
def get_lines(segmentations:list|np.ndarray) -> list:
    """
    Expects list of arrays where each array is 2D segmentation with only 2 classes. This function gets the border between those two
    so it could be plotted. Used with qim3d.processing.segment_layers

    Args:
        segmentations: list of arrays where each array is 2D segmentation with only 2 classes

    Returns:
        segmentation_lines: list of 1D numpy arrays
    """
    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)

Works on 2D and 3D data. Light one function wrapper around slgbuilder https://github.com/Skielex/slgbuilder to do layer segmentation Now uses only MaxflowBuilder for solving.

Parameters:

Name Type Description Default
data ndarray

2D or 3D array on which it will be computed

 required
inverted bool

if True, it will invert the brightness of the image

 False
n_layers int

How many layers are we looking for (result in a layer and background)

 1
delta float

Smoothness parameter

 1
min_margin int

If we want more layers, we have to have a margin otherwise they are all going to be exactly the same

 10
max_margin int

Maximum margin between layers

 None
wrap bool

If True, starting and ending point of the border between layers are at the same level

 False

Returns:

Name Type Description
segmentations

list of numpy arrays, even if n_layers == 1, each array is only 0s and 1s, 1s segmenting this specific layer

Raises:

Type Description
TypeError

If Data is not np.array, if n_layers is not integer.
ValueError

If n_layers is less than 1, if delta is negative or zero
Example

Example is only shown on 2D image, but segment_layers can also take 3D structures. 

import qim3d

layers_image = qim3d.io.load('layers3d.tif')[:,:,0]
layers = qim3d.processing.segment_layers(layers_image, n_layers = 2)
layer_lines = qim3d.processing.get_lines(layers)

from matplotlib import pyplot as plt

plt.imshow(layers_image, cmap='gray')
plt.axis('off')
for layer_line in layer_lines:
    plt.plot(layer_line, linewidth = 3)
layer_segmentation layer_segmentation

Source code in qim3d/processing/layers2d.py 
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
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):
    """
    Works on 2D and 3D data.
    Light one function wrapper around slgbuilder https://github.com/Skielex/slgbuilder to do layer segmentation
    Now uses only MaxflowBuilder for solving.

    Args:
        data: 2D or 3D array on which it will be computed
        inverted: if True, it will invert the brightness of the image
        n_layers: How many layers are we looking for (result in a layer and background)
        delta: Smoothness parameter
        min_margin: If we want more layers, we have to have a margin otherwise they are all going to be exactly the same
        max_margin: Maximum margin between layers
        wrap: If True, starting and ending point of the border between layers are at the same level

    Returns:
        segmentations: list of numpy arrays, even if n_layers == 1, each array is only 0s and 1s, 1s segmenting this specific layer

    Raises:
        TypeError: If Data is not np.array, if n_layers is not integer.
        ValueError: If n_layers is less than 1, if delta is negative or zero

    Example:
        Example is only shown on 2D image, but segment_layers can also take 3D structures.
        ```python
        import qim3d

        layers_image = qim3d.io.load('layers3d.tif')[:,:,0]
        layers = qim3d.processing.segment_layers(layers_image, n_layers = 2)
        layer_lines = qim3d.processing.get_lines(layers)

        from matplotlib import pyplot as plt

        plt.imshow(layers_image, cmap='gray')
        plt.axis('off')
        for layer_line in layer_lines:
            plt.plot(layer_line, linewidth = 3)
        ```
        ![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

qim3d.processing.create_mesh 

create_mesh(volume, level=None, step_size=1, allow_degenerate=False, padding=(2, 2, 2), **kwargs)

Convert a volume to a mesh using the Marching Cubes algorithm, with optional thresholding and padding.

Parameters:

Name Type Description Default
volume ndarray

The 3D numpy array representing the volume.

 required
level float

The threshold value for Marching Cubes. If None, Otsu's method is used.

 None
step_size int

The step size for the Marching Cubes algorithm.

 1
allow_degenerate bool

Whether to allow degenerate (i.e. zero-area) triangles in the end-result.

 False
padding tuple of int

Padding to add around the volume.

 (2, 2, 2)
**kwargs Any

Additional keyword arguments to pass to skimage.measure.marching_cubes.

 {}

Returns:

Name Type Description
trimesh Trimesh

The generated mesh.
Example 
import qim3d
vol = qim3d.generate.blob(base_shape=(128,128,128),
                          final_shape=(128,128,128),
                          noise_scale=0.03,
                          order=1,
                          gamma=1,
                          max_value=255,
                          threshold=0.5,
                          dtype='uint8'
                          )
mesh = qim3d.processing.create_mesh(vol, step_size=3)
qim3d.viz.mesh(mesh.vertices, mesh.faces)
Source code in qim3d/processing/mesh.py 
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
def create_mesh(
    volume: np.ndarray,
    level: float = None,
    step_size=1,
    allow_degenerate=False,
    padding: Tuple[int, int, int] = (2, 2, 2),
    **kwargs: Any,
) -> trimesh.Trimesh:
    """
    Convert a volume to a mesh using the Marching Cubes algorithm, with optional thresholding and padding.

    Args:
        volume (np.ndarray): The 3D numpy array representing the volume.
        level (float, optional): The threshold value for Marching Cubes. If None, Otsu's method is used.
        step_size (int, optional): The step size for the Marching Cubes algorithm.
        allow_degenerate (bool, optional): Whether to allow degenerate (i.e. zero-area) triangles in the end-result.
        If False, degenerate triangles are removed, at the cost of making the algorithm slower. Default False.
        padding (tuple of int, optional): Padding to add around the volume.
        **kwargs: Additional keyword arguments to pass to `skimage.measure.marching_cubes`.

    Returns:
        trimesh: The generated mesh.

    Example:
        ```python
        import qim3d
        vol = qim3d.generate.blob(base_shape=(128,128,128),
                                  final_shape=(128,128,128),
                                  noise_scale=0.03,
                                  order=1,
                                  gamma=1,
                                  max_value=255,
                                  threshold=0.5,
                                  dtype='uint8'
                                  )
        mesh = qim3d.processing.create_mesh(vol, step_size=3)
        qim3d.viz.mesh(mesh.vertices, mesh.faces)
        ```
        <iframe src="https://platform.qim.dk/k3d/mesh_visualization.html" width="100%" height="500" frameborder="0"></iframe>
    """
    if volume.ndim != 3:
        raise ValueError("The input volume must be a 3D numpy array.")

    # Compute the threshold level if not provided
    if level is None:
        level = filters.threshold_otsu(volume)
        log.info(f"Computed level using Otsu's method: {level}")

    # Apply padding to the volume
    if padding is not None:
        pad_z, pad_y, pad_x = padding
        padding_value = np.min(volume)
        volume = np.pad(
            volume,
            ((pad_z, pad_z), (pad_y, pad_y), (pad_x, pad_x)),
            mode="constant",
            constant_values=padding_value,
        )
        log.info(f"Padded volume with {padding} to shape: {volume.shape}")

    # Call skimage.measure.marching_cubes with user-provided kwargs
    verts, faces, normals, values = measure.marching_cubes(
        volume, level=level, step_size=step_size, allow_degenerate=allow_degenerate, **kwargs
    )

    # Create the Trimesh object
    mesh = trimesh.Trimesh(vertices=verts, faces=faces)

    # Fix face orientation to ensure normals point outwards
    trimesh.repair.fix_inversion(mesh, multibody=True) 

    return mesh

qim3d.processing.Pipeline

Example 

import qim3d
from qim3d.processing import Pipeline, Median, Gaussian, Maximum, Minimum

# Get data
vol = qim3d.examples.fly_150x256x256

# Show original
qim3d.viz.slices(vol, axis=0, show=True)

# Create filter pipeline
pipeline = Pipeline(
    Median(size=5),
    Gaussian(sigma=3, dask = True)
)

# Append a third filter to the pipeline
pipeline.append(Maximum(size=3))

# Apply filter pipeline
vol_filtered = pipeline(vol)

# Show filtered
qim3d.viz.slices(vol_filtered, axis=0)
original volume filtered volume

Source code in qim3d/processing/filters.py 
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
class Pipeline:
    """
    Example:
        ```python
        import qim3d
        from qim3d.processing import Pipeline, Median, Gaussian, Maximum, Minimum

        # Get data
        vol = qim3d.examples.fly_150x256x256

        # Show original
        qim3d.viz.slices(vol, axis=0, show=True)

        # Create filter pipeline
        pipeline = Pipeline(
            Median(size=5),
            Gaussian(sigma=3, dask = True)
        )

        # Append a third filter to the pipeline
        pipeline.append(Maximum(size=3))

        # Apply filter pipeline
        vol_filtered = pipeline(vol)

        # Show filtered
        qim3d.viz.slices(vol_filtered, axis=0)
        ```
        ![original volume](assets/screenshots/filter_original.png)
        ![filtered volume](assets/screenshots/filter_processed.png)

        """
    def __init__(self, *args: Type[FilterBase]):
        """
        Represents a sequence of image filters.

        Args:
            *args: Variable number of filter instances to be applied sequentially.
        """
        self.filters = {}

        for idx, fn in enumerate(args):
            self._add_filter(str(idx), fn)

    def _add_filter(self, name: str, fn: Type[FilterBase]):
        """
        Adds a filter to the sequence.

        Args:
            name: A string representing the name or identifier of the filter.
            fn: An instance of a FilterBase subclass.

        Raises:
            AssertionError: If `fn` is not an instance of the FilterBase class.
        """
        if not isinstance(fn, FilterBase):
            filter_names = [
                subclass.__name__ for subclass in FilterBase.__subclasses__()
            ]
            raise AssertionError(
                f"filters should be instances of one of the following classes: {filter_names}"
            )
        self.filters[name] = fn

    def append(self, fn: Type[FilterBase]):
        """
        Appends a filter to the end of the sequence.

        Args:
            fn: An instance of a FilterBase subclass to be appended.

        Example:
            ```python
            import qim3d
            from qim3d.processing import Pipeline, Maximum, Median

            # Create filter pipeline
            pipeline = Pipeline(
                Maximum(size=3, dask=True),
            )

            # Append a second filter to the pipeline
            pipeline.append(Median(size=5))
            ```
        """
        self._add_filter(str(len(self.filters)), fn)

    def __call__(self, input):
        """
        Applies the sequential filters to the input in order.

        Args:
            input: The input image or volume.

        Returns:
            The filtered image or volume after applying all sequential filters.
        """
        for fn in self.filters.values():
            input = fn(input)
        return input

qim3d.processing.Pipeline.append 

append(fn)

Appends a filter to the end of the sequence.

Parameters:

Name Type Description Default
fn Type[FilterBase]

An instance of a FilterBase subclass to be appended.

 required
Example 
import qim3d
from qim3d.processing import Pipeline, Maximum, Median

# Create filter pipeline
pipeline = Pipeline(
    Maximum(size=3, dask=True),
)

# Append a second filter to the pipeline
pipeline.append(Median(size=5))
Source code in qim3d/processing/filters.py 
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
def append(self, fn: Type[FilterBase]):
    """
    Appends a filter to the end of the sequence.

    Args:
        fn: An instance of a FilterBase subclass to be appended.

    Example:
        ```python
        import qim3d
        from qim3d.processing import Pipeline, Maximum, Median

        # Create filter pipeline
        pipeline = Pipeline(
            Maximum(size=3, dask=True),
        )

        # Append a second filter to the pipeline
        pipeline.append(Median(size=5))
        ```
    """
    self._add_filter(str(len(self.filters)), fn)

qim3d.processing.operations

qim3d.processing.operations.remove_background 

remove_background(vol, median_filter_size=2, min_object_radius=3, background='dark', **median_kwargs)

Remove background from a volume using a qim3d filters.

Parameters:

Name Type Description Default
vol ndarray

The volume to remove background from.

 required
median_filter_size int

The size of the median filter. Defaults to 2.

 2
min_object_radius int

The radius of the structuring element for the tophat filter. Defaults to 3.

 3
background str

The background type. Can be 'dark' or 'bright'. Defaults to 'dark'.

 'dark'
**median_kwargs

Additional keyword arguments for the Median filter.

 {}

Returns:

Type Description
ndarray

np.ndarray: The volume with background removed.
Example 

import qim3d

vol = qim3d.examples.cement_128x128x128
qim3d.viz.slices(vol, vmin=0, vmax=255)
operations-remove_background_before

vol_filtered  = qim3d.processing.operations.remove_background(vol,
                                                      min_object_radius=3,
                                                      background="bright")
qim3d.viz.slices(vol_filtered, vmin=0, vmax=255)
operations-remove_background_after

Source code in qim3d/processing/operations.py 
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
def remove_background(
    vol: np.ndarray,
    median_filter_size: int = 2,
    min_object_radius: int = 3,
    background: str = "dark",
    **median_kwargs,
) -> np.ndarray:
    """
    Remove background from a volume using a qim3d filters.

    Args:
        vol (np.ndarray): The volume to remove background from.
        median_filter_size (int, optional): The size of the median filter. Defaults to 2.
        min_object_radius (int, optional): The radius of the structuring element for the tophat filter. Defaults to 3.
        background (str, optional): The background type. Can be 'dark' or 'bright'. Defaults to 'dark'.
        **median_kwargs: Additional keyword arguments for the Median filter.

    Returns:
        np.ndarray: The volume with background removed.


    Example:
        ```python
        import qim3d

        vol = qim3d.examples.cement_128x128x128
        qim3d.viz.slices(vol, vmin=0, vmax=255)
        ```
        ![operations-remove_background_before](assets/screenshots/operations-remove_background_before.png)

        ```python
        vol_filtered  = qim3d.processing.operations.remove_background(vol,
                                                              min_object_radius=3,
                                                              background="bright")
        qim3d.viz.slices(vol_filtered, vmin=0, vmax=255)
        ```
        ![operations-remove_background_after](assets/screenshots/operations-remove_background_after.png)
    """

    # Create a pipeline with a median filter and a tophat filter
    pipeline = filters.Pipeline(
        filters.Median(size=median_filter_size, **median_kwargs),
        filters.Tophat(radius=min_object_radius, background=background),
    )

    # Apply the pipeline to the volume
    return pipeline(vol)

qim3d.processing.operations.watershed 

watershed(bin_vol, min_distance=5)

Apply watershed segmentation to a binary volume.

Parameters:

Name Type Description Default
bin_vol ndarray

Binary volume to segment. The input should be a 3D binary image where non-zero elements represent the objects to be segmented.

 required
min_distance int

Minimum number of pixels separating peaks in the distance transform. Peaks that are too close will be merged, affecting the number of segmented objects. Default is 5.

 5

Returns:

Type Description
tuple[ndarray, int]

tuple[np.ndarray, int]: - Labeled volume (np.ndarray): A 3D array of the same shape as the input bin_vol, where each segmented object is assigned a unique integer label. - num_labels (int): The total number of unique objects found in the labeled volume.
Example 

import qim3d

vol = qim3d.examples.cement_128x128x128
binary = qim3d.processing.filters.gaussian(vol, sigma = 2)<60

qim3d.viz.slices(binary, axis=1)
operations-watershed_before

labeled_volume, num_labels = qim3d.processing.operations.watershed(binary)

cmap = qim3d.viz.colormaps.objects(num_labels)
qim3d.viz.slices(labeled_volume, axis = 1, cmap = cmap)
operations-watershed_after

Source code in qim3d/processing/operations.py 
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
def watershed(bin_vol: np.ndarray, min_distance: int = 5) -> tuple[np.ndarray, int]:
    """
    Apply watershed segmentation to a binary volume.

    Args:
        bin_vol (np.ndarray): Binary volume to segment. The input should be a 3D binary image where non-zero elements 
                              represent the objects to be segmented.
        min_distance (int): Minimum number of pixels separating peaks in the distance transform. Peaks that are 
                            too close will be merged, affecting the number of segmented objects. Default is 5.

    Returns:
        tuple[np.ndarray, int]: 
            - Labeled volume (np.ndarray): A 3D array of the same shape as the input `bin_vol`, where each segmented object
              is assigned a unique integer label.
            - num_labels (int): The total number of unique objects found in the labeled volume.

    Example:
        ```python
        import qim3d

        vol = qim3d.examples.cement_128x128x128
        binary = qim3d.processing.filters.gaussian(vol, sigma = 2)<60

        qim3d.viz.slices(binary, axis=1)
        ```
        ![operations-watershed_before](assets/screenshots/operations-watershed_before.png)

        ```python
        labeled_volume, num_labels = qim3d.processing.operations.watershed(binary)

        cmap = qim3d.viz.colormaps.objects(num_labels)
        qim3d.viz.slices(labeled_volume, axis = 1, cmap = cmap)
        ```
        ![operations-watershed_after](assets/screenshots/operations-watershed_after.png)

    """
    import skimage
    import scipy

    # Compute distance transform of binary volume
    distance = scipy.ndimage.distance_transform_edt(bin_vol)

    # Find peak coordinates in distance transform
    coords = skimage.feature.peak_local_max(
        distance, min_distance=min_distance, labels=bin_vol
    )

    # Create a mask with peak coordinates
    mask = np.zeros(distance.shape, dtype=bool)
    mask[tuple(coords.T)] = True

    # Label peaks
    markers, _ = scipy.ndimage.label(mask)

    # Apply watershed segmentation
    labeled_volume = skimage.segmentation.watershed(
        -distance, markers=markers, mask=bin_vol
    )

    # Extract number of objects found
    num_labels = len(np.unique(labeled_volume)) - 1
    log.info(f"Total number of objects found: {num_labels}")

    return labeled_volume, num_labels

qim3d.processing.operations.fade_mask 

fade_mask(vol, decay_rate=10, ratio=0.5, geometry='spherical', invert=False, axis=0)

Apply edge fading to a volume.

Parameters:

Name Type Description Default
vol ndarray

The volume to apply edge fading to.

 required
decay_rate float

The decay rate of the fading. Defaults to 10.

 10
ratio float

The ratio of the volume to fade. Defaults to 0.5.

 0.5
geometric str

The geometric shape of the fading. Can be 'spherical' or 'cylindrical'. Defaults to 'spherical'.

 required
axis int

The axis along which to apply the fading. Defaults to 0.

 0

Returns:

Name Type Description
vol_faded ndarray

The volume with edge fading applied.
Example 

import qim3d
qim3d.viz.vol(vol)
Image before edge fading has visible artifacts from the support. Which obscures the object of interest. operations-edge_fade_before

import qim3d
vol_faded = qim3d.processing.operations.edge_fade(vol, decay_rate=4, ratio=0.45, geometric='cylindrical')
qim3d.viz.vol(vol_faded)
Afterwards the artifacts are faded out, making the object of interest more visible for visualization purposes. operations-edge_fade_after

Source code in qim3d/processing/operations.py 
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
def fade_mask(
    vol: np.ndarray,
    decay_rate: float = 10,
    ratio: float = 0.5,
    geometry: str = "spherical",
    invert: bool = False,
    axis: int = 0,
) -> np.ndarray:
    """
    Apply edge fading to a volume.

    Args:
        vol (np.ndarray): The volume to apply edge fading to.
        decay_rate (float, optional): The decay rate of the fading. Defaults to 10.
        ratio (float, optional): The ratio of the volume to fade. Defaults to 0.5.
        geometric (str, optional): The geometric shape of the fading. Can be 'spherical' or 'cylindrical'. Defaults to 'spherical'.
        axis (int, optional): The axis along which to apply the fading. Defaults to 0.

    Returns:
        vol_faded (np.ndarray): The volume with edge fading applied.

    Example:
        ```python
        import qim3d
        qim3d.viz.vol(vol)
        ```
        Image before edge fading has visible artifacts from the support. Which obscures the object of interest.
        ![operations-edge_fade_before](assets/screenshots/operations-edge_fade_before.png)

        ```python
        import qim3d
        vol_faded = qim3d.processing.operations.edge_fade(vol, decay_rate=4, ratio=0.45, geometric='cylindrical')
        qim3d.viz.vol(vol_faded)
        ```
        Afterwards the artifacts are faded out, making the object of interest more visible for visualization purposes.
        ![operations-edge_fade_after](assets/screenshots/operations-edge_fade_after.png)

    """
    if 0 > axis or axis >= vol.ndim:
        raise ValueError(
            "Axis must be between 0 and the number of dimensions of the volume"
        )

    # Generate the coordinates of each point in the array
    shape = vol.shape
    z, y, x = np.indices(shape)

    # Calculate the center of the array
    center = np.array([(s - 1) / 2 for s in shape])

    # Calculate the distance of each point from the center
    if geometry == "spherical":
        distance = np.linalg.norm([z - center[0], y - center[1], x - center[2]], axis=0)
    elif geometry == "cylindrical":
        distance_list = np.array([z - center[0], y - center[1], x - center[2]])
        # remove the axis along which the fading is not applied
        distance_list = np.delete(distance_list, axis, axis=0)
        distance = np.linalg.norm(distance_list, axis=0)
    else:
        raise ValueError("geometric must be 'spherical' or 'cylindrical'")

    # Normalize the distances so that they go from 0 at the center to 1 at the farthest point
    max_distance = np.linalg.norm(center)
    normalized_distance = distance / (max_distance * ratio)

    # Apply the decay rate
    faded_distance = normalized_distance**decay_rate

    # Invert the distances to have 1 at the center and 0 at the edges
    fade_array = 1 - faded_distance
    fade_array[fade_array <= 0] = 0

    if invert:
        fade_array = -(fade_array - 1)

    # Apply the fading to the volume
    vol_faded = vol * fade_array

    return vol_faded

qim3d.processing.operations.overlay_rgb_images 

overlay_rgb_images(background, foreground, alpha=0.5, hide_black=True)

Overlay an RGB foreground onto an RGB background using alpha blending.

Parameters:

Name Type Description Default
background ndarray

The background RGB image.

 required
foreground ndarray

The foreground RGB image (usually masks).

 required
alpha float

The alpha value for blending. Defaults to 0.5.

 0.5
hide_black bool

If True, black pixels will have alpha value 0, so the black won't be visible. Used for segmentation where we don't care about background. Defaults to True.

 True

Returns:

Name Type Description
composite ndarray

The composite RGB image with overlaid foreground.

Raises:

Type Description
ValueError

If input images have different shapes.
Note
  • The function performs alpha blending to overlay the foreground onto the background.
  • It ensures that the background and foreground have the same first two dimensions (image size matches).
  • It can handle greyscale images, values from 0 to 1, raw values which are negative or bigger than 255.
  • It calculates the maximum projection of the foreground and blends them onto the background.
Source code in qim3d/processing/operations.py 
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
def overlay_rgb_images(
    background: np.ndarray, foreground: np.ndarray, alpha: float = 0.5, hide_black:bool = True,
) -> np.ndarray:
    """
    Overlay an RGB foreground onto an RGB background using alpha blending.

    Args:
        background (numpy.ndarray): The background RGB image.
        foreground (numpy.ndarray): The foreground RGB image (usually masks).
        alpha (float, optional): The alpha value for blending. Defaults to 0.5.
        hide_black (bool, optional): If True, black pixels will have alpha value 0, so the black won't be visible. Used for segmentation where we don't care about background. Defaults to True.

    Returns:
        composite (numpy.ndarray): The composite RGB image with overlaid foreground.

    Raises:
        ValueError: If input images have different shapes.

    Note:
        - The function performs alpha blending to overlay the foreground onto the background.
        - It ensures that the background and foreground have the same first two dimensions (image size matches).
        - It can handle greyscale images, values from 0 to 1, raw values which are negative or bigger than 255.
        - It calculates the maximum projection of the foreground and blends them onto the background.
    """

    def to_uint8(image:np.ndarray):
        if np.min(image) < 0:
            image = image - np.min(image)

        maxim = np.max(image)
        if maxim > 255:
            image = (image / maxim)*255
        elif maxim <= 1:
            image = image*255

        if image.ndim == 2:
            image = np.repeat(image[..., None], 3, -1)
        elif image.ndim == 3:
            image = image[..., :3] # Ignoring alpha channel
        else:
            raise ValueError(F'Input image can not have higher dimension than 3. Yours have {image.ndim}')

        return image.astype(np.uint8)

    background = to_uint8(background)
    foreground = to_uint8(foreground)

    # Ensure both images have the same shape
    if background.shape != foreground.shape:
        raise ValueError(F"Input images must have the same first two dimensions. But background is of shape {background.shape} and foreground is of shape {foreground.shape}")

    # Perform alpha blending
    foreground_max_projection = np.amax(foreground, axis=2)
    foreground_max_projection = np.stack((foreground_max_projection,) * 3, axis=-1)

    # Normalize if we have something
    if np.max(foreground_max_projection) > 0:
        foreground_max_projection = foreground_max_projection / np.max(
            foreground_max_projection
        )
    # Check alpha validity
    if alpha < 0:
        raise ValueError(F'Alpha has to be positive number. You used {alpha}')
    elif alpha > 1:
        alpha = 1

    # If the pixel is black, its alpha value is set to 0, so it has no effect on the image
    if hide_black:
        alpha = np.full((background.shape[0], background.shape[1],1), alpha)
        alpha[np.apply_along_axis(lambda x: (x == [0,0,0]).all(), axis = 2, arr = foreground)] = 0

    composite = background * (1 - alpha) + foreground * alpha
    composite = np.clip(composite, 0, 255).astype("uint8")

    return composite.astype("uint8")

qim3d.processing.features

qim3d.processing.features.area 

area(obj, **mesh_kwargs)

Compute the surface area of a 3D volume or mesh.

Parameters:

Name Type Description Default
obj

Either a np.ndarray volume or a mesh object of type trimesh.Trimesh.

 required
**mesh_kwargs

Additional arguments for mesh creation if the input is a volume.

 {}

Returns:

Name Type Description
area float

The surface area of the object.
Example

Compute area from a mesh: 

import qim3d

# Load a mesh from a file
mesh = qim3d.io.load_mesh('path/to/mesh.obj')

# Compute the surface area of the mesh
area = qim3d.processing.area(mesh)
print(f"Area: {area}")

Compute area from a np.ndarray: 

import qim3d

# Generate a 3D blob
synthetic_blob = qim3d.generate.blob(noise_scale = 0.015)

# Compute the surface area of the blob
volume = qim3d.processing.area(synthetic_blob, level=0.5)
print('Area:', volume)

Source code in qim3d/processing/features.py 
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
def area(obj, **mesh_kwargs) -> float:
    """
    Compute the surface area of a 3D volume or mesh.

    Args:
        obj: Either a np.ndarray volume or a mesh object of type trimesh.Trimesh.
        **mesh_kwargs: Additional arguments for mesh creation if the input is a volume.

    Returns:
        area: The surface area of the object.

    Example:
        Compute area from a mesh:
        ```python
        import qim3d

        # Load a mesh from a file
        mesh = qim3d.io.load_mesh('path/to/mesh.obj')

        # Compute the surface area of the mesh
        area = qim3d.processing.area(mesh)
        print(f"Area: {area}")
        ```

        Compute area from a np.ndarray:
        ```python
        import qim3d

        # Generate a 3D blob
        synthetic_blob = qim3d.generate.blob(noise_scale = 0.015)

        # Compute the surface area of the blob
        volume = qim3d.processing.area(synthetic_blob, level=0.5)
        print('Area:', volume)
        ```
    """
    if isinstance(obj, np.ndarray):
        log.info("Converting volume to mesh.")
        obj = qim3d.processing.create_mesh(obj, **mesh_kwargs)

    return obj.area

qim3d.processing.features.volume 

volume(obj, **mesh_kwargs)

Compute the volume of a 3D volume or mesh.

Parameters:

Name Type Description Default
obj

Either a np.ndarray volume or a mesh object of type trimesh.Trimesh.

 required
**mesh_kwargs

Additional arguments for mesh creation if the input is a volume.

 {}

Returns:

Name Type Description
volume float

The volume of the object.
Example

Compute volume from a mesh: 

import qim3d

# Load a mesh from a file
mesh = qim3d.io.load_mesh('path/to/mesh.obj')

# Compute the volume of the mesh
volume = qim3d.processing.volume(mesh)
print('Volume:', volume)

Compute volume from a np.ndarray: 

import qim3d

# Generate a 3D blob
synthetic_blob = qim3d.generate.blob(noise_scale = 0.015)

# Compute the volume of the blob
volume = qim3d.processing.volume(synthetic_blob, level=0.5)
print('Volume:', volume)

Source code in qim3d/processing/features.py 
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
def volume(obj, **mesh_kwargs) -> float:
    """
    Compute the volume of a 3D volume or mesh.

    Args:
        obj: Either a np.ndarray volume or a mesh object of type trimesh.Trimesh.
        **mesh_kwargs: Additional arguments for mesh creation if the input is a volume.

    Returns:
        volume: The volume of the object.

    Example:
        Compute volume from a mesh:
        ```python
        import qim3d

        # Load a mesh from a file
        mesh = qim3d.io.load_mesh('path/to/mesh.obj')

        # Compute the volume of the mesh
        volume = qim3d.processing.volume(mesh)
        print('Volume:', volume)
        ```

        Compute volume from a np.ndarray:
        ```python
        import qim3d

        # Generate a 3D blob
        synthetic_blob = qim3d.generate.blob(noise_scale = 0.015)

        # Compute the volume of the blob
        volume = qim3d.processing.volume(synthetic_blob, level=0.5)
        print('Volume:', volume)
        ```

    """
    if isinstance(obj, np.ndarray):
        log.info("Converting volume to mesh.")
        obj = qim3d.processing.create_mesh(obj, **mesh_kwargs)

    return obj.volume

qim3d.processing.features.sphericity 

sphericity(obj, **mesh_kwargs)

Compute the sphericity of a 3D volume or mesh.

Sphericity is a measure of how spherical an object is. It is defined as the ratio of the surface area of a sphere with the same volume as the object to the object's actual surface area.

Parameters:

Name Type Description Default
obj

Either a np.ndarray volume or a mesh object of type trimesh.Trimesh.

 required
**mesh_kwargs

Additional arguments for mesh creation if the input is a volume.

 {}

Returns:

Name Type Description
sphericity float

A float value representing the sphericity of the object.
Example

Compute sphericity from a mesh: 

import qim3d

# Load a mesh from a file
mesh = qim3d.io.load_mesh('path/to/mesh.obj')

# Compute the sphericity of the mesh
sphericity = qim3d.processing.sphericity(mesh)

Compute sphericity from a np.ndarray: 

import qim3d

# Generate a 3D blob
synthetic_blob = qim3d.generate.blob(noise_scale = 0.015)

# Compute the sphericity of the blob
sphericity = qim3d.processing.sphericity(synthetic_blob, level=0.5)

Limitations due to pixelation

Sphericity is particularly sensitive to the resolution of the mesh, as it directly impacts the accuracy of surface area and volume calculations. Since the mesh is generated from voxel-based 3D volume data, the discrete nature of the voxels leads to pixelation effects that reduce the precision of sphericity measurements. Higher resolution meshes may mitigate these errors but often at the cost of increased computational demands.

Source code in qim3d/processing/features.py 
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
def sphericity(obj, **mesh_kwargs) -> float:
    """
    Compute the sphericity of a 3D volume or mesh.

    Sphericity is a measure of how spherical an object is. It is defined as the ratio
    of the surface area of a sphere with the same volume as the object to the object's
    actual surface area.

    Args:
        obj: Either a np.ndarray volume or a mesh object of type trimesh.Trimesh.
        **mesh_kwargs: Additional arguments for mesh creation if the input is a volume.

    Returns:
        sphericity: A float value representing the sphericity of the object.

    Example:
        Compute sphericity from a mesh:
        ```python
        import qim3d

        # Load a mesh from a file
        mesh = qim3d.io.load_mesh('path/to/mesh.obj')

        # Compute the sphericity of the mesh
        sphericity = qim3d.processing.sphericity(mesh)
        ```

        Compute sphericity from a np.ndarray:
        ```python
        import qim3d

        # Generate a 3D blob
        synthetic_blob = qim3d.generate.blob(noise_scale = 0.015)

        # Compute the sphericity of the blob
        sphericity = qim3d.processing.sphericity(synthetic_blob, level=0.5)
        ```

    !!! info "Limitations due to pixelation"
        Sphericity is particularly sensitive to the resolution of the mesh, as it directly impacts the accuracy of surface area and volume calculations.
        Since the mesh is generated from voxel-based 3D volume data, the discrete nature of the voxels leads to pixelation effects that reduce the precision of sphericity measurements.
        Higher resolution meshes may mitigate these errors but often at the cost of increased computational demands.
    """
    if isinstance(obj, np.ndarray):
        log.info("Converting volume to mesh.")
        obj = qim3d.processing.create_mesh(obj, **mesh_kwargs)

    volume = qim3d.processing.volume(obj)
    area = qim3d.processing.area(obj)

    if area == 0:
        log.warning("Surface area is zero, sphericity is undefined.")
        return np.nan

    sphericity = (np.pi ** (1 / 3) * (6 * volume) ** (2 / 3)) / area
    log.info(f"Sphericity: {sphericity}")
    return sphericity