diff --git a/field.py b/field.py
index 5b053e1..1a02636 100644
--- a/field.py
+++ b/field.py
@@ -762,9 +762,15 @@ class Features3d:
         # the volume computation, an arbitrary component of the normal has to be chosen which
         # defaults to the z-component and is set by 'cellvol_normal_component'.
         t = time()
-        self._cell_areas   = contour_.compute_cell_sizes(length=False,area=True,volume=False)["Area"]
-        self._cell_volumes = contour_.cell_normals[:,cellvol_normal_component]*\
-                             contour_.cell_centers().points[:,cellvol_normal_component]*self._cell_areas
+        faces  = contour_.faces.reshape(contour_.n_faces,4)[:,:]
+        points = contour_.points
+        X = points[faces[:,1],:]
+        Y = points[faces[:,2],:]
+        Z = points[faces[:,3],:]
+        cn = np.cross(X-Z,Y-X)
+        cc = (X+Y+Z)/3
+        area = 0.5*np.sqrt(np.square(cn).sum(axis=1))
+        vol  = 0.5*cn[:,cellvol_normal_component]*cc[:,cellvol_normal_component]
         print('CELLPROPERTIES:',time()-t)
         # Now the label is known per cell. We only need to find all cells with the same label
         # and group them. Internally we will store the points in one array and the cells in 
@@ -774,9 +780,9 @@ class Features3d:
         ind    = np.argsort(cell_val)
         self._offset   = offset
         self._faces    = contour_.faces.reshape(contour_.n_faces,4)[ind,:]
-        self._points   = contour_.points
-        self._cell_areas   = self._cell_areas[ind]
-        self._cell_volumes = self._cell_volumes[ind]
+        self._points   = points
+        self._cell_areas   = area[ind]
+        self._cell_volumes = vol[ind]
         print('FINALIZING:',time()-t)
         return