diff --git a/qim3d/features/_common_features_methods.py b/qim3d/features/_common_features_methods.py
index b448b1d300b4985231262b0428eee158611cd453..62a7611e63dc391891ba72437d0f0ca64fbf575c 100644
--- a/qim3d/features/_common_features_methods.py
+++ b/qim3d/features/_common_features_methods.py
@@ -3,6 +3,7 @@ import qim3d.processing
 from qim3d.utils._logger import log
 import trimesh
 import qim3d
+from pygel3d import hmesh
 
 
 def volume(obj: np.ndarray|trimesh.Trimesh, 
@@ -89,15 +90,14 @@ def area(obj: np.ndarray|trimesh.Trimesh,
         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, 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)
-        obj = qim3d.mesh.from_volume(obj, **mesh_kwargs)
 
     return obj.area
 
@@ -153,7 +153,6 @@ def sphericity(obj: np.ndarray|trimesh.Trimesh,
     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)
 
     volume = qim3d.features.volume(obj)
     area = qim3d.features.area(obj)
@@ -167,3 +166,63 @@ def sphericity(obj: np.ndarray|trimesh.Trimesh,
     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:
+    """
+    Compute the volume of a 3D mesh using the Pygel3D library.
+
+    Args:
+        obj: Either a np.ndarray volume or a mesh object of type hmesh.Manifold.
+
+    Returns:
+        volume (float): The volume of the object. 
+    """
+
+    if isinstance(obj, np.ndarray):
+        log.info("Converting volume to mesh.")
+        obj = qim3d.mesh.from_volume_pygel3d(obj)
+
+    return hmesh.volume(obj)
+
+def area_pygel3d(obj: np.ndarray|hmesh.Manifold) -> float:
+    """
+    Compute the surface area of a 3D mesh using the Pygel3D library.
+
+    Args:
+        obj: Either a np.ndarray volume or a mesh object of type hmesh.Manifold.
+
+    Returns:
+        area (float): The surface area of the object. 
+    """
+
+    if isinstance(obj, np.ndarray):
+        log.info("Converting volume to mesh.")
+        obj = qim3d.mesh.from_volume_pygel3d(obj)
+
+    return hmesh.area(obj)
+
+def sphericity_pygel3d(obj: np.ndarray|hmesh.Manifold) -> float:
+    """
+    Compute the sphericity of a 3D mesh using the Pygel3D library.
+
+    Args:
+        obj: Either a np.ndarray volume or a mesh object of type hmesh.Manifold.
+
+    Returns:
+        sphericity (float): The sphericity of the object. 
+    """
+
+    if isinstance(obj, np.ndarray):
+        log.info("Converting volume to mesh.")
+        obj = qim3d.mesh.from_volume_pygel3d(obj)
+
+    volume = volume_pygel3d(obj)
+    area = area_pygel3d(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
\ No newline at end of file
diff --git a/qim3d/mesh/__init__.py b/qim3d/mesh/__init__.py
index 932dfd9e1f859491652cfde541a71ba6adc28f24..ffa1fd604fe366973eaab01beeee9dd7c23e9ca6 100644
--- a/qim3d/mesh/__init__.py
+++ b/qim3d/mesh/__init__.py
@@ -1 +1 @@
-from ._common_mesh_methods import from_volume
+from ._common_mesh_methods import from_volume, from_volume_pygel3d
diff --git a/qim3d/mesh/_common_mesh_methods.py b/qim3d/mesh/_common_mesh_methods.py
index a7e358ea7eed9bade283bd411f33da6f64ffa3a0..c6b1266a82bfa8ec180357d6c3bdda21ca05b441 100644
--- a/qim3d/mesh/_common_mesh_methods.py
+++ b/qim3d/mesh/_common_mesh_methods.py
@@ -1,6 +1,7 @@
 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
 
@@ -76,3 +77,15 @@ def from_volume(
     trimesh.repair.fix_inversion(mesh, multibody=True) 
 
     return mesh
+
+
+def from_volume_pygel3d(
+        volume: np.ndarray,
+        **Kwargs
+) -> hmesh.Manifold:
+    
+    if volume.ndim != 3:
+        raise ValueError("The input volume must be a 3D numpy array.")
+    
+    m = hmesh.volumetric_isocontour(volume)
+    return m
\ No newline at end of file
diff --git a/test_mesh.ipynb b/test_mesh.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..c32c3fca598dffdf80d630ee86e3dc025252420f
Binary files /dev/null and b/test_mesh.ipynb differ