Plot and render

Visualization of volumetric data.

qim3d.viz.export_rotation

export_rotation(
    path,
    volume,
    degrees=360,
    n_frames=180,
    fps=30,
    image_size=(256, 256),
    colormap='magma',
    camera_height=2.0,
    camera_distance='auto',
    camera_focus='center',
    show=False,
)

Exports a 360-degree turntable animation of the volume to a video or GIF.

Generates a spinning orbit visualization of the 3D data, perfect for presentations, reports, or sharing results on the web. It renders the volume from a rotating camera perspective and saves the output as a movie file (.mp4, .webm, .avi) or an animated .gif.

Key Features:

• Presentation Ready: Creates smooth, professional animations of your data.
• Flexible Output: Supports common video formats and high-quality GIFs.
• Customizable Camera: Control the height, distance, and focus point of the rotation.

Parameters:

Name	Type	Description	Default
path	str	The destination file path. Must end with .gif, .avi, .mp4, or .webm. If no extension is provided, defaults to .gif.	required
volume	ndarray	The 3D input volume to be animated.	required
degrees	int	Total rotation angle in degrees (e.g., 360 for a full spin).	360
n_frames	int	Total number of frames to generate. Higher values create smoother/slower animations at fixed FPS.	180
fps	int	Frames per second. Controls the playback speed.	30
image_size	tuple[int, int] or None	Resolution (width, height) of the output frames.	(256, 256)
colormap	str	Matplotlib colormap name for the volume rendering.	'magma'
camera_height	float	Vertical position of the camera relative to the volume's Z-axis height.	2.0
camera_distance	float or str	Distance from the camera to the focus point. str: Use 'auto' to automatically calculate a fitting distance. float: Specific distance in voxel units.	'auto'
camera_focus	list or str	The point the camera rotates around. str: Use 'center' to rotate around the volume center. list: A list of 3 integers [z, y, x] specifying the voxel coordinate.	'center'
show	bool	If True, displays the generated animation inline in the notebook.	False

Raises:

Type	Description
TypeError	If camera_focus or camera_distance is invalid.
ValueError	If path contains an unsupported file extension.

Example

Creation of .gif file with default parameters of a generated volume.

import qim3d
vol = qim3d.generate.volume()

qim3d.viz.export_rotation('test.gif', vol, show=True)

Example

Creation of a .webm file with specified parameters of a generated volume in the shape of a tube.

import qim3d

vol = qim3d.generate.volume(shape='tube')

qim3d.viz.export_rotation('test.webm', vol,
                          degrees = 360,
                          n_frames = 120,
                          fps = 30,
                          image_size = (512,512),
                          camera_height = 3.0,
                          camera_distance = 'auto',
                          camera_focus = 'center',
                          show = True)

Source code in qim3d/viz/_data_exploration.py

def export_rotation(
    path: str,
    volume: np.ndarray,
    degrees: int = 360,
    n_frames: int = 180,
    fps: int = 30,
    image_size: tuple[int, int] | None = (256, 256),
    colormap: str = 'magma',
    camera_height: float = 2.0,
    camera_distance: float | str = 'auto',
    camera_focus: list | str = 'center',
    show: bool = False,
) -> None:
    """
    Exports a 360-degree turntable animation of the volume to a video or GIF.

    Generates a spinning orbit visualization of the 3D data, perfect for presentations, reports, or sharing results on the web. It renders the volume from a rotating camera perspective and saves the output as a movie file (.mp4, .webm, .avi) or an animated .gif.

    **Key Features:**

    * **Presentation Ready:** Creates smooth, professional animations of your data.
    * **Flexible Output:** Supports common video formats and high-quality GIFs.
    * **Customizable Camera:** Control the height, distance, and focus point of the rotation.

    Args:
        path (str): The destination file path. Must end with .gif, .avi, .mp4, or .webm. If no extension is provided, defaults to .gif.
        volume (numpy.ndarray): The 3D input volume to be animated.
        degrees (int, optional): Total rotation angle in degrees (e.g., 360 for a full spin).
        n_frames (int, optional): Total number of frames to generate. Higher values create smoother/slower animations at fixed FPS.
        fps (int, optional): Frames per second. Controls the playback speed.
        image_size (tuple[int, int] or None, optional): Resolution (width, height) of the output frames.
        colormap (str, optional): Matplotlib colormap name for the volume rendering.
        camera_height (float, optional): Vertical position of the camera relative to the volume's Z-axis height.
        camera_distance (float or str, optional): Distance from the camera to the focus point.

            * **str**: Use `'auto'` to automatically calculate a fitting distance.
            * **float**: Specific distance in voxel units.

        camera_focus (list or str, optional): The point the camera rotates around.

            * **str**: Use `'center'` to rotate around the volume center.
            * **list**: A list of 3 integers `[z, y, x]` specifying the voxel coordinate.

        show (bool, optional): If `True`, displays the generated animation inline in the notebook.

    Raises:
        TypeError: If `camera_focus` or `camera_distance` is invalid.
        ValueError: If `path` contains an unsupported file extension.

    Example:
        Creation of .gif file with default parameters of a generated volume.
        ```python
        import qim3d
        vol = qim3d.generate.volume()

        qim3d.viz.export_rotation('test.gif', vol, show=True)
        ```
        ![export_rotation_defaults](../../assets/screenshots/export_rotation_defaults.gif)

    Example:
        Creation of a .webm file with specified parameters of a generated volume in the shape of a tube.
        ```python
        import qim3d

        vol = qim3d.generate.volume(shape='tube')

        qim3d.viz.export_rotation('test.webm', vol,
                                  degrees = 360,
                                  n_frames = 120,
                                  fps = 30,
                                  image_size = (512,512),
                                  camera_height = 3.0,
                                  camera_distance = 'auto',
                                  camera_focus = 'center',
                                  show = True)
        ```
        ![export_rotation_video](../../assets/screenshots/export_rotation_video.gif)
    """
    if not (
        camera_focus == 'center'
        or (
            isinstance(camera_focus, list | np.ndarray)
            and not isinstance(camera_focus, str)
            and len(camera_focus) == 3
        )
    ):
        msg = f'Value "{camera_focus}" for camera focus is invalid. Use "center" or a list of three values.'
        raise TypeError(msg)
    if not (isinstance(camera_distance, float) or camera_distance == 'auto'):
        msg = f'Value "{camera_distance}" for camera distance is invalid. Use "auto" or a float value.'
        raise TypeError(msg)

    if Path(path).suffix == '':
        print(f'Input path: "{path}" does not have a filetype. Defaulting to .gif.')
        path += '.gif'

    # Handle img in (xyz) instead of (zyx) (due to rendering issues with the up-vector, ensure that z=y, such that we now have (x,z,y))
    vol = np.transpose(volume, (2, 0, 1))

    # Create a uniform grid
    grid = pv.ImageData()
    grid.dimensions = np.array(vol.shape) + 1  # PyVista dims are +1 from volume shape
    grid.spacing = (1, 1, 1)
    grid.origin = (0, 0, 0)
    grid.cell_data['values'] = vol.flatten(order='F')  # Fortran order

    # Initialize plotter
    plotter = pv.Plotter(off_screen=True)
    plotter.add_volume(grid, opacity='linear', cmap=colormap)
    plotter.remove_scalar_bar()  # Remove colorbar

    frames = []
    camera_height = vol.shape[1] * camera_height

    if camera_distance == 'auto':
        bounds = np.array(plotter.bounds)  # (xmin, xmax, ymin, ymax, zmin, zmax)
        diag = np.linalg.norm(
            [bounds[1] - bounds[0], bounds[3] - bounds[2], bounds[5] - bounds[4]]
        )
        camera_distance = diag * 2.0

    if camera_focus == 'center':
        _, center, _ = plotter.camera_position
    else:
        center = camera_focus

    center = np.array(center)

    angle_per_frame = degrees / n_frames
    radians_per_frame = np.radians(angle_per_frame)

    # Set up orbit radius and fixed up
    radius = camera_distance
    fixed_up = [0, 1, 0]
    for i in tqdm(range(n_frames), desc='Rendering'):
        theta = radians_per_frame * i
        x = radius * np.sin(theta)
        z = radius * np.cos(theta)
        y = camera_height  # fixed height

        eye = center + np.array([x, y, z])
        plotter.camera_position = [eye.tolist(), center.tolist(), fixed_up]

        plotter.render()
        img = plotter.screenshot(return_img=True, window_size=image_size)
        frames.append(img)

    if path[-4:] == '.gif':
        imageio.mimsave(path, frames, fps=fps, loop=0)

    elif path[-4:] == '.avi' or path[-4:] == '.mp4':
        writer = imageio.get_writer(path, fps=fps)
        for frame in frames:
            writer.append_data(frame)
        writer.close()

    elif path[-5:] == '.webm':
        writer = imageio.get_writer(
            path, fps=fps, codec='vp9', ffmpeg_params=['-crf', '32']
        )
        for frame in frames:
            writer.append_data(frame)
        writer.close()

    else:
        msg = 'Invalid file extension. Please use .gif, .avi, .mp4 or .webm'
        raise ValueError(msg)

    path = _get_save_path(path)
    log.info('File saved to ' + str(path.resolve()))

    if show:
        if path.suffix == '.gif':
            display(Image(filename=path))
        elif path.suffix in ['.avi', '.mp4', '.webm']:
            display(Video(filename=path, html_attributes='controls autoplay loop'))

qim3d.viz.circles

circles(
    blobs, volume, alpha=0.5, color='#ff9900', **kwargs
)

Visualizes detected blobs as circles overlaid on the volume slices.

This function is primarily used to verify the results of blob detection algorithms. It takes a list of detected features (defined by their center coordinates and radius) and projects them onto the 2D slices of the volume. As you scroll through the stack, the circles dynamically resize to represent the cross-section of the 3D spherical blobs at that specific depth, providing an intuitive check for detection accuracy.

Parameters:

Name	Type	Description	Default
blobs	ndarray	A list or array of detected blobs. Each blob is expected to be a 4-tuple or array (z, y, x, radius). This is typically the output from qim3d.detection.blobs.	required
volume	ndarray	The 3D volume (image stack) on which the blobs were detected.	required
alpha	float	The transparency level of the filled circles (0.0 to 1.0). Defaults to 0.5.	0.5
color	str	The color of the circles, capable of accepting hex codes or standard color names. Defaults to "#ff9900" (orange).	'#ff9900'
**kwargs	Any	Additional keyword arguments passed to the underlying qim3d.viz.slices_grid function (e.g., vmin, vmax).	{}

Returns:

Name	Type	Description
slicer_obj	interactive	An interactive widget with a slider to navigate through slices, showing the overlay of detected blobs.

Example

import qim3d
import qim3d.detection

# Get data
vol = qim3d.examples.cement_128x128x128

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

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

Source code in qim3d/viz/_detection.py

def circles(
    blobs: tuple[float, float, float, float],
    volume: np.ndarray,
    alpha: float = 0.5,
    color: str = '#ff9900',
    **kwargs,
) -> widgets.interactive:
    """
    Visualizes detected blobs as circles overlaid on the volume slices.

    This function is primarily used to verify the results of blob detection algorithms. It takes a list of detected features (defined by their center coordinates and radius) and projects them onto the 2D slices of the volume. As you scroll through the stack, the circles dynamically resize to represent the cross-section of the 3D spherical blobs at that specific depth, providing an intuitive check for detection accuracy.

    Args:
        blobs (np.ndarray): A list or array of detected blobs. Each blob is expected to be a 4-tuple or array `(z, y, x, radius)`. This is typically the output from `qim3d.detection.blobs`.
        volume (np.ndarray): The 3D volume (image stack) on which the blobs were detected.
        alpha (float, optional): The transparency level of the filled circles (0.0 to 1.0). Defaults to 0.5.
        color (str, optional): The color of the circles, capable of accepting hex codes or standard color names. Defaults to "#ff9900" (orange).
        **kwargs (Any): Additional keyword arguments passed to the underlying `qim3d.viz.slices_grid` function (e.g., `vmin`, `vmax`).

    Returns:
        slicer_obj (widgets.interactive):
            An interactive widget with a slider to navigate through slices, showing the overlay of detected blobs.

    Example:
        ```python
        import qim3d
        import qim3d.detection

        # Get data
        vol = qim3d.examples.cement_128x128x128

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

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

    def _slicer(z_slice):
        clear_output(wait=True)
        fig = qim3d.viz.slices_grid(
            volume[z_slice : z_slice + 1],
            n_slices=1,
            colormap='gray',
            display_figure=False,
            display_positions=False,
            **kwargs,
        )
        # Add circles from deteced blobs
        for detected in blobs:
            z, y, x, s = detected
            if abs(z - z_slice) < s:  # The blob is in the slice
                # Adjust the radius based on the distance from the center of the sphere
                distance_from_center = abs(z - z_slice)
                angle = (
                    np.pi / 2 * (distance_from_center / s)
                )  # Angle varies from 0 at the center to pi/2 at the edge
                adjusted_radius = s * np.cos(angle)  # Radius follows a cosine curve

                if adjusted_radius > 0.5:
                    c = plt.Circle(
                        (x, y),
                        adjusted_radius,
                        color=color,
                        linewidth=0,
                        fill=True,
                        alpha=alpha,
                    )
                    fig.get_axes()[0].add_patch(c)

        display(fig)
        return fig

    position_slider = widgets.IntSlider(
        value=volume.shape[0] // 2,
        min=0,
        max=volume.shape[0] - 1,
        description='Slice',
        continuous_update=True,
    )
    slicer_obj = widgets.interactive(_slicer, z_slice=position_slider)
    slicer_obj.layout = widgets.Layout(align_items='flex-start')

    return slicer_obj

qim3d.viz.local_thickness

local_thickness(
    image,
    image_lt,
    max_projection=False,
    axis=0,
    slice_index=None,
    show=False,
    figsize=(15, 5),
)

Visualizes a local thickness map alongside the original image and a statistics histogram.

This function provides a comprehensive view of structure width or pore size distribution. It displays a side-by-side comparison of the original data and the computed local thickness (heat map), where color intensity represents the diameter of the largest sphere that fits inside the structure at that point. It also includes a histogram to quantify the distribution of thickness values.

For 3D volumes, the output can be either an interactive slice viewer or a static Maximum Intensity Projection (MIP).

Parameters:

Name	Type	Description	Default
image	ndarray	The original 2D or 3D input data (binary or grayscale).	required
image_lt	ndarray	The computed local thickness map (must have the same shape as image). This is typically the output of qim3d.processing.local_thickness.	required
max_projection	bool	If True (and input is 3D), collapses the volume along the specified axis using maximum projection before plotting. Results in a static 2D figure. Defaults to False.	False
axis	int	The axis along which to slice or project the volume. Defaults to 0.	0
slice_index	int or float	The initial slice to display for 3D volumes. int: The exact index of the slice. float: A fraction between 0.0 and 1.0 (e.g., 0.5 for the middle). None: Defaults to the middle slice.	None
show	bool	If True, explicitly calls plt.show() to render the plot immediately. Defaults to False.	False
figsize	tuple[int, int]	The width and height of the figure in inches. Defaults to (15, 5).	(15, 5)

Returns:

Name	Type	Description
object	interactive or Figure	The visualization object, depending on the input and parameters: widgets.interactive: Returned if the input is 3D and max_projection=False. Contains a slider for slice navigation. matplotlib.figure.Figure: Returned if the input is 2D or if max_projection=True.

Raises:

Type	Description
ValueError	If slice_index is a float outside the range [0, 1].

Example

import qim3d

fly = qim3d.examples.fly_150x256x256
lt_fly = qim3d.processing.local_thickness(fly)
qim3d.viz.local_thickness(fly, lt_fly, axis=0)

Source code in qim3d/viz/_local_thickness.py

def local_thickness(
    image: np.ndarray,
    image_lt: np.ndarray,
    max_projection: bool = False,
    axis: int = 0,
    slice_index: Optional[Union[int, float]] = None,
    show: bool = False,
    figsize: Tuple[int, int] = (15, 5),
) -> Union[plt.Figure, widgets.interactive]:
    """
    Visualizes a local thickness map alongside the original image and a statistics histogram.

    This function provides a comprehensive view of structure width or pore size distribution. It displays a side-by-side comparison of the original data and the computed local thickness (heat map), where color intensity represents the diameter of the largest sphere that fits inside the structure at that point. It also includes a histogram to quantify the distribution of thickness values.

    For 3D volumes, the output can be either an interactive slice viewer or a static Maximum Intensity Projection (MIP).

    Args:
        image (np.ndarray): The original 2D or 3D input data (binary or grayscale).
        image_lt (np.ndarray): The computed local thickness map (must have the same shape as `image`). This is typically the output of `qim3d.processing.local_thickness`.
        max_projection (bool, optional): If `True` (and input is 3D), collapses the volume along the specified axis using maximum projection before plotting. Results in a static 2D figure. Defaults to `False`.
        axis (int, optional): The axis along which to slice or project the volume. Defaults to 0.
        slice_index (int or float, optional): The initial slice to display for 3D volumes.

            * **int**: The exact index of the slice.
            * **float**: A fraction between 0.0 and 1.0 (e.g., 0.5 for the middle).
            * **None**: Defaults to the middle slice.

        show (bool, optional): If `True`, explicitly calls `plt.show()` to render the plot immediately. Defaults to `False`.
        figsize (tuple[int, int], optional): The width and height of the figure in inches. Defaults to (15, 5).

    Returns:
        object (widgets.interactive or matplotlib.figure.Figure):
            The visualization object, depending on the input and parameters:

            * **widgets.interactive**: Returned if the input is 3D and `max_projection=False`. Contains a slider for slice navigation.
            * **matplotlib.figure.Figure**: Returned if the input is 2D or if `max_projection=True`.

    Raises:
        ValueError: If `slice_index` is a float outside the range [0, 1].

    Example:
        ```python
        import qim3d

        fly = qim3d.examples.fly_150x256x256
        lt_fly = qim3d.processing.local_thickness(fly)
        qim3d.viz.local_thickness(fly, lt_fly, axis=0)
        ```
        ![local thickness 3d](../../assets/screenshots/local_thickness_3d.gif)
    """

    def _local_thickness(image, image_lt, show, figsize, axis=None, slice_index=None):
        if slice_index is not None:
            image = image.take(slice_index, axis=axis)
            image_lt = image_lt.take(slice_index, axis=axis)

        fig, axs = plt.subplots(1, 3, figsize=figsize, layout='constrained')

        axs[0].imshow(image, cmap='gray')
        axs[0].set_title('Original image')
        axs[0].axis('off')

        axs[1].imshow(image_lt, cmap='viridis')
        axs[1].set_title('Local thickness')
        axs[1].axis('off')

        plt.colorbar(
            axs[1].imshow(image_lt, cmap='viridis'), ax=axs[1], orientation='vertical'
        )

        axs[2].hist(image_lt[image_lt > 0].ravel(), bins=32, edgecolor='black')
        axs[2].set_title('Local thickness histogram')
        axs[2].set_xlabel('Local thickness')
        axs[2].set_ylabel('Count')

        if show:
            plt.show()

        plt.close()

        return fig

    # Get the middle slice if the input is 3D
    if len(image.shape) == 3:
        if max_projection:
            if slice_index is not None:
                log.warning(
                    'slice_index is not used for max_projection. It will be ignored.'
                )
            image = image.max(axis=axis)
            image_lt = image_lt.max(axis=axis)
            return _local_thickness(image, image_lt, show, figsize)
        else:
            if slice_index is None:
                slice_index = image.shape[axis] // 2
            elif isinstance(slice_index, float):
                if slice_index < 0 or slice_index > 1:
                    raise ValueError(
                        'Values of slice_index of float type must be between 0 and 1.'
                    )
                slice_index = int(slice_index * image.shape[0]) - 1
            slice_index_slider = widgets.IntSlider(
                min=0,
                max=image.shape[axis] - 1,
                step=1,
                value=slice_index,
                description='Slice index',
                layout=widgets.Layout(width='450px'),
            )
            widget_obj = widgets.interactive(
                _local_thickness,
                image=widgets.fixed(image),
                image_lt=widgets.fixed(image_lt),
                show=widgets.fixed(True),
                figsize=widgets.fixed(figsize),
                axis=widgets.fixed(axis),
                slice_index=slice_index_slider,
            )
            widget_obj.layout = widgets.Layout(align_items='center')
            if show:
                display(widget_obj)
            return widget_obj
    else:
        if max_projection:
            log.warning(
                'max_projection is only used for 3D images. It will be ignored.'
            )
        if slice_index is not None:
            log.warning('slice_index is only used for 3D images. It will be ignored.')
        return _local_thickness(image, image_lt, show, figsize)

qim3d.viz.vectors

vectors(
    volume,
    vectors,
    axis=0,
    volume_colormap='grey',
    min_value=None,
    max_value=None,
    slice_index=None,
    grid_size=10,
    interactive=True,
    figsize=(10, 5),
    show=False,
)

Visualizes the local orientation of structures using structure tensor eigenvectors.

This function is designed to analyze anisotropy and fiber direction in 3D volumes. It generates a three-panel visualization to interpret the structure tensor output:

1. Quiver Plot: Overlays vector arrows on the image slice to show the dominant direction in the 2D plane.
2. Orientation Histogram: displays the distribution of angles, helping to identify if structures are aligned or random.
3. Color Map: Renders the slice using an HSV scheme where Hue represents the in-plane angle and Saturation represents the out-of-plane alignment (vectors pointing out of the screen appear desaturated/gray).

Parameters:

Name	Type	Description	Default
volume	ndarray	The 3D input volume (scalar field).	required
vectors	ndarray	The eigenvectors of the structure tensor, typically shape (3, Z, Y, X).	required
axis	int	The axis along which to slice the volume for visualization (0, 1, or 2). Defaults to 0.	0
volume_colormap	str	The colormap for the background volume slice in the quiver plot. Defaults to 'grey'.	'grey'
min_value	float	Minimum value for the volume display contrast.	None
max_value	float	Maximum value for the volume display contrast.	None
slice_index	int or float	The initial slice to display. int: The exact index of the slice. float: A fraction between 0.0 and 1.0. None: Defaults to the middle slice.	None
grid_size	int	The spacing between arrows in the quiver plot. Lower values result in denser vector fields. Defaults to 10.	10
interactive	bool	If True, returns a widget with sliders for slicing and grid density. If False, returns a static figure. Defaults to True.	True
figsize	tuple[int, int]	The width and height of the figure in inches. Defaults to (10, 5).	(10, 5)
show	bool	If True, calls plt.show() to display the plot immediately. Defaults to False.	False

Returns:

Name	Type	Description
object	interactive or Figure	The visualization object. widgets.interactive: Returned if interactive=True. matplotlib.figure.Figure: Returned if interactive=False.

Raises:

Type	Description
ValueError	If axis is invalid or slice_index is out of bounds.

Example

import qim3d

vol = qim3d.examples.NT_128x128x128
val, vec = qim3d.processing.structure_tensor(vol)

# Visualize the structure tensor
qim3d.viz.vectors(vol, vec, axis = 2, interactive = True)

Source code in qim3d/viz/_structure_tensor.py

def vectors(
    volume: np.ndarray,
    vectors: np.ndarray,
    axis: int = 0,
    volume_colormap: str = 'grey',
    min_value: float | None = None,
    max_value: float | None = None,
    slice_index: int | float | None = None,
    grid_size: int = 10,
    interactive: bool = True,
    figsize: tuple[int, int] = (10, 5),
    show: bool = False,
) -> plt.Figure | widgets.interactive:
    """
    Visualizes the local orientation of structures using structure tensor eigenvectors.

    This function is designed to analyze anisotropy and fiber direction in 3D volumes. It generates a three-panel visualization to interpret the structure tensor output:

    1.  **Quiver Plot:** Overlays vector arrows on the image slice to show the dominant direction in the 2D plane.
    2.  **Orientation Histogram:** displays the distribution of angles, helping to identify if structures are aligned or random.
    3.  **Color Map:** Renders the slice using an HSV scheme where Hue represents the in-plane angle and Saturation represents the out-of-plane alignment (vectors pointing out of the screen appear desaturated/gray).

    Args:
        volume (np.ndarray): The 3D input volume (scalar field).
        vectors (np.ndarray): The eigenvectors of the structure tensor, typically shape `(3, Z, Y, X)`.
        axis (int, optional): The axis along which to slice the volume for visualization (0, 1, or 2). Defaults to 0.
        volume_colormap (str, optional): The colormap for the background volume slice in the quiver plot. Defaults to 'grey'.
        min_value (float, optional): Minimum value for the volume display contrast.
        max_value (float, optional): Maximum value for the volume display contrast.
        slice_index (int or float, optional): The initial slice to display.

            * **int**: The exact index of the slice.
            * **float**: A fraction between 0.0 and 1.0.
            * **None**: Defaults to the middle slice.

        grid_size (int, optional): The spacing between arrows in the quiver plot. Lower values result in denser vector fields. Defaults to 10.
        interactive (bool, optional): If `True`, returns a widget with sliders for slicing and grid density. If `False`, returns a static figure. Defaults to `True`.
        figsize (tuple[int, int], optional): The width and height of the figure in inches. Defaults to (10, 5).
        show (bool, optional): If `True`, calls `plt.show()` to display the plot immediately. Defaults to `False`.

    Returns:
        object (widgets.interactive or matplotlib.figure.Figure):
            The visualization object.

            * **widgets.interactive**: Returned if `interactive=True`.
            * **matplotlib.figure.Figure**: Returned if `interactive=False`.

    Raises:
        ValueError: If `axis` is invalid or `slice_index` is out of bounds.

    Example:
        ```python
        import qim3d

        vol = qim3d.examples.NT_128x128x128
        val, vec = qim3d.processing.structure_tensor(vol)

        # Visualize the structure tensor
        qim3d.viz.vectors(vol, vec, axis = 2, interactive = True)
        ```
        ![structure tensor](../../assets/screenshots/structure_tensor_visualization.gif)

    """

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

    # Normalize the volume if needed (i.e. if values are in [0, 255])
    if volume.max() > 1.0:
        volume = volume / 255.0

    # Define grid size limits
    min_grid_size = max(1, volume.shape[axis] // 50)
    max_grid_size = max(1, volume.shape[axis] // 10)
    if max_grid_size <= min_grid_size:
        max_grid_size = min_grid_size * 5

    if not grid_size:
        grid_size = (min_grid_size + max_grid_size) // 2

    # Testing
    if grid_size < min_grid_size or grid_size > max_grid_size:
        # Adjust grid size as little as possible to be within the limits
        grid_size = min(max(min_grid_size, grid_size), max_grid_size)
        log.warning(f'Adjusting grid size to {grid_size} as it is out of bounds.')

    def _structure_tensor(volume, vectors, axis, slice_index, grid_size, figsize, show):
        # Choose the appropriate slice based on the specified dimension
        if axis == 0:
            data_slice = volume[slice_index, :, :]
            vectors_slice_x = vectors[0, slice_index, :, :]
            vectors_slice_y = vectors[1, slice_index, :, :]
            vectors_slice_z = vectors[2, slice_index, :, :]

        elif axis == 1:
            data_slice = volume[:, slice_index, :]
            vectors_slice_x = vectors[0, :, slice_index, :]
            vectors_slice_y = vectors[2, :, slice_index, :]
            vectors_slice_z = vectors[1, :, slice_index, :]

        elif axis == 2:
            data_slice = volume[:, :, slice_index]
            vectors_slice_x = vectors[1, :, :, slice_index]
            vectors_slice_y = vectors[2, :, :, slice_index]
            vectors_slice_z = vectors[0, :, :, slice_index]

        else:
            msg = 'Invalid dimension. Use 0 for Z, 1 for Y, or 2 for X.'
            raise ValueError(msg)

        # Create three subplots
        fig, ax = plt.subplots(1, 3, figsize=figsize, layout='constrained')

        blend_hue_saturation = (
            lambda hue, sat: hue * (1 - sat) + 0.5 * sat
        )  # Function for blending hue and saturation
        blend_slice_colors = lambda slice, colors: 0.5 * (
            slice + colors
        )  # Function for blending image slice with orientation colors

        # ----- Subplot 1: Image slice with orientation vectors ----- #
        # Create meshgrid with the correct dimensions
        xmesh, ymesh = np.mgrid[0 : data_slice.shape[0], 0 : data_slice.shape[1]]

        # Create a slice object for selecting the grid points
        g = slice(grid_size // 2, None, grid_size)  # noqa: A002

        # Angles from 0 to pi
        angles_quiver = np.mod(
            np.arctan2(
                vectors_slice_y[g, g],
                vectors_slice_x[g, g],
            ),
            np.pi,
        )

        # Calculate z-component (saturation)
        saturation_quiver = (vectors_slice_z[g, g] ** 2)[:, :, np.newaxis]

        # Calculate hue
        hue_quiver = plt.cm.hsv(angles_quiver / np.pi)

        # Blend hue and saturation
        rgba_quiver = blend_hue_saturation(hue_quiver, saturation_quiver)
        rgba_quiver = np.clip(
            rgba_quiver, 0, 1
        )  # Ensure rgba values are values within [0, 1]
        rgba_quiver_flat = rgba_quiver.reshape(
            (rgba_quiver.shape[0] * rgba_quiver.shape[1], 4)
        )  # Flatten array for quiver plot

        # Plot vectors
        ax[0].quiver(
            ymesh[g, g],
            xmesh[g, g],
            vectors_slice_x[g, g],
            vectors_slice_y[g, g],
            color=rgba_quiver_flat,
            angles='xy',
        )
        ax[0].quiver(
            ymesh[g, g],
            xmesh[g, g],
            -vectors_slice_x[g, g],
            -vectors_slice_y[g, g],
            color=rgba_quiver_flat,
            angles='xy',
        )

        ax[0].imshow(data_slice, cmap=volume_colormap, vmin=min_value, vmax=max_value)
        ax[0].set_title(
            f'Orientation vectors (slice {slice_index})'
            if not interactive
            else 'Orientation vectors'
        )
        ax[0].set_axis_off()

        # ----- Subplot 2: Orientation histogram ----- #
        nbins = 36

        # Angles from 0 to pi
        angles = np.mod(np.arctan2(vectors_slice_y, vectors_slice_x), np.pi)

        # Orientation histogram over angles
        distribution, bin_edges = np.histogram(angles, bins=nbins, range=(0.0, np.pi))

        # Half circle (180 deg)
        bin_centers = (np.arange(nbins) + 0.5) * np.pi / nbins

        # Calculate z-component (saturation) for each bin
        bins = np.digitize(angles.ravel(), bin_edges)
        saturation_bin = np.array(
            [
                (
                    np.mean((vectors_slice_z**2).ravel()[bins == i])
                    if np.sum(bins == i) > 0
                    else 0
                )
                for i in range(1, len(bin_edges))
            ]
        )

        # Calculate hue for each bin
        hue_bin = plt.cm.hsv(bin_centers / np.pi)

        # Blend hue and saturation
        rgba_bin = hue_bin.copy()
        rgba_bin[:, :3] = blend_hue_saturation(
            hue_bin[:, :3], saturation_bin[:, np.newaxis]
        )

        ax[1].bar(bin_centers, distribution, width=np.pi / nbins, color=rgba_bin)
        ax[1].set_xlabel('Angle [radians]')
        ax[1].set_xlim([0, np.pi])
        ax[1].set_aspect(np.pi / ax[1].get_ylim()[1])
        ax[1].set_xticks([0, np.pi / 2, np.pi])
        ax[1].set_xticklabels(['0', '$\\frac{\\pi}{2}$', '$\\pi$'])
        ax[1].set_yticks([])
        ax[1].set_ylabel('Frequency')
        ax[1].set_title('Histogram over orientation angles')

        # ----- Subplot 3: Image slice colored according to orientation ----- #
        # Calculate z-component (saturation)
        saturation = (vectors_slice_z**2)[:, :, np.newaxis]

        # Calculate hue
        hue = plt.cm.hsv(angles / np.pi)

        # Blend hue and saturation
        rgba = blend_hue_saturation(hue, saturation)

        # Grayscale image slice blended with orientation colors
        data_slice_orientation_colored = (
            blend_slice_colors(plt.cm.gray(data_slice), rgba) * 255
        ).astype('uint8')

        ax[2].imshow(data_slice_orientation_colored)
        ax[2].set_title(
            f'Colored orientations (slice {slice_index})'
            if not interactive
            else 'Colored orientations'
        )
        ax[2].set_axis_off()

        if show:
            plt.show()

        plt.close()

        return fig

    if vectors.ndim == 5:
        vectors = vectors[0, ...]
        log.warning(
            'Eigenvector array is full. Only the eigenvectors corresponding to the first eigenvalue will be used.'
        )

    if slice_index is None:
        slice_index = volume.shape[axis] // 2

    elif isinstance(slice_index, float):
        if slice_index < 0 or slice_index > 1:
            raise ValueError(
                'Values of slice_index of float type must be between 0 and 1.'
            )
        slice_index = int(slice_index * volume.shape[0]) - 1

    if interactive:
        slice_index_slider = widgets.IntSlider(
            min=0,
            max=volume.shape[axis] - 1,
            step=1,
            value=slice_index,
            description='Slice index',
            layout=widgets.Layout(width='450px'),
        )

        grid_size_slider = widgets.IntSlider(
            min=min_grid_size,
            max=max_grid_size,
            step=1,
            value=grid_size,
            description='Grid size',
            layout=widgets.Layout(width='450px'),
        )

        widget_obj = widgets.interactive(
            _structure_tensor,
            volume=widgets.fixed(volume),
            vectors=widgets.fixed(vectors),
            axis=widgets.fixed(axis),
            slice_index=slice_index_slider,
            grid_size=grid_size_slider,
            figsize=widgets.fixed(figsize),
            show=widgets.fixed(True),
        )
        # Arrange sliders horizontally
        sliders_box = widgets.HBox([slice_index_slider, grid_size_slider])
        widget_obj = widgets.VBox([sliders_box, widget_obj.children[-1]])
        widget_obj.layout.align_items = 'center'

        if show:
            display(widget_obj)

        return widget_obj

    else:
        return _structure_tensor(
            volume, vectors, axis, slice_index, grid_size, figsize, show
        )

qim3d.viz.vector_field_3d

vector_field_3d(
    vec,
    val,
    select_eigen='smallest',
    sampling_step=4,
    max_cones=50000,
    cone_size=1,
    verbose=True,
    colormap='Portland',
    cmin=None,
    cmax=None,
    **kwargs,
)

Generates an interactive 3D quiver plot (vector field) using cones to visualize local orientation.

This function is designed to visualize the output of structure tensor analysis. It draws 3D cones representing the eigenvectors at various points in the volume. This is widely used to analyze material anisotropy, such as identifying fiber orientations (linear structures) or surface normals (planar structures).

To handle large 3D datasets efficiently, the function downsamples the volume by averaging vectors within grid blocks (sampling_step) and filters to show only the most significant regions (max_cones).

Parameters:

Name	Type	Description	Default
vec	ndarray	The eigenvectors from the structure tensor. Expected shapes: (3, Z, Y, X): Contains a single eigenvector per voxel (e.g., result of structure_tensor(..., smallest=True)). (3, 3, Z, Y, X): Contains all three eigenvectors per voxel.	required
val	ndarray	The eigenvalues from the structure tensor, with shape (3, Z, Y, X).	required
select_eigen	str	Determines which eigenvector to visualize (only used if vec contains all three). 'smallest': Visualizes the vector associated with the smallest eigenvalue (typically the surface normal in planar structures). 'largest': Visualizes the vector associated with the largest eigenvalue (typically the fiber direction in linear structures). 'middle': Visualizes the intermediate vector.	'smallest'
sampling_step	int	The stride for the sampling grid. A value of 4 means every 4th voxel in each dimension is sampled (averaging the region). Higher values improve performance but reduce resolution. Defaults to 4.	4
max_cones	int	The maximum number of cones to render. If the sampled grid exceeds this, only the points with the highest eigenvalues (strongest features) are kept. Defaults to 50000.	50000
cone_size	float	A scaling factor for the size of the cones. Defaults to 1.	1
verbose	bool	If True, prints progress and grid statistics. Defaults to True.	True
colormap	str	The name of the Plotly colorscale. Defaults to 'Portland'.	'Portland'
cmin	float	The minimum value for the color mapping. If None, uses the data minimum.	None
cmax	float	The maximum value for the color mapping. If None, uses the data maximum.	None
**kwargs		Additional keyword arguments passed to plotly.graph_objects.Cone. See the Plotly Cone documentation for full customization options.	{}

Returns:

Name	Type	Description
fig	Figure	The interactive Plotly figure containing the 3D vector field.

Raises:

Type	Description
ValueError	If select_eigen is not one of 'smallest', 'largest', or 'middle'.

Example

vol = qim3d.examples.fibers_150x150x150
val, vec = qim3d.processing.structure_tensor(vol, smallest = True)
qim3d.viz.vector_field_3d(vec, val, sampling_step=12, max_cones=5000, cone_size = 2)

Notes

Interpreting Eigenvalues and Eigenvectors

The structure tensor yields three eigenvalues (λ₁ ≤ λ₂ ≤ λ₃) and corresponding eigenvectors. The smallest eigenvector (v₁) indicates the predominant orientation, the direction of minimum intensity variation.

Structure Type	Eigenvalue Pattern	Interpretation
Linear (fibers, tubes)	λ₁ ≪ λ₂ ≈ λ₃	Follow v₁ to trace fiber
Planar (sheets, boundaries)	λ₁ ≈ λ₂ ≪ λ₃	v₁ tangent to surface
Isotropic (noise, blobs)	λ₁ ≈ λ₂ ≈ λ₃	No predominant direction

Structure Tensor Notes

Source code in qim3d/viz/_structure_tensor.py

def vector_field_3d(
    vec: np.ndarray,
    val: np.ndarray,
    select_eigen: Literal['smallest', 'largest', 'middle'] = 'smallest',
    sampling_step: int = 4,
    max_cones: int = 50000,
    cone_size: float = 1,
    verbose: bool = True,
    colormap: str = 'Portland',
    cmin: float = None,
    cmax: float = None,
    **kwargs,
) -> go.Figure:
    """
    Generates an interactive 3D quiver plot (vector field) using cones to visualize local orientation.

    This function is designed to visualize the output of structure tensor analysis. It draws 3D cones representing the eigenvectors at various points in the volume. This is widely used to analyze material anisotropy, such as identifying fiber orientations (linear structures) or surface normals (planar structures).

    To handle large 3D datasets efficiently, the function downsamples the volume by averaging vectors within grid blocks (`sampling_step`) and filters to show only the most significant regions (`max_cones`).

    Args:
        vec (np.ndarray): The eigenvectors from the structure tensor.
            Expected shapes:

            * `(3, Z, Y, X)`: Contains a single eigenvector per voxel (e.g., result of `structure_tensor(..., smallest=True)`).
            * `(3, 3, Z, Y, X)`: Contains all three eigenvectors per voxel.

        val (np.ndarray): The eigenvalues from the structure tensor, with shape `(3, Z, Y, X)`.
        select_eigen (str, optional): Determines which eigenvector to visualize (only used if `vec` contains all three).

            * `'smallest'`: Visualizes the vector associated with the smallest eigenvalue (typically the **surface normal** in planar structures).
            * `'largest'`: Visualizes the vector associated with the largest eigenvalue (typically the **fiber direction** in linear structures).
            * `'middle'`: Visualizes the intermediate vector.

        sampling_step (int, optional): The stride for the sampling grid. A value of `4` means every 4th voxel in each dimension is sampled (averaging the region). Higher values improve performance but reduce resolution. Defaults to 4.
        max_cones (int, optional): The maximum number of cones to render. If the sampled grid exceeds this, only the points with the highest eigenvalues (strongest features) are kept. Defaults to 50000.
        cone_size (float, optional): A scaling factor for the size of the cones. Defaults to 1.
        verbose (bool, optional): If `True`, prints progress and grid statistics. Defaults to `True`.
        colormap (str, optional): The name of the Plotly colorscale. Defaults to 'Portland'.
        cmin (float, optional): The minimum value for the color mapping. If `None`, uses the data minimum.
        cmax (float, optional): The maximum value for the color mapping. If `None`, uses the data maximum.
        **kwargs: Additional keyword arguments passed to `plotly.graph_objects.Cone`. See the [Plotly Cone documentation](https://plotly.com/python-api-reference/generated/plotly.graph_objects.Cone.html)
            for full customization options.

    Returns:
        fig (plotly.graph_objects.Figure):
            The interactive Plotly figure containing the 3D vector field.

    Raises:
        ValueError: If `select_eigen` is not one of 'smallest', 'largest', or 'middle'.

    Example:
        ```python
        vol = qim3d.examples.fibers_150x150x150
        val, vec = qim3d.processing.structure_tensor(vol, smallest = True)
        qim3d.viz.vector_field_3d(vec, val, sampling_step=12, max_cones=5000, cone_size = 2)
        ```
        ![vector field](../../assets/screenshots/viz-vector_field.png)


    Notes:
        **Interpreting Eigenvalues and Eigenvectors**

        The structure tensor yields three eigenvalues (λ₁ ≤ λ₂ ≤ λ₃) and corresponding eigenvectors.
        The smallest eigenvector (v₁) indicates the predominant orientation, the direction of minimum
        intensity variation.

        | Structure Type | Eigenvalue Pattern | Interpretation |
        | :--- | :--- | :--- |
        | **Linear** (fibers, tubes) | λ₁ ≪ λ₂ ≈ λ₃ | Follow v₁ to trace fiber |
        | **Planar** (sheets, boundaries) | λ₁ ≈ λ₂ ≪ λ₃ | v₁ tangent to surface  |
        | **Isotropic** (noise, blobs) | λ₁ ≈ λ₂ ≈ λ₃ | No predominant direction |

        [Structure Tensor Notes](https://people.compute.dtu.dk/vand/notes/ST_intro.pdf)

    """
    if vec.ndim == 4:
        val = val[0]
    elif vec.ndim == 5:
        if select_eigen == 'largest':
            val = val[2]
            vec = vec[2, :, ...]
        elif select_eigen == 'smallest':
            val = val[0]
            vec = vec[0, :, ...]
        elif select_eigen == 'middle':
            val = val[1]
            vec = vec[1, :, ...]
        else:
            msg = f'Invalid select_eigen value: {select_eigen}. Choose between "smallest", "largest", or "middle".'
            raise ValueError(msg)
    vec = np.transpose(vec, (1, 2, 3, 0))

    nx, ny, nz, _ = vec.shape
    if verbose:
        log.info(f'Original number of grid points: {nx * ny * nz}')
    half = sampling_step // 2

    # Sampling grid
    grid_x = np.arange(0, nx, sampling_step)
    grid_y = np.arange(0, ny, sampling_step)
    grid_z = np.arange(0, nz, sampling_step)

    points, vectors, values = [], [], []

    # Average vectors and eigenvalues in each sampling cube
    for px in grid_x:
        for py in grid_y:
            for pz in grid_z:
                x0, x1 = max(px - half, 0), min(px + half + 1, nx)
                y0, y1 = max(py - half, 0), min(py + half + 1, ny)
                z0, z1 = max(pz - half, 0), min(pz + half + 1, nz)

                region_vecs = vec[x0:x1, y0:y1, z0:z1, :]
                region_vals = val[x0:x1, y0:y1, z0:z1]

                avg_vec = region_vecs.mean(axis=(0, 1, 2))
                avg_val = region_vals.mean()

                points.append((px, py, pz))
                vectors.append(avg_vec)
                values.append(avg_val)

    points = np.array(points)
    vectors = np.array(vectors)
    values = np.array(values)

    # Select top N highest eigenvalue locations
    idx_top = np.argsort(values)[::-1][:max_cones]
    points_top = points[idx_top]
    vectors_top = vectors[idx_top]
    values_top = values[idx_top]

    if verbose:
        log.info(f'Number of grid points sampled: {len(values)}')
        log.info(f'Number of cones actually plotted: {len(points_top)}')

    # Normalize vectors and scale by eigenvalue magnitude
    norms = np.linalg.norm(vectors_top, axis=1, keepdims=True)
    norms[norms == 0] = 1
    unit_vecs = vectors_top / norms

    # Apply decay to downscale weak directions
    # scaled_strength = values_top * np.exp(
    #     -decay_rate * (1 - values_top / values_top.max())
    # )
    # scaled_strength = np.log(values_top - values_top.min() + 1)

    scaled_strength = np.log1p(values_top)

    u = unit_vecs[:, 0] * scaled_strength
    v = unit_vecs[:, 1] * scaled_strength
    w = unit_vecs[:, 2] * scaled_strength

    # Compute magnitude for color scaling if needed
    magnitude = np.sqrt(u**2 + v**2 + w**2)
    min_mag = magnitude.min()
    max_mag = magnitude.max()
    if verbose:
        log.info(f'Min magnitude: {min_mag:.4f}, Max magnitude: {max_mag:.4f}')

    if cmin is None:
        cmin = min_mag
    if cmax is None:
        cmax = max_mag

    # Use 'raw' to display sizes in actual vector length.
    fig = go.Figure(
        data=go.Cone(
            x=points_top[:, 0],
            y=points_top[:, 1],
            z=points_top[:, 2],
            u=u,
            v=v,
            w=w,
            sizemode='scaled',
            sizeref=cone_size,
            colorscale=colormap,
            colorbar_title='Orientation strength',
            cmin=cmin,
            cmax=cmax,
            **kwargs,
        ),
        layout={'width': 900, 'height': 700},
    )

    return fig