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.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. Defaults to 1.0.

 1.0
rho float

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

 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 Any

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: bool = 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. Defaults to 1.0.
        rho (float, optional): An integration scale giving the size over the neighborhood in which the orientation is to be analysed. Defaults to 6.0.
        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 (Any): 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 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 or None

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 Any

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

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

import qim3d

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

# Extract one slice to show that localthickness works on 2D slices too
slice = vol[:,:,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
def local_thickness(
    image: np.ndarray,
    scale: float = 1,
    mask: Optional[np.ndarray] = None,
    visualize: bool = 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 or None, 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 (Any): 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

        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.noise_object_collection(num_objects=15)

        # Extract one slice to show that localthickness 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 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_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 of arrays

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/_layers.py 
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
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): List of arrays where each array is 2D segmentation with only 2 classes

    Returns:
        segmentation_lines (list): 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. Defaults to False

 False
n_layers int

Determines amount of layers to look for (result in a layer and background). Defaults to 1.

 1
delta float

Patameter determining smoothness. Defaults to 1.

 1
min_margin int or None

Parameter for minimum margin. If more layers are wanted, a margin is necessary to avoid layers being identical. Defaults to None.

 10
max_margin int or None

Parameter for maximum margin. If more layers are wanted, a margin is necessary to avoid layers being identical. Defaults to None.

 None
wrap bool

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

 False

Returns:

Name Type Description
segmentations list[ndarray]

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)

import matplotlib.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/_layers.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
84
85
86
87
88
89
90
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:
    """
    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 (np.ndarray): 2D or 3D array on which it will be computed
        inverted (bool): If True, it will invert the brightness of the image. Defaults to False
        n_layers (int): Determines amount of layers to look for (result in a layer and background). Defaults to 1.
        delta (float): Patameter determining smoothness. Defaults to 1.
        min_margin (int or None): Parameter for minimum margin. If more layers are wanted, a margin is necessary to avoid layers being identical. Defaults to None.
        max_margin (int or None): Parameter for maximum margin. If more layers are wanted, a margin is necessary to avoid layers being identical. Defaults to None.
        wrap (bool): If True, starting and ending point of the border between layers are at the same level. Defaults to False.

    Returns:
        segmentations (list[np.ndarray]): 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)

        import matplotlib.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