docs/io.md

@@ -7,5 +7,9 @@ Currently, it is possible to directly load `tiff`, `h5`, `nii`,`txm`, `vol` and
     options:
         members:
             - load
+            - load_mesh
             - save
+            - save_mesh
             - Downloader
+            - export_ome_zarr
+            - import_ome_zarr
\ No newline at end of file

docs/processing.md

@@ -14,6 +14,7 @@ Here, we provide functionalities designed specifically for 3D image analysis and
             - maximum
             - minimum
             - tophat
+            - create_mesh
 
 ::: qim3d.processing.Pipeline
     options:

docs/releases.md

@@ -9,6 +9,17 @@ As the library is still in its early development stages, **there may be breaking
 
 And remember to keep your pip installation [up to date](/qim3d/#get-the-latest-version) so that you have the latest features!
 
+### v0.4.3 (02/10/2024)
+
+- Updated requirements
+- Introduction of mesh generation, visualization, saving and loading 🎉  ![Mesh generation](assets/screenshots/releases/mesh_generation.png)
+
+### v0.4.2 (30/09/2024)
+
+- Export and import is now possible in the OME-Zarr standard, including multi-scale datasets.
+- Filters now have the option to use Dask when available ![Filters with dask](assets/screenshots/releases/filters_with_dask.png)
+
+
 ### v0.4.1 (30/07/2024)
 
 - Fixed issue with example volumes not being loaded

docs/viz.md

@@ -8,6 +8,7 @@ The `qim3d` library aims to provide easy ways to explore and get insights from v
             - slicer
             - orthogonal
             - vol
+            - mesh
             - local_thickness
             - vectors
             - plot_cc

qim3d/__init__.py

@@ -8,7 +8,7 @@ Documentation available at https://platform.qim.dk/qim3d/
 
 """
 
-__version__ = "0.4.1"
+__version__ = "0.4.3"
 
 
 import importlib as _importlib

qim3d/io/__init__.py

-from .loading import DataLoader, load
+from .loading import DataLoader, load, load_mesh
 from .downloader import Downloader
-from .saving import DataSaver, save 
+from .saving import DataSaver, save, save_mesh
 from .sync import Sync
 from .convert import convert
 from ..utils import logger
+from .ome_zarr import export_ome_zarr, import_ome_zarr

qim3d/io/loading.py

@@ -27,6 +27,7 @@ from qim3d.utils.logger import log
 from qim3d.utils.misc import get_file_size, sizeof, stringify_path
 from qim3d.utils.system import Memory
 from qim3d.utils import ProgressBar
+import trimesh
 
 dask.config.set(scheduler="processes") 
 
@@ -861,3 +862,22 @@ def load(
 
     return data
 
+def load_mesh(filename):
+    """
+    Load a mesh from an .obj file using trimesh.
+
+    Args:
+        filename: The path to the .obj file.
+
+    Returns:
+        mesh: A trimesh object containing the mesh data (vertices and faces).
+
+    Example:
+        ```python
+        import qim3d
+
+        mesh = qim3d.io.load_mesh("path/to/mesh.obj")
+        ```
+    """
+    mesh = trimesh.load(filename)
+    return mesh

qim3d/io/ome_zarr.py

+"""
+Exporting data to different formats.
+"""
+
+import os
+import numpy as np
+import zarr
+from ome_zarr.io import parse_url
+from ome_zarr.writer import (
+    write_image,
+    _create_mip,
+    write_multiscale,
+    CurrentFormat,
+    Format,
+)
+from ome_zarr.scale import dask_resize
+from ome_zarr.reader import Reader
+from ome_zarr import scale
+
+import math
+import shutil
+from qim3d.utils.logger import log
+from scipy.ndimage import zoom
+from typing import Any, Callable, Iterator, List, Tuple, Union
+import dask.array as da
+from skimage.transform import (
+    resize,
+)
+
+
+ListOfArrayLike = Union[List[da.Array], List[np.ndarray]]
+ArrayLike = Union[da.Array, np.ndarray]
+
+
+class OMEScaler(
+    scale.Scaler,
+):
+    """Scaler in the style of OME-Zarr.
+    This is needed because their current zoom implementation is broken."""
+
+    def __init__(self, order=0, downscale=2, max_layer=5, method="scaleZYXdask"):
+        self.order = order
+        self.downscale = downscale
+        self.max_layer = max_layer
+        self.method = method
+
+    def scaleZYX(self, base):
+        """Downsample using :func:`scipy.ndimage.zoom`."""
+        rv = [base]
+        log.info(f"- Scale 0: {rv[-1].shape}")
+
+        for i in range(self.max_layer):
+            downscale_ratio = (1 / self.downscale) ** (i + 1)
+            rv.append(zoom(base, zoom=downscale_ratio, order=self.order))
+            log.info(f"- Scale {i+1}: {rv[-1].shape}")
+        return list(rv)
+
+    def scaleZYXdask(self, base):
+        """Downsample using :func:`scipy.ndimage.zoom`."""
+        rv = [base]
+        log.info(f"- Scale 0: {rv[-1].shape}")
+
+        for i in range(self.max_layer):
+
+            scaled_shape = tuple(
+                base.shape[j] // (self.downscale ** (i + 1)) for j in range(3)
+            )
+            rv.append(dask_resize(base, scaled_shape, order=self.order))
+
+            log.info(f"- Scale {i+1}: {rv[-1].shape}")
+        return list(rv)
+
+
+def export_ome_zarr(
+    path,
+    data,
+    chunk_size=100,
+    downsample_rate=2,
+    order=0,
+    replace=False,
+    method="scaleZYX",
+):
+    """
+    Export image data to OME-Zarr format with pyramidal downsampling.
+
+    Automatically calculates the number of downsampled scales such that the smallest scale fits within the specified `chunk_size`.
+
+    Args:
+        path (str): The directory where the OME-Zarr data will be stored.
+        data (np.ndarray): The image data to be exported.
+        chunk_size (int, optional): The size of the chunks for storing data. Defaults to 100.
+        downsample_rate (int, optional): Factor by which to downsample the data for each scale. Must be greater than 1. Defaults to 2.
+        order (int, optional): Interpolation order to use when downsampling. Defaults to 0 (nearest-neighbor).
+        replace (bool, optional): Whether to replace the existing directory if it already exists. Defaults to False.
+
+    Raises:
+        ValueError: If the directory already exists and `replace` is False.
+        ValueError: If `downsample_rate` is less than or equal to 1.
+
+    Example:
+        ```python
+        import qim3d
+
+        downloader = qim3d.io.Downloader()
+        data = downloader.Snail.Escargot(load_file=True)
+
+        qim3d.io.export_ome_zarr("Escargot.zarr", data, chunk_size=100, downsample_rate=2)
+
+        ```
+
+    """
+
+    # Check if directory exists
+    if os.path.exists(path):
+        if replace:
+            shutil.rmtree(path)
+        else:
+            raise ValueError(
+                f"Directory {path} already exists. Use replace=True to overwrite."
+            )
+
+    # Check if downsample_rate is valid
+    if downsample_rate <= 1:
+        raise ValueError("Downsample rate must be greater than 1.")
+
+    log.info(f"Exporting data to OME-Zarr format at {path}")
+
+    # Get the number of scales
+    min_dim = np.max(np.shape(data))
+    nscales = math.ceil(math.log(min_dim / chunk_size) / math.log(downsample_rate))
+    log.info(f"Number of scales: {nscales + 1}")
+
+    # Create scaler
+    scaler = OMEScaler(
+        downscale=downsample_rate, max_layer=nscales, method=method, order=order
+    )
+
+    # write the image data
+    os.mkdir(path)
+    store = parse_url(path, mode="w").store
+    root = zarr.group(store=store)
+
+    fmt = CurrentFormat()
+    log.info("Creating a multi-scale pyramid")
+    mip, axes = _create_mip(image=data, fmt=fmt, scaler=scaler, axes="zyx")
+
+    log.info("Writing data to disk")
+    write_multiscale(
+        mip,
+        group=root,
+        fmt=fmt,
+        axes=axes,
+        name=None,
+        compute=True,
+    )
+
+    log.info("All done!")
+    return
+
+
+def import_ome_zarr(path, scale=0, load=True):
+    """
+    Import image data from an OME-Zarr file.
+
+    This function reads OME-Zarr formatted volumetric image data and returns the specified scale.
+    The image data can be lazily loaded (as Dask arrays) or fully computed into memory.
+
+    Args:
+        path (str): The file path to the OME-Zarr data.
+        scale (int or str, optional): The scale level to load.
+            If 'highest', loads the finest scale (scale 0).
+            If 'lowest', loads the coarsest scale (last available scale). Defaults to 0.
+        load (bool, optional): Whether to compute the selected scale into memory.
+            If False, returns a lazy Dask array. Defaults to True.
+
+    Returns:
+        np.ndarray or dask.array.Array: The requested image data, either as a NumPy array if `load=True`,
+        or a Dask array if `load=False`.
+
+    Raises:
+        ValueError: If the requested `scale` does not exist in the data.
+
+    Example:
+        ```python
+        import qim3d
+
+        data = qim3d.io.import_ome_zarr("Escargot.zarr", scale=0, load=True)
+
+        ```
+
+    """
+
+    # read the image data
+    # store = parse_url(path, mode="r").store
+
+    reader = Reader(parse_url(path))
+    nodes = list(reader())
+    image_node = nodes[0]
+    dask_data = image_node.data
+
+    log.info(f"Data contains {len(dask_data)} scales:")
+    for i in np.arange(len(dask_data)):
+        log.info(f"- Scale {i}: {dask_data[i].shape}")
+
+    if scale == "highest":
+        scale = 0
+
+    if scale == "lowest":
+        scale = len(dask_data) - 1
+
+    if scale >= len(dask_data):
+        raise ValueError(
+            f"Scale {scale} does not exist in the data. Please choose a scale between 0 and {len(dask_data)-1}."
+        )
+
+    log.info(f"\nLoading scale {scale} with shape {dask_data[scale].shape}")
+
+    if load:
+        vol = dask_data[scale].compute()
+    else:
+        vol = dask_data[scale]
+
+    return vol

qim3d/io/saving.py

@@ -24,7 +24,6 @@ Example:
 
 import datetime
 import os
-
 import dask.array as da
 import h5py
 import nibabel as nib
@@ -34,6 +33,7 @@ import tifffile
 import zarr
 from pydicom.dataset import FileDataset, FileMetaDataset
 from pydicom.uid import UID
+import trimesh
 
 from qim3d.utils.logger import log
 from qim3d.utils.misc import sizeof, stringify_path
@@ -77,6 +77,7 @@ class DataSaver:
         self.compression = kwargs.get("compression", False)
         self.basename = kwargs.get("basename", None)
         self.sliced_dim = kwargs.get("sliced_dim", 0)
+        self.chunk_shape = kwargs.get("chunk_shape", "auto")
 
     def save_tiff(self, path, data):
         """Save data to a TIFF file to the given path.
@@ -267,11 +268,24 @@ class DataSaver:
         Returns:
             zarr.core.Array: The Zarr array saved on disk.
         """
-        assert isinstance(data, da.Array), 'data must be a dask array'
 
-        # forces compute when saving to zarr
-        da.to_zarr(data, path, compute=True, overwrite=self.replace, compressor=zarr.Blosc(cname='zstd', clevel=3, shuffle=2))
+        if isinstance(data, da.Array):
+            # If the data is a Dask array, save using dask
+            if self.chunk_shape:
+                log.info("Rechunking data to shape %s", self.chunk_shape)
+                data = data.rechunk(self.chunk_shape)
+            log.info("Saving Dask array to Zarr array on disk")
+            da.to_zarr(data, path, overwrite=self.replace)
 
+        else:
+            zarr_array = zarr.open(
+                path,
+                mode="w",
+                shape=data.shape,
+                chunks=self.chunk_shape,
+                dtype=data.dtype,
+            )
+            zarr_array[:] = data
 
     def save_PIL(self, path, data):
         """Save data to a PIL file to the given path.
@@ -319,7 +333,15 @@ class DataSaver:
 
         # If path is an existing directory
         if isdir:
-            # If basename is provided
+            # Check if this is a Zarr directory
+            if ".zarr" in path:
+                if self.replace:
+                    return self.save_to_zarr(path, data)
+                if not self.replace:
+                    raise ValueError(
+                        "A Zarr directory with the provided path already exists. To replace it set 'replace=True'"
+                    )
+            # If basename is provided, user wants to save as tiff stack
             if self.basename:
                 # Save as tiff stack
                 return self.save_tiff_stack(path, data)
@@ -388,7 +410,14 @@ class DataSaver:
 
 
 def save(
-    path, data, replace=False, compression=False, basename=None, sliced_dim=0, **kwargs
+    path,
+    data,
+    replace=False,
+    compression=False,
+    basename=None,
+    sliced_dim=0,
+    chunk_shape="auto",
+    **kwargs,
 ):
     """Save data to a specified file path.
 
@@ -419,5 +448,33 @@ def save(
         compression=compression,
         basename=basename,
         sliced_dim=sliced_dim,
+        chunk_shape=chunk_shape,
         **kwargs,
     ).save(path, data)
+
+
+def save_mesh(filename, mesh):
+    """
+    Save a trimesh object to an .obj file.
+
+    Args:
+        filename: The name of the file to save the mesh.
+        mesh: A trimesh.Trimesh object representing the mesh.
+
+    Example:
+        ```python
+        import qim3d
+
+        vol = qim3d.generate.blob(base_shape=(32, 32, 32),
+                                  final_shape=(32, 32, 32),
+                                  noise_scale=0.05,
+                                  order=1,
+                                  gamma=1.0,
+                                  max_value=255,
+                                  threshold=0.5)
+        mesh = qim3d.processing.create_mesh(vol)
+        qim3d.io.save_mesh("mesh.obj", mesh)
+        ```
+    """
+    # Export the mesh to the specified filename
+    mesh.export(filename)
\ No newline at end of file

qim3d/processing/__init__.py

@@ -4,3 +4,4 @@ from .detection import blob_detection
 from .filters import *
 from .operations import *
 from .cc import get_3d_cc
+from .mesh import create_mesh

qim3d/processing/filters.py

@@ -5,6 +5,8 @@ from typing import Type, Union
 import numpy as np
 from scipy import ndimage
 from skimage import morphology
+import dask.array as da
+import dask_image.ndfilters as dask_ndfilters
 
 from qim3d.utils.logger import log
 
@@ -24,7 +26,7 @@ __all__ = [
 
 
 class FilterBase:
-    def __init__(self, *args, **kwargs):
+    def __init__(self, dask=False, chunks="auto", *args, **kwargs):
         """
         Base class for image filters.
 
@@ -33,9 +35,10 @@ class FilterBase:
             **kwargs: Additional keyword arguments for filter initialization.
         """
         self.args = args
+        self.dask = dask
+        self.chunks = chunks
         self.kwargs = kwargs
 
-
 class Gaussian(FilterBase):
     def __call__(self, input):
         """
@@ -47,7 +50,7 @@ class Gaussian(FilterBase):
         Returns:
             The filtered image or volume.
         """
-        return gaussian(input, *self.args, **self.kwargs)
+        return gaussian(input, dask=self.dask, chunks=self.chunks, *self.args, **self.kwargs)
 
 
 class Median(FilterBase):
@@ -61,7 +64,7 @@ class Median(FilterBase):
         Returns:
             The filtered image or volume.
         """
-        return median(input, **self.kwargs)
+        return median(input, dask=self.dask, chunks=self.chunks, **self.kwargs)
 
 
 class Maximum(FilterBase):
@@ -75,7 +78,7 @@ class Maximum(FilterBase):
         Returns:
             The filtered image or volume.
         """
-        return maximum(input, **self.kwargs)
+        return maximum(input, dask=self.dask, chunks=self.chunks, **self.kwargs)
 
 
 class Minimum(FilterBase):
@@ -89,7 +92,7 @@ class Minimum(FilterBase):
         Returns:
             The filtered image or volume.
         """
-        return minimum(input, **self.kwargs)
+        return minimum(input, dask=self.dask, chunks=self.chunks, **self.kwargs)
 
 class Tophat(FilterBase):
     def __call__(self, input):
@@ -102,7 +105,7 @@ class Tophat(FilterBase):
         Returns:
             The filtered image or volume.
         """
-        return tophat(input, **self.kwargs)
+        return tophat(input, dask=self.dask, chunks=self.chunks, **self.kwargs)
 
 
 class Pipeline:
@@ -121,7 +124,7 @@ class Pipeline:
         # Create filter pipeline
         pipeline = Pipeline(
             Median(size=5),
-            Gaussian(sigma=3)
+            Gaussian(sigma=3, dask = True)
         )
 
         # Append a third filter to the pipeline
@@ -183,7 +186,7 @@ class Pipeline:
 
             # Create filter pipeline
             pipeline = Pipeline(
-                Maximum(size=3)
+                Maximum(size=3, dask=True),
             )
 
             # Append a second filter to the pipeline
@@ -207,77 +210,125 @@ class Pipeline:
         return input
 
 
-def gaussian(vol, *args, **kwargs):
+def gaussian(vol, dask=False, chunks='auto', *args, **kwargs):
     """
-    Applies a Gaussian filter to the input volume using scipy.ndimage.gaussian_filter.
+    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.
     """
-    return ndimage.gaussian_filter(vol, *args, **kwargs)
+    
+    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
 
 
-def median(vol, **kwargs):
+def median(vol, dask=False, chunks='auto', **kwargs):
     """
-    Applies a median filter to the input volume using scipy.ndimage.median_filter.
+    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.
     """
-    return ndimage.median_filter(vol, **kwargs)
+    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
 
 
-def maximum(vol, **kwargs):
+def maximum(vol, dask=False, chunks='auto', **kwargs):
     """
-    Applies a maximum filter to the input volume using scipy.ndimage.maximum_filter.
+    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.
     """
-    return ndimage.maximum_filter(vol, **kwargs)
+    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
 
 
-def minimum(vol, **kwargs):
+def minimum(vol, dask=False, chunks='auto', **kwargs):
     """
-    Applies a minimum filter to the input volume using scipy.ndimage.mainimum_filter.
+    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.
     """
-    return ndimage.minimum_filter(vol, **kwargs)
+    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
 
-def tophat(vol, **kwargs):
+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.
+        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
+        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)

qim3d/processing/mesh.py

+import numpy as np
+from skimage import measure, filters
+import trimesh
+from typing import Tuple, Any
+
+
+def create_mesh(
+    volume: np.ndarray,
+    level: float = None,
+    step_size=1,
+    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.
+        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)
+        ```
+
+    """
+    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)
+        print(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,
+        )
+        print(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, **kwargs
+    )
+
+    print(len(verts))
+
+    # Create the Trimesh object
+    mesh = trimesh.Trimesh(vertices=verts, faces=faces)
+
+    return mesh

qim3d/processing/operations.py

@@ -8,7 +8,7 @@ def remove_background(
     median_filter_size: int = 2,
     min_object_radius: int = 3,
     background: str = "dark",
-    **median_kwargs
+    **median_kwargs,
 ) -> np.ndarray:
     """
     Remove background from a volume using a qim3d filters.
@@ -51,17 +51,22 @@ def remove_background(
     # Apply the pipeline to the volume
     return pipeline(vol)
 
-def watershed(
-        bin_vol: np.ndarray
-) -> tuple[np.ndarray, int]:
+
+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.
+        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 after segmentation, number of objects found.
+        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
@@ -90,7 +95,9 @@ def watershed(
     distance = scipy.ndimage.distance_transform_edt(bin_vol)
 
     # Find peak coordinates in distance transform
-    coords = skimage.feature.peak_local_max(distance, labels=bin_vol)
+    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)
@@ -100,7 +107,9 @@ def watershed(
     markers, _ = scipy.ndimage.label(mask)
 
     # Apply watershed segmentation
-    labeled_volume = skimage.segmentation.watershed(-distance, markers=markers, mask=bin_vol)
+    labeled_volume = skimage.segmentation.watershed(
+        -distance, markers=markers, mask=bin_vol
+    )
 
     # Extract number of objects found
     num_labels = len(np.unique(labeled_volume)) - 1
@@ -108,6 +117,7 @@ def watershed(
 
     return labeled_volume, num_labels
 
+
 def fade_mask(
     vol: np.ndarray,
     decay_rate: float = 10,
@@ -147,7 +157,9 @@ def fade_mask(
 
     """
     if 0 > axis or axis >= vol.ndim:
-        raise ValueError("Axis must be between 0 and the number of dimensions of the volume")
+        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
@@ -186,10 +198,9 @@ def fade_mask(
 
     return vol_faded
 
+
 def overlay_rgb_images(
-    background: np.ndarray, 
-    foreground: np.ndarray, 
-    alpha: float = 0.5
+    background: np.ndarray, foreground: np.ndarray, alpha: float = 0.5
 ) -> np.ndarray:
     """
     Overlay an RGB foreground onto an RGB background using alpha blending.

qim3d/utils/misc.py

@@ -223,7 +223,7 @@ def downscale_img(img, max_voxels=512**3):
     zoom_factor = (max_voxels / total_voxels) ** (1 / 3)
 
     # Downscale image
-    return zoom(img, zoom_factor)
+    return zoom(img, zoom_factor, order=0)
 
 
 def scale_to_float16(arr: np.ndarray):

qim3d/viz/__init__.py

@@ -7,7 +7,7 @@ from .explore import (
     slicer,
     slices,
 )
-from .k3d import vol
+from .k3d import vol, mesh
 from .local_thickness_ import local_thickness
 from .structure_tensor import vectors
 from .metrics import plot_metrics, grid_overview, grid_pred, vol_masked

qim3d/viz/k3d.py

@@ -22,7 +22,7 @@ def vol(
     grid_visible=False,
     cmap=None,
     samples="auto",
-    max_voxels=412**3,
+    max_voxels=512**3,
     data_type="scaled_float16",
     **kwargs,
 ):
@@ -97,7 +97,7 @@ def vol(
 
     if original_shape != new_shape:
         log.warning(
-            f"Downsampled image for visualization. From {original_shape} to {new_shape}"
+            f"Downsampled image for visualization, from {original_shape} to {new_shape}"
         )
 
     # Scale the image to float16 if needed
@@ -141,3 +141,84 @@ def vol(
         plot.display()
     else:
         return plot  
+
+
+def mesh(
+    verts,
+    faces,
+    wireframe=True,
+    flat_shading=True,
+    grid_visible=False,
+    show=True,
+    save=False,
+    **kwargs,
+):
+    """
+    Visualizes a 3D mesh using K3D.
+
+    Args:
+        verts (numpy.ndarray): A 2D array (Nx3) containing the vertices of the mesh.
+        faces (numpy.ndarray): A 2D array (Mx3) containing the indices of the mesh faces.
+        wireframe (bool, optional): If True, the mesh is rendered as a wireframe. Defaults to True.
+        flat_shading (bool, optional): If True, flat shading is applied to the mesh. Defaults to True.
+        grid_visible (bool, optional): If True, the grid is visible in the plot. Defaults to False.
+        show (bool, optional): If True, displays the visualization inline. Defaults to True.
+        save (bool or str, optional): If True, saves the visualization as an HTML file.
+            If a string is provided, it's interpreted as the file path where the HTML
+            file will be saved. Defaults to False.
+        **kwargs: Additional keyword arguments to be passed to the `k3d.plot` function.
+
+    Returns:
+        plot (k3d.plot): If `show=False`, returns the K3D plot object.
+
+    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)
+        ```
+    """
+    import k3d
+
+    # Validate the inputs
+    if verts.shape[1] != 3:
+        raise ValueError("Vertices array must have shape (N, 3)")
+    if faces.shape[1] != 3:
+        raise ValueError("Faces array must have shape (M, 3)")
+    
+    # Ensure the correct data types and memory layout
+    verts = np.ascontiguousarray(verts.astype(np.float32))  # Cast and ensure C-contiguous layout
+    faces = np.ascontiguousarray(faces.astype(np.uint32))    # Cast and ensure C-contiguous layout
+
+
+    # Create the mesh plot
+    plt_mesh = k3d.mesh(
+        vertices=verts,
+        indices=faces,
+        wireframe=wireframe,
+        flat_shading=flat_shading,
+    )
+
+    # Create plot
+    plot = k3d.plot(grid_visible=grid_visible, **kwargs)
+    plot += plt_mesh
+
+    if save:
+        # Save html to disk
+        with open(str(save), "w", encoding="utf-8") as fp:
+            fp.write(plot.get_snapshot())
+
+    if show:
+        plot.display()
+    else:
+        return plot

requirements.txt

@@ -26,3 +26,6 @@ psutil>=5.9.0
 structure-tensor>=0.2.1
 noise>=1.2.2
 zarr>=2.18.2
+ome_zarr>=0.9.0
+dask-image>=2024.5.3
+trimesh>=4.4.9