diff --git a/docs/assets/screenshots/pygel3d_visualization.png b/docs/assets/screenshots/pygel3d_visualization.png
new file mode 100644
index 0000000000000000000000000000000000000000..93d87dde129906bfaa1944ded8f5e0df7783251c
Binary files /dev/null and b/docs/assets/screenshots/pygel3d_visualization.png differ
diff --git a/docs/doc/cli/cli.md b/docs/doc/cli/cli.md
index 58fae7e2253830a9ab3156712fe7fac14d2c3f24..123fee26161db2f28737d6302cccb10520b9d88a 100644
--- a/docs/doc/cli/cli.md
+++ b/docs/doc/cli/cli.md
@@ -44,7 +44,7 @@ The command line interface allows you to run graphical user interfaces directly
 
 !!! Example
 
-    Here's an example of how to open the [Data Explorer](gui.md#qim3d.gui.data_explorer)
+    Here's an example of how to open the [Data Explorer](../gui/gui.md#qim3d.gui.data_explorer)
 
     ``` title="Command"
     qim3d gui --data-explorer
diff --git a/docs/doc/data_handling/io.md b/docs/doc/data_handling/io.md
index 82c984bbced233b2a35836bdd072e3c73aa8bf8b..607d81ede891e5700eb593ab7ddfa876a6037ac0 100644
--- a/docs/doc/data_handling/io.md
+++ b/docs/doc/data_handling/io.md
@@ -8,5 +8,5 @@
             - Downloader
             - export_ome_zarr
             - import_ome_zarr
-            - save_mesh
-            - load_mesh
\ No newline at end of file
+            - load_mesh
+            - save_mesh
\ No newline at end of file
diff --git a/docs/doc/gui/gui.md b/docs/doc/gui/gui.md
index 681041fc4f2515ce6c89112b07c12d1d407f44f3..4bc2de99a9e81559313a703df31a1ef40be4c811 100644
--- a/docs/doc/gui/gui.md
+++ b/docs/doc/gui/gui.md
@@ -21,7 +21,7 @@ The `qim3d` library provides a set of custom made GUIs that ease the interaction
     ```
 
 In general, the GUIs can be launched directly from the command line. 
-For details see [here](cli.md#qim3d-gui).
+For details see [here](../cli/cli.md#qim3d-gui).
 
 ::: qim3d.gui.data_explorer
     options:
diff --git a/docs/doc/releases/releases.md b/docs/doc/releases/releases.md
index 8174ea97c706b522de00460a71f1cdf7c59da2f8..bdb5c2a8007d46a3622a5215bfb2231752d45dbc 100644
--- a/docs/doc/releases/releases.md
+++ b/docs/doc/releases/releases.md
@@ -11,7 +11,7 @@ hide:
 
 Below, you'll find details about the version history of `qim3d`.
 
-Remember to keep your pip installation [up to date](index.md/#get-the-latest-version) so that you have the latest features!
+Remember to keep your pip installation [up to date](../../index.md/#get-the-latest-version) so that you have the latest features!
 
 ### v1.0.0 (21/01/2025)
 
diff --git a/qim3d/__init__.py b/qim3d/__init__.py
index 51d18ccba10a49d30b257dac43d37a73da247adf..5857fad7f1ed0bb1493b9a622be948cac277aea2 100644
--- a/qim3d/__init__.py
+++ b/qim3d/__init__.py
@@ -9,7 +9,7 @@ Documentation available at https://platform.qim.dk/qim3d/
 
 """
 
-__version__ = '1.0.0'
+__version__ = '1.1.0'
 
 
 import importlib as _importlib
diff --git a/qim3d/features/_common_features_methods.py b/qim3d/features/_common_features_methods.py
index 3f797ee998e0c278a1674b581c8a3b6898e76d02..4a547b6c30c730e7c0d0ab583f3f9995d91df747 100644
--- a/qim3d/features/_common_features_methods.py
+++ b/qim3d/features/_common_features_methods.py
@@ -1,22 +1,19 @@
 import numpy as np
-import trimesh
-
 import qim3d
-import qim3d.processing
 from qim3d.utils._logger import log
+import qim3d
+from pygel3d import hmesh
 
-
-def volume(obj: np.ndarray | trimesh.Trimesh, **mesh_kwargs) -> float:
+def volume(obj: np.ndarray|hmesh.Manifold) -> float:
     """
-    Compute the volume of a 3D volume or mesh.
+    Compute the volume of a 3D mesh using the Pygel3D library.
 
     Args:
-        obj (np.ndarray or trimesh.Trimesh): Either a np.ndarray volume or a mesh object of type trimesh.Trimesh.
-        **mesh_kwargs (Any): Additional arguments for mesh creation if the input is a volume.
+        obj (numpy.ndarray or pygel3d.hmesh.Manifold): Either a np.ndarray volume or a mesh object of type pygel3d.hmesh.Manifold.
 
     Returns:
         volume (float): The volume of the object.
-
+    
     Example:
         Compute volume from a mesh:
         ```python
@@ -26,8 +23,8 @@ def volume(obj: np.ndarray | trimesh.Trimesh, **mesh_kwargs) -> float:
         mesh = qim3d.io.load_mesh('path/to/mesh.obj')
 
         # Compute the volume of the mesh
-        vol = qim3d.features.volume(mesh)
-        print('Volume:', vol)
+        volume = qim3d.features.volume(mesh)
+        print(f'Volume: {volume}')
         ```
 
         Compute volume from a np.ndarray:
@@ -36,33 +33,30 @@ def volume(obj: np.ndarray | trimesh.Trimesh, **mesh_kwargs) -> float:
 
         # Generate a 3D blob
         synthetic_blob = qim3d.generate.noise_object(noise_scale = 0.015)
-        synthetic_blob = qim3d.generate.noise_object(noise_scale = 0.015)
 
         # Compute the volume of the blob
-        volume = qim3d.features.volume(synthetic_blob, level=0.5)
-        volume = qim3d.features.volume(synthetic_blob, level=0.5)
-        print('Volume:', volume)
+        volume = qim3d.features.volume(synthetic_blob)
+        print(f'Volume: {volume}')
         ```
 
     """
-    if isinstance(obj, np.ndarray):
-        log.info('Converting volume to mesh.')
-        obj = qim3d.mesh.from_volume(obj, **mesh_kwargs)
 
-    return obj.volume
+    if isinstance(obj, np.ndarray):
+        log.info("Converting volume to mesh.")
+        obj = qim3d.mesh.from_volume(obj)
 
+    return hmesh.volume(obj)
 
-def area(obj: np.ndarray | trimesh.Trimesh, **mesh_kwargs) -> float:
+def area(obj: np.ndarray|hmesh.Manifold) -> float:
     """
-    Compute the surface area of a 3D volume or mesh.
+    Compute the surface area of a 3D mesh using the Pygel3D library.
 
     Args:
-        obj (np.ndarray or trimesh.Trimesh): Either a np.ndarray volume or a mesh object of type trimesh.Trimesh.
-        **mesh_kwargs (Any): Additional arguments for mesh creation if the input is a volume.
+        obj (numpy.ndarray or pygel3d.hmesh.Manifold): Either a np.ndarray volume or a mesh object of type pygel3d.hmesh.Manifold.
 
     Returns:
-        area (float): The surface area of the object.
-
+        area (float): The surface area of the object. 
+    
     Example:
         Compute area from a mesh:
         ```python
@@ -73,8 +67,7 @@ def area(obj: np.ndarray | trimesh.Trimesh, **mesh_kwargs) -> float:
 
         # Compute the surface area of the mesh
         area = qim3d.features.area(mesh)
-        area = qim3d.features.area(mesh)
-        print(f"Area: {area}")
+        print(f'Area: {area}')
         ```
 
         Compute area from a np.ndarray:
@@ -83,38 +76,30 @@ def area(obj: np.ndarray | trimesh.Trimesh, **mesh_kwargs) -> float:
 
         # Generate a 3D blob
         synthetic_blob = qim3d.generate.noise_object(noise_scale = 0.015)
-        synthetic_blob = qim3d.generate.noise_object(noise_scale = 0.015)
 
         # Compute the surface area of the blob
-        volume = qim3d.features.area(synthetic_blob, level=0.5)
-        volume = qim3d.features.area(synthetic_blob, level=0.5)
-        print('Area:', volume)
+        area = qim3d.features.area(synthetic_blob)
+        print(f'Area: {area}')
         ```
-
+    
     """
-    if isinstance(obj, np.ndarray):
-        log.info('Converting volume to mesh.')
-        obj = qim3d.mesh.from_volume(obj, **mesh_kwargs)
-        obj = qim3d.mesh.from_volume(obj, **mesh_kwargs)
 
-    return obj.area
+    if isinstance(obj, np.ndarray):
+        log.info("Converting volume to mesh.")
+        obj = qim3d.mesh.from_volume(obj)
 
+    return hmesh.area(obj)
 
-def sphericity(obj: np.ndarray | trimesh.Trimesh, **mesh_kwargs) -> float:
+def sphericity(obj: np.ndarray|hmesh.Manifold) -> float:
     """
-    Compute the sphericity of a 3D volume or mesh.
-
-    Sphericity is a measure of how spherical an object is. It is defined as the ratio
-    of the surface area of a sphere with the same volume as the object to the object's
-    actual surface area.
+    Compute the sphericity of a 3D mesh using the Pygel3D library.
 
     Args:
-        obj (np.ndarray or trimesh.Trimesh): Either a np.ndarray volume or a mesh object of type trimesh.Trimesh.
-        **mesh_kwargs (Any): Additional arguments for mesh creation if the input is a volume.
+        obj (numpy.ndarray or pygel3d.hmesh.Manifold): Either a np.ndarray volume or a mesh object of type pygel3d.hmesh.Manifold.
 
     Returns:
-        sphericity (float): A float value representing the sphericity of the object.
-
+        sphericity (float): The sphericity of the object. 
+    
     Example:
         Compute sphericity from a mesh:
         ```python
@@ -125,7 +110,7 @@ def sphericity(obj: np.ndarray | trimesh.Trimesh, **mesh_kwargs) -> float:
 
         # Compute the sphericity of the mesh
         sphericity = qim3d.features.sphericity(mesh)
-        sphericity = qim3d.features.sphericity(mesh)
+        print(f'Sphericity: {sphericity}')
         ```
 
         Compute sphericity from a np.ndarray:
@@ -134,33 +119,25 @@ def sphericity(obj: np.ndarray | trimesh.Trimesh, **mesh_kwargs) -> float:
 
         # Generate a 3D blob
         synthetic_blob = qim3d.generate.noise_object(noise_scale = 0.015)
-        synthetic_blob = qim3d.generate.noise_object(noise_scale = 0.015)
 
         # Compute the sphericity of the blob
-        sphericity = qim3d.features.sphericity(synthetic_blob, level=0.5)
-        sphericity = qim3d.features.sphericity(synthetic_blob, level=0.5)
+        sphericity = qim3d.features.sphericity(synthetic_blob)
+        print(f'Sphericity: {sphericity}')
         ```
 
-    !!! info "Limitations due to pixelation"
-        Sphericity is particularly sensitive to the resolution of the mesh, as it directly impacts the accuracy of surface area and volume calculations.
-        Since the mesh is generated from voxel-based 3D volume data, the discrete nature of the voxels leads to pixelation effects that reduce the precision of sphericity measurements.
-        Higher resolution meshes may mitigate these errors but often at the cost of increased computational demands.
-
     """
+
     if isinstance(obj, np.ndarray):
-        log.info('Converting volume to mesh.')
-        obj = qim3d.mesh.from_volume(obj, **mesh_kwargs)
-        obj = qim3d.mesh.from_volume(obj, **mesh_kwargs)
+        log.info("Converting volume to mesh.")
+        obj = qim3d.mesh.from_volume(obj)
 
     volume = qim3d.features.volume(obj)
     area = qim3d.features.area(obj)
-    volume = qim3d.features.volume(obj)
-    area = qim3d.features.area(obj)
 
     if area == 0:
         log.warning('Surface area is zero, sphericity is undefined.')
         return np.nan
 
     sphericity = (np.pi ** (1 / 3) * (6 * volume) ** (2 / 3)) / area
-    log.info(f'Sphericity: {sphericity}')
-    return sphericity
+    # log.info(f"Sphericity: {sphericity}")
+    return sphericity
\ No newline at end of file
diff --git a/qim3d/filters/_common_filter_methods.py b/qim3d/filters/_common_filter_methods.py
index 548d6b73f8840c33a9395ca1fbded83e811a8d35..c550fef0b1ca7a2f1ffc74cd3dd19af5c678baa6 100644
--- a/qim3d/filters/_common_filter_methods.py
+++ b/qim3d/filters/_common_filter_methods.py
@@ -334,7 +334,6 @@ def gaussian(
         sigma (float or sequence of floats): The standard deviations of the Gaussian filter are given for each axis as a sequence, or as a single number, in which case it is equal for all axes.
         dask (bool, optional): Whether to use Dask for the Gaussian filter.
         chunks (int or tuple or "'auto'", optional): 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 (Any): Additional positional arguments for the Gaussian filter.
         **kwargs (Any): Additional keyword arguments for the Gaussian filter.
 
     Returns:
diff --git a/qim3d/io/_loading.py b/qim3d/io/_loading.py
index b221aef80829ca62c27cc75faf48c11bd2cdeb73..cd041a1f03d348134eaa39f80451c4e1a00d120d 100644
--- a/qim3d/io/_loading.py
+++ b/qim3d/io/_loading.py
@@ -27,8 +27,11 @@ import qim3d
 from qim3d.utils import Memory, log
 from qim3d.utils._misc import get_file_size, sizeof, stringify_path
 from qim3d.utils._progress_bar import FileLoadingProgressBar
+import trimesh
+from pygel3d import hmesh
+from typing import Optional, Dict
 
-dask.config.set(scheduler='processes')
+dask.config.set(scheduler="processes") 
 
 
 class DataLoader:
@@ -891,15 +894,23 @@ def load(
     return data
 
 
-def load_mesh(filename: str) -> trimesh.Trimesh:
+def load_mesh(filename: str) -> hmesh.Manifold:
     """
-    Load a mesh from an .obj file using trimesh.
+    Load a mesh from a specific file.
+    This function is based on the [PyGEL3D library's loading function implementation](https://www2.compute.dtu.dk/projects/GEL/PyGEL/pygel3d/hmesh.html#load).
+
+    Supported formats:
+
+    - `X3D`
+    - `OBJ`
+    - `OFF`
+    - `PLY`
 
     Args:
-        filename (str or os.PathLike): The path to the .obj file.
+        filename (str or os.PathLike): The path to the file.
 
     Returns:
-        mesh (trimesh.Trimesh): A trimesh object containing the mesh data (vertices and faces).
+        mesh (hmesh.Manifold or None): A hmesh object containing the mesh data or None if loading failed.
 
     Example:
         ```python
@@ -909,5 +920,6 @@ def load_mesh(filename: str) -> trimesh.Trimesh:
         ```
 
     """
-    mesh = trimesh.load(filename)
-    return mesh
+    mesh = hmesh.load(filename)
+
+    return mesh
\ No newline at end of file
diff --git a/qim3d/io/_saving.py b/qim3d/io/_saving.py
index 6e5a8dabd240da19a8ed7acc63334224bbd2f1be..cde1c358f041b65826572c7dc1f5d0a69896b506 100644
--- a/qim3d/io/_saving.py
+++ b/qim3d/io/_saving.py
@@ -37,7 +37,8 @@ import trimesh
 import zarr
 from pydicom.dataset import FileDataset, FileMetaDataset
 from pydicom.uid import UID
-
+import trimesh
+from pygel3d import hmesh
 from qim3d.utils import log
 from qim3d.utils._misc import stringify_path
 
@@ -484,29 +485,54 @@ def save(
     ).save(path, data)
 
 
-def save_mesh(filename: str, mesh: trimesh.Trimesh) -> None:
+# def save_mesh(
+#         filename: str, 
+#         mesh: trimesh.Trimesh
+#         ) -> None:
+#     """
+#     Save a trimesh object to an .obj file.
+
+#     Args:
+#         filename (str or os.PathLike): The name of the file to save the mesh.
+#         mesh (trimesh.Trimesh): A trimesh.Trimesh object representing the mesh.
+
+#     Example:
+#         ```python
+#         import qim3d
+
+#         vol = qim3d.generate.noise_object(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.mesh.from_volume(vol)
+#         qim3d.io.save_mesh("mesh.obj", mesh)
+#         ```
+#     """
+#     # Export the mesh to the specified filename
+#     mesh.export(filename)
+    
+def save_mesh(filename: str, mesh: hmesh.Manifold) -> None:
     """
-    Save a trimesh object to an .obj file.
+    Save a mesh object to a specific file.
+    This function is based on the [PyGEL3D library's saving function implementation](https://www2.compute.dtu.dk/projects/GEL/PyGEL/pygel3d/hmesh.html#save).
+
 
     Args:
-        filename (str or os.PathLike): The name of the file to save the mesh.
-        mesh (trimesh.Trimesh): A trimesh.Trimesh object representing the mesh.
+        filename (str or os.PathLike): The path to save file to. File format is chosen based on the extension. Supported extensions are: '.x3d', '.obj', '.off'.
+        mesh (pygel3d.hmesh.Manifold): A hmesh.Manifold object representing the mesh.
 
     Example:
         ```python
         import qim3d
 
-        vol = qim3d.generate.noise_object(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.mesh.from_volume(vol)
-        qim3d.io.save_mesh("mesh.obj", mesh)
+        synthetic_blob = qim3d.generate.noise_object(noise_scale = 0.015)
+        mesh = qim3d.mesh.from_volume(synthetic_blob)
+        qim3d.io.save_mesh("mesh.obj", mesh) 
         ```
 
     """
     # Export the mesh to the specified filename
-    mesh.export(filename)
+    hmesh.save(filename, mesh)
\ No newline at end of file
diff --git a/qim3d/mesh/_common_mesh_methods.py b/qim3d/mesh/_common_mesh_methods.py
index 883629d98ff74d7a03e8faf7b15bc04cbada4261..9e1c5e814b1704ec6f0a95764e1f7b16285b9570 100644
--- a/qim3d/mesh/_common_mesh_methods.py
+++ b/qim3d/mesh/_common_mesh_methods.py
@@ -1,85 +1,46 @@
 from typing import Any, Tuple
 
 import numpy as np
-import trimesh
-from skimage import filters, measure
-
+from skimage import measure, filters
+from pygel3d import hmesh
+from typing import Tuple, Any
 from qim3d.utils._logger import log
 
 
 def from_volume(
     volume: np.ndarray,
-    level: float = None,
-    step_size: int = 1,
-    allow_degenerate: bool = False,
-    padding: Tuple[int, int, int] = (2, 2, 2),
-    **kwargs: Any,
-) -> trimesh.Trimesh:
-    """
-    Convert a volume to a mesh using the Marching Cubes algorithm, with optional thresholding and padding.
+    **kwargs: any
+) -> hmesh.Manifold:
+    """ Convert a 3D numpy array to a mesh object using the [volumetric_isocontour](https://www2.compute.dtu.dk/projects/GEL/PyGEL/pygel3d/hmesh.html#volumetric_isocontour) function from Pygel3D.
 
     Args:
-        volume (np.ndarray): The 3D numpy array representing the volume.
-        level (float, optional): The threshold value for Marching Cubes. If None, Otsu's method is used.
-        step_size (int, optional): The step size for the Marching Cubes algorithm.
-        allow_degenerate (bool, optional): Whether to allow degenerate (i.e. zero-area) triangles in the end-result. If False, degenerate triangles are removed, at the cost of making the algorithm slower. Default False.
-        padding (tuple of ints, optional): Padding to add around the volume.
-        **kwargs: Additional keyword arguments to pass to `skimage.measure.marching_cubes`.
+        volume (np.ndarray): A 3D numpy array representing a volume.
+        **kwargs: Additional arguments to pass to the Pygel3D volumetric_isocontour function.
+
+    Raises:
+        ValueError: If the input volume is not a 3D numpy array or if the input volume is empty.
 
     Returns:
-        mesh (trimesh.Trimesh): The generated mesh.
+        hmesh.Manifold: A Pygel3D mesh object representing the input volume.
 
     Example:
+        Convert a 3D numpy array to a Pygel3D mesh object:
         ```python
         import qim3d
-        vol = qim3d.generate.noise_object(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.mesh.from_volume(vol, step_size=3)
-        qim3d.viz.mesh(mesh.vertices, mesh.faces)
-        ```
-        <iframe src="https://platform.qim.dk/k3d/mesh_visualization.html" width="100%" height="500" frameborder="0"></iframe>
 
+        # Generate a 3D blob
+        synthetic_blob = qim3d.generate.noise_object(noise_scale = 0.015)
+
+        # Convert the 3D numpy array to a Pygel3D mesh object
+        mesh = qim3d.mesh.from_volume(synthetic_blob)
+        ```
     """
+    
     if volume.ndim != 3:
-        raise ValueError('The input volume must be a 3D numpy array.')
-
-    # Compute the threshold level if not provided
-    if level is None:
-        level = filters.threshold_otsu(volume)
-        log.info(f"Computed level using Otsu's method: {level}")
-
-    # Apply padding to the volume
-    if padding is not None:
-        pad_z, pad_y, pad_x = padding
-        padding_value = np.min(volume)
-        volume = np.pad(
-            volume,
-            ((pad_z, pad_z), (pad_y, pad_y), (pad_x, pad_x)),
-            mode='constant',
-            constant_values=padding_value,
-        )
-        log.info(f'Padded volume with {padding} to shape: {volume.shape}')
-
-    # Call skimage.measure.marching_cubes with user-provided kwargs
-    verts, faces, normals, values = measure.marching_cubes(
-        volume,
-        level=level,
-        step_size=step_size,
-        allow_degenerate=allow_degenerate,
-        **kwargs,
-    )
-
-    # Create the Trimesh object
-    mesh = trimesh.Trimesh(vertices=verts, faces=faces)
-
-    # Fix face orientation to ensure normals point outwards
-    trimesh.repair.fix_inversion(mesh, multibody=True)
+        raise ValueError("The input volume must be a 3D numpy array.")
+    
+    if volume.size == 0:
+        raise ValueError("The input volume must not be empty.")
 
-    return mesh
+    mesh = hmesh.volumetric_isocontour(volume, **kwargs)
+    return mesh
\ No newline at end of file
diff --git a/qim3d/tests/mesh/test_mesh.py b/qim3d/tests/mesh/test_mesh.py
new file mode 100644
index 0000000000000000000000000000000000000000..9e19475bfd7adbbd0d11dd4e96c03a15ca2e970a
--- /dev/null
+++ b/qim3d/tests/mesh/test_mesh.py
@@ -0,0 +1,42 @@
+import pytest
+import numpy as np
+from pygel3d import hmesh
+import qim3d
+
+def test_from_volume_valid_input():
+    """Test that from_volume returns a hmesh.Manifold object for a valid 3D input."""
+    volume = np.random.rand(50, 50, 50).astype(np.float32)  # Generate a random 3D volume
+    mesh = qim3d.mesh.from_volume(volume)
+    assert isinstance(mesh, hmesh.Manifold)  # Check if output is a Manifold object
+
+def test_from_volume_invalid_input():
+    """Test that from_volume raises ValueError for non-3D input."""
+    volume = np.random.rand(50, 50)  # A 2D array
+    with pytest.raises(ValueError, match="The input volume must be a 3D numpy array."):
+        qim3d.mesh.from_volume(volume)
+
+def test_from_volume_empty_array():
+    """Test how from_volume handles an empty 3D array."""
+    volume = np.empty((0, 0, 0))  # Empty 3D array
+    with pytest.raises(ValueError):  # It should fail because it doesn't make sense to generate a mesh from empty data
+        qim3d.mesh.from_volume(volume)
+
+def test_from_volume_with_kwargs():
+    """Test that from_volume correctly passes kwargs."""
+    volume = np.random.rand(50, 50, 50).astype(np.float32)
+
+    # Mock volumetric_isocontour to check if kwargs are passed
+    def mock_volumetric_isocontour(vol, **kwargs):
+        assert "isovalue" in kwargs
+        assert kwargs["isovalue"] == 0.5
+        return hmesh.Manifold()
+
+    # Replace the function temporarily
+    original_function = hmesh.volumetric_isocontour
+    hmesh.volumetric_isocontour = mock_volumetric_isocontour
+
+    try:
+        qim3d.mesh.from_volume(volume, isovalue=0.5)
+    finally:
+        hmesh.volumetric_isocontour = original_function  # Restore original function
+
diff --git a/qim3d/viz/_k3d.py b/qim3d/viz/_k3d.py
index 448a79b927bbd91e732089b826691c0f312e7054..5d1ecb5532187fbabe2217e2b5eb826a03e7e843 100644
--- a/qim3d/viz/_k3d.py
+++ b/qim3d/viz/_k3d.py
@@ -13,7 +13,10 @@ from matplotlib.colors import Colormap
 
 from qim3d.utils._logger import log
 from qim3d.utils._misc import downscale_img, scale_to_float16
-
+from pygel3d import hmesh
+from pygel3d import jupyter_display as jd
+import k3d
+from typing import Optional
 
 def volumetric(
     img: np.ndarray,
@@ -82,7 +85,6 @@ def volumetric(
         ```
 
     """
-    import k3d
 
     pixel_count = img.shape[0] * img.shape[1] * img.shape[2]
     # target is 60fps on m1 macbook pro, using test volume: https://data.qim.dk/pages/foam.html
@@ -172,88 +174,106 @@ def volumetric(
     else:
         return plot
 
-
 def mesh(
-    verts: np.ndarray,
-    faces: np.ndarray,
+    mesh,
+    backend: str = "pygel3d",
     wireframe: bool = True,
     flat_shading: bool = True,
     grid_visible: bool = False,
     show: bool = True,
     save: bool = False,
     **kwargs,
-):
-    """
-    Visualizes a 3D mesh using K3D.
-
+)-> Optional[k3d.Plot]:
+    """Visualize a 3D mesh using `pygel3d` or `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.
+        mesh (pygel3d.hmesh.Manifold): The input mesh object.
+        backend (str, optional): The visualization backend to use. 
+            Choose between `pygel3d` (default) and `k3d`.
+        wireframe (bool, optional): If True, displays the mesh as a wireframe.
+            Works both with `pygel3d` and `k3d`. Defaults to True.
+        flat_shading (bool, optional): If True, applies flat shading to the mesh.
+            Works only with `k3d`. Defaults to True.
+        grid_visible (bool, optional): If True, shows a grid in the visualization.
+            Works only with `k3d`. Defaults to False.
+        show (bool, optional): If True, displays the visualization inline.
+            Works only with `k3d`. 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 (Any): Additional keyword arguments to be passed to the `k3d.plot` function.
-
+            file will be saved. Works only with `k3d`. Defaults to False.
+        **kwargs (Any): Additional keyword arguments specific to the chosen backend:
+            
+            - `k3d.plot` kwargs: Arguments that customize the [`k3d.plot`](https://k3d-jupyter.org/reference/factory.plot.html) visualization. 
+            
+            - `pygel3d.display` kwargs: Arguments that customize the [`pygel3d.display`](https://www2.compute.dtu.dk/projects/GEL/PyGEL/pygel3d/jupyter_display.html#display) visualization.
+    
     Returns:
-        plot (k3d.plot): If `show=False`, returns the K3D plot object.
+        k3d.Plot or None:
+        
+            - If `backend="k3d"`, returns a `k3d.Plot` object.
+            - If `backend="pygel3d"`, the function displays the mesh but does not return a plot object.
+    
+    Raises:
+        ValueError: If `backend` is not `pygel3d` or `k3d`.
 
     Example:
         ```python
         import qim3d
-
-        vol = qim3d.generate.noise_object(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.mesh.from_volume(vol, step_size=3)
-        qim3d.viz.mesh(mesh.vertices, mesh.faces)
+        synthetic_blob = qim3d.generate.noise_object(noise_scale = 0.015)
+        mesh = qim3d.mesh.from_volume(synthetic_blob)
+        qim3d.viz.mesh(mesh, backend="pygel3d") # or qim3d.viz.mesh(mesh, backend="k3d")
         ```
-        <iframe src="https://platform.qim.dk/k3d/mesh_visualization.html" width="100%" height="500" frameborder="0"></iframe>
-
+    ![pygel3d_visualization](../../assets/screenshots/pygel3d_visualization.png)
     """
-    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 backend not in ["k3d", "pygel3d"]:
+        raise ValueError("Invalid backend. Choose 'pygel3d' or 'k3d'.")
 
-    if show:
-        plot.display()
-    else:
-        return plot
+    # Extract vertex positions and face indices
+    face_indices = list(mesh.faces())
+    vertices_array = np.array(mesh.positions())
+
+    # Extract face vertex indices
+    face_vertices = [
+        list(mesh.circulate_face(int(fid), mode="v"))[:3] for fid in face_indices
+    ]
+    face_vertices = np.array(face_vertices, dtype=np.uint32)
+
+    # Validate the mesh structure
+    if vertices_array.shape[1] != 3 or face_vertices.shape[1] != 3:
+        raise ValueError("Vertices must have shape (N, 3) and faces (M, 3)")
+
+    # Separate valid kwargs for each backend
+    valid_k3d_kwargs = {k: v for k, v in kwargs.items() if k not in ["smooth", "data"]}
+    valid_pygel_kwargs = {k: v for k, v in kwargs.items() if k in ["smooth", "data"]}
+
+    if backend == "k3d":
+        vertices_array = np.ascontiguousarray(vertices_array.astype(np.float32))
+        face_vertices = np.ascontiguousarray(face_vertices)
+
+        mesh_plot = k3d.mesh(
+            vertices=vertices_array,
+            indices=face_vertices,
+            wireframe=wireframe,
+            flat_shading=flat_shading,
+        )
+
+        # Create plot
+        plot = k3d.plot(grid_visible=grid_visible, **valid_k3d_kwargs)
+        plot += mesh_plot
+
+        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
+
+
+    elif backend == "pygel3d":
+        jd.set_export_mode(True)
+        return jd.display(mesh, wireframe=wireframe, **valid_pygel_kwargs)
diff --git a/requirements.txt b/requirements.txt
index 7ec31e7cb5fa53cc6ada37eccbe79fbce9f4b1c2..c17889959dbab43dd7ada03e9adac084c65f86b0 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -31,3 +31,4 @@ ome_zarr>=0.9.0
 dask-image>=2024.5.3
 trimesh>=4.4.9
 slgbuilder>=0.2.1
+PyGEL3D>=0.5.2
\ No newline at end of file
diff --git a/setup.py b/setup.py
index 883da55c9b9caa42d808ee75874f58a24e341b9c..1da74828d168f805e4ad3d530be16c5865b87244 100644
--- a/setup.py
+++ b/setup.py
@@ -71,7 +71,8 @@ setup(
         "ome_zarr>=0.9.0",
         "dask-image>=2024.5.3",
         "scikit-image>=0.24.0",
-        "trimesh>=4.4.9"
+        "trimesh>=4.4.9",
+        "PyGEL3D>=0.5.2"
     ],
     extras_require={
     "deep-learning": [