diff --git a/qim3d/features/_common_features_methods.py b/qim3d/features/_common_features_methods.py
index 62a7611e63dc391891ba72437d0f0ca64fbf575c..aef19d81461d91a7b950847f2e2772d5a4b4cfb2 100644
--- a/qim3d/features/_common_features_methods.py
+++ b/qim3d/features/_common_features_methods.py
@@ -1,173 +1,10 @@
 import numpy as np
 import qim3d.processing
 from qim3d.utils._logger import log
-import trimesh
 import qim3d
 from pygel3d import hmesh
 
-
-def volume(obj: np.ndarray|trimesh.Trimesh, 
-           **mesh_kwargs
-           ) -> float:
-    """
-    Compute the volume of a 3D volume or mesh.
-
-    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.
-
-    Returns:
-        volume (float): The volume of the object.
-
-    Example:
-        Compute volume from a mesh:
-        ```python
-        import qim3d
-
-        # Load a mesh from a file
-        mesh = qim3d.io.load_mesh('path/to/mesh.obj')
-
-        # Compute the volume of the mesh
-        vol = qim3d.features.volume(mesh)
-        print('Volume:', vol)
-        ```
-
-        Compute volume from a np.ndarray:
-        ```python
-        import qim3d
-
-        # 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)
-        ```
-
-    """
-    if isinstance(obj, np.ndarray):
-        log.info("Converting volume to mesh.")
-        obj = qim3d.mesh.from_volume(obj, **mesh_kwargs)
-
-    return obj.volume
-
-
-def area(obj: np.ndarray|trimesh.Trimesh, 
-         **mesh_kwargs
-         ) -> float:
-    """
-    Compute the surface area of a 3D volume or mesh.
-
-    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.
-
-    Returns:
-        area (float): The surface area of the object.
-
-    Example:
-        Compute area from a mesh:
-        ```python
-        import qim3d
-
-        # Load a mesh from a file
-        mesh = qim3d.io.load_mesh('path/to/mesh.obj')
-
-        # Compute the surface area of the mesh
-        area = qim3d.features.area(mesh)
-        area = qim3d.features.area(mesh)
-        print(f"Area: {area}")
-        ```
-
-        Compute area from a np.ndarray:
-        ```python
-        import qim3d
-
-        # 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
-        area = qim3d.features.area(synthetic_blob, level=0.5)
-        area = qim3d.features.area(synthetic_blob, level=0.5)
-        print('Area:', area)
-        ```
-    """
-    if isinstance(obj, np.ndarray):
-        log.info("Converting volume to mesh.")
-        obj = qim3d.mesh.from_volume(obj, **mesh_kwargs)
-
-    return obj.area
-
-
-def sphericity(obj: np.ndarray|trimesh.Trimesh, 
-               **mesh_kwargs
-               ) -> float:
-    """
-    Compute the sphericity of a 3D volume or mesh.
-
-    Sphericity is a measure of how spherical an object is. It is defined as the ratio
-    of the surface area of a sphere with the same volume as the object to the object's
-    actual surface area.
-
-    Args:
-        obj (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.
-
-    Returns:
-        sphericity (float): A float value representing the sphericity of the object.
-
-    Example:
-        Compute sphericity from a mesh:
-        ```python
-        import qim3d
-
-        # Load a mesh from a file
-        mesh = qim3d.io.load_mesh('path/to/mesh.obj')
-
-        # Compute the sphericity of the mesh
-        sphericity = qim3d.features.sphericity(mesh)
-        sphericity = qim3d.features.sphericity(mesh)
-        ```
-
-        Compute sphericity from a np.ndarray:
-        ```python
-        import qim3d
-
-        # 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)
-        ```
-
-    !!! 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)
-
-    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
-
-def volume_pygel3d(obj: np.ndarray|hmesh.Manifold) -> float:
+def volume(obj: np.ndarray|hmesh.Manifold) -> float:
     """
     Compute the volume of a 3D mesh using the Pygel3D library.
 
@@ -180,11 +17,11 @@ def volume_pygel3d(obj: np.ndarray|hmesh.Manifold) -> float:
 
     if isinstance(obj, np.ndarray):
         log.info("Converting volume to mesh.")
-        obj = qim3d.mesh.from_volume_pygel3d(obj)
+        obj = qim3d.mesh.from_volume(obj)
 
     return hmesh.volume(obj)
 
-def area_pygel3d(obj: np.ndarray|hmesh.Manifold) -> float:
+def area(obj: np.ndarray|hmesh.Manifold) -> float:
     """
     Compute the surface area of a 3D mesh using the Pygel3D library.
 
@@ -197,7 +34,7 @@ def area_pygel3d(obj: np.ndarray|hmesh.Manifold) -> float:
 
     if isinstance(obj, np.ndarray):
         log.info("Converting volume to mesh.")
-        obj = qim3d.mesh.from_volume_pygel3d(obj)
+        obj = qim3d.mesh.from_volume(obj)
 
     return hmesh.area(obj)
 
@@ -216,8 +53,8 @@ def sphericity_pygel3d(obj: np.ndarray|hmesh.Manifold) -> float:
         log.info("Converting volume to mesh.")
         obj = qim3d.mesh.from_volume_pygel3d(obj)
 
-    volume = volume_pygel3d(obj)
-    area = area_pygel3d(obj)
+    volume = volume(obj)
+    area = area(obj)
 
     if area == 0:
         log.warning("Surface area is zero, sphericity is undefined.")
diff --git a/qim3d/io/_loading.py b/qim3d/io/_loading.py
index df8320c9aefbeb25c9f02a3f176450bd58229747..0ed9ae83294edfe53bc04a0137f118a1e43dab04 100644
--- a/qim3d/io/_loading.py
+++ b/qim3d/io/_loading.py
@@ -28,7 +28,7 @@ from qim3d.utils._misc import get_file_size, sizeof, stringify_path
 from qim3d.utils import Memory
 from qim3d.utils._progress_bar import FileLoadingProgressBar
 import trimesh
-
+from pygel3d import hmesh
 from typing import Optional, Dict
 
 dask.config.set(scheduler="processes") 
@@ -865,15 +865,16 @@ 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 an an X3D/OBJ/OFF/PLY file.
 
     Args:
         filename (str or os.PathLike): The path to the .obj 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
@@ -882,5 +883,6 @@ def load_mesh(filename: str) -> trimesh.Trimesh:
         mesh = qim3d.io.load_mesh("path/to/mesh.obj")
         ```
     """
-    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 0a1fae81cbb4de7357b9b9c57d627ed183d04813..6dd770b59a85f435c5ba6fc695e174b1a92d19ba 100644
--- a/qim3d/io/_saving.py
+++ b/qim3d/io/_saving.py
@@ -36,7 +36,7 @@ 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 sizeof, stringify_path
 
@@ -464,16 +464,41 @@ def save(
     ).save(path, data)
 
 
-def save_mesh(
-        filename: str, 
-        mesh: trimesh.Trimesh
-        ) -> None:
-    """
-    Save a trimesh object to an .obj file.
+# 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 mesh object to an X3D/OBJ/OFF file. The file format is determined by the file extension.
 
     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.
+        mesh (hmesh.Manifold): A hmesh.Manifold object representing the mesh.
 
     Example:
         ```python
@@ -481,7 +506,7 @@ def save_mesh(
 
         vol = qim3d.generate.noise_object(base_shape=(32, 32, 32),
                                   final_shape=(32, 32, 32),
-                                  noise_scale=0.05,
+                                  noise_scale=0.05,mesh.export(filename)
                                   order=1,
                                   gamma=1.0,
                                   max_value=255,
@@ -491,4 +516,4 @@ def save_mesh(
         ```
     """
     # Export the mesh to the specified filename
-    mesh.export(filename)
\ No newline at end of file
+    hmesh.save(filename, mesh)
\ No newline at end of file
diff --git a/qim3d/mesh/__init__.py b/qim3d/mesh/__init__.py
index ffa1fd604fe366973eaab01beeee9dd7c23e9ca6..932dfd9e1f859491652cfde541a71ba6adc28f24 100644
--- a/qim3d/mesh/__init__.py
+++ b/qim3d/mesh/__init__.py
@@ -1 +1 @@
-from ._common_mesh_methods import from_volume, from_volume_pygel3d
+from ._common_mesh_methods import from_volume
diff --git a/qim3d/mesh/_common_mesh_methods.py b/qim3d/mesh/_common_mesh_methods.py
index c6b1266a82bfa8ec180357d6c3bdda21ca05b441..101539ab7262efe28cd92b2660ec05d298c6dcd6 100644
--- a/qim3d/mesh/_common_mesh_methods.py
+++ b/qim3d/mesh/_common_mesh_methods.py
@@ -1,6 +1,5 @@
 import numpy as np
 from skimage import measure, filters
-import trimesh
 from pygel3d import hmesh
 from typing import Tuple, Any
 from qim3d.utils._logger import log
@@ -8,80 +7,7 @@ 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.
-
-    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`.
-
-    Returns:
-        mesh (trimesh.Trimesh): The generated mesh.
-
-    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)
-        ```
-        <iframe src="https://platform.qim.dk/k3d/mesh_visualization.html" width="100%" height="500" frameborder="0"></iframe>
-    """
-    if volume.ndim != 3:
-        raise ValueError("The input volume must be a 3D numpy array.")
-
-    # Compute the threshold level if not provided
-    if level is None:
-        level = filters.threshold_otsu(volume)
-        log.info(f"Computed level using Otsu's method: {level}")
-
-    # Apply padding to the volume
-    if padding is not None:
-        pad_z, pad_y, pad_x = padding
-        padding_value = np.min(volume)
-        volume = np.pad(
-            volume,
-            ((pad_z, pad_z), (pad_y, pad_y), (pad_x, pad_x)),
-            mode="constant",
-            constant_values=padding_value,
-        )
-        log.info(f"Padded volume with {padding} to shape: {volume.shape}")
-
-    # Call skimage.measure.marching_cubes with user-provided kwargs
-    verts, faces, normals, values = measure.marching_cubes(
-        volume, level=level, step_size=step_size, allow_degenerate=allow_degenerate, **kwargs
-    )
-
-    # Create the Trimesh object
-    mesh = trimesh.Trimesh(vertices=verts, faces=faces)
-
-    # Fix face orientation to ensure normals point outwards
-    trimesh.repair.fix_inversion(mesh, multibody=True) 
-
-    return mesh
-
-
-def from_volume_pygel3d(
-        volume: np.ndarray,
-        **Kwargs
+    **Kwargs
 ) -> hmesh.Manifold:
     
     if volume.ndim != 3:
diff --git a/qim3d/viz/_k3d.py b/qim3d/viz/_k3d.py
index 6fa107d4a21cebd432bf49e49f18a52db4830788..2bf4aa54af77ad0cff8fe57c68fc7c3009be7fcd 100644
--- a/qim3d/viz/_k3d.py
+++ b/qim3d/viz/_k3d.py
@@ -12,7 +12,9 @@ import matplotlib.pyplot as plt
 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
 
 def volumetric(
     img: np.ndarray,
@@ -81,7 +83,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,10 +173,9 @@ def volumetric(
     else:
         return plot
 
-
 def mesh(
-    verts: np.ndarray,
-    faces: np.ndarray,
+    mesh,
+    backend: str = "k3d",
     wireframe: bool = True,
     flat_shading: bool = True,
     grid_visible: bool = False,
@@ -183,76 +183,81 @@ def mesh(
     save: bool = False,
     **kwargs,
 ):
-    """
-    Visualizes a 3D mesh using K3D.
-
+    """Visualize a 3D mesh using `k3d` or `pygel3d`.
+    
     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.HMesh): The input mesh object.
+        backend (str, optional): The visualization backend to use. 
+            Choose between `"k3d"` (default) and `"pygel3d"`.
+        wireframe (bool, optional): If True, displays the mesh as a wireframe.
+            Works both with `"k3d"` and `"pygel3d"`. 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 for both `"k3d"` and `"pygel3d"`. 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` visualization. 
+              See full reference: https://k3d-jupyter.org/reference/factory.plot.html
+            
+            - `pygel3d.display` kwargs: Arguments for `pygel3d` visualization, such as:
+                - `smooth` (bool, default=True): Enables smooth shading.
+                - `data` (optional): Allows embedding custom data in the visualization.
+              See full reference: https://www2.compute.dtu.dk/projects/GEL/PyGEL/pygel3d/jupyter_display.html#display
 
     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.
+    """
 
-    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)
-        ```
-        <iframe src="https://platform.qim.dk/k3d/mesh_visualization.html" width="100%" height="500" frameborder="0"></iframe>
-    """
-    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,
-    )
+    if backend not in ["k3d", "pygel3d"]:
+        raise ValueError("Invalid backend. Choose 'k3d' or 'pygel3d'.")
 
-    # Create plot
-    plot = k3d.plot(grid_visible=grid_visible, **kwargs)
-    plot += plt_mesh
+    # Extract vertex positions and face indices
+    face_indices = list(mesh.faces())
+    vertices_array = np.array(mesh.positions())
 
-    if save:
-        # Save html to disk
-        with open(str(save), "w", encoding="utf-8") as fp:
-            fp.write(plot.get_snapshot())
+    # 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)
 
-    if show:
-        plot.display()
-    else:
-        return plot
+    # 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,
+        )
+
+        plot = k3d.plot(grid_visible=grid_visible, **valid_k3d_kwargs)
+        plot += mesh_plot
+
+        if save:
+            with open(str(save), "w", encoding="utf-8") as fp:
+                fp.write(plot.get_snapshot())
+
+        return plot.display() if show else plot
+
+    elif backend == "pygel3d":
+        jd.set_export_mode(True)
+        return jd.display(mesh, **valid_pygel_kwargs)
diff --git a/test_mesh.ipynb b/test_mesh.ipynb
index c32c3fca598dffdf80d630ee86e3dc025252420f..6893bcf0e769a2931e9a729b9d32f589ae463a1a 100644
Binary files a/test_mesh.ipynb and b/test_mesh.ipynb differ