diff --git a/field.py b/field.py
index 1007b46..d36d985 100644
--- a/field.py
+++ b/field.py
@@ -790,20 +790,46 @@ class BinaryFieldNd:
         the volume of the corresponding region is returned.
         Note: it is more efficient to retrieve all volumes at once than querying
         single labels.'''
+        from scipy import ndimage
         if self.labels is None:
             self.label()
         if label is None:
             return np.bincount(self.labels.ravel())
-            # TBD: compare performance with
-            # sizes = ndimage.sum(mask, label_im, range(nb_labels + 1))
-            # http://scipy-lectures.org/advanced/image_processing/auto_examples/plot_find_object.html
+            # return ndimage.sum_labels(self.data,labels=self.labels,index=range(0,self.nlabels+1)) # performs much worse
         else:
             return np.sum(self.labels==label)
+
     
     def volume_domain(self):
         '''Returns volume of entire domain. Should be equal to sum(volume_feature()).'''
         return np.prod(self._dim)
 
+    def area_feature(self,ignore_curved=False,ignore_pathological=False):
+        from scipy import ndimage
+        from time import time
+        # Pad array in periodic directions
+        pw = tuple((1,1) if x else (0,0) for x in self.periodicity)
+        sl_pad = tuple(slice(1,-1) if x else slice(None) for x in self.periodicity)
+        data_ = np.pad(self.data,pw,mode='wrap')
+        #
+        area = np.zeros(self.nlabels+1)
+        stypes = list(range(1,10))
+        if ignore_pathological: del stypes[6:9]
+        if ignore_curved: del stypes[3:6]
+        print(stypes,ignore_curved,ignore_pathological)
+        for stype in stypes:
+            print(stype)
+            svs_,w_ = BinaryFieldNd.surface_voxel_structure(stype)
+            hit_ = np.zeros_like(data_)
+            t = time()
+            ii=0
+            for s_ in svs_: # Loop over permutations
+                ii+=1
+                np.logical_or(hit_,ndimage.binary_hit_or_miss(data_,structure1=s_),out=hit_) # Check if this works inplace
+                print(ii,time()-t,np.sum(hit_))
+            area += w_*np.bincount(self.labels[hit_[sl_pad]],minlength=self.nlabels+1)
+        return area
+
     def feature_labels_by_volume(self,descending=True):
         '''Returns labels of connected regions sorted by volume.'''
         labels = np.argsort(self.volume_feature()[1:])+1
@@ -923,40 +949,42 @@ class BinaryFieldNd:
             raise ValueError('Invalid input. Accepting int,list,tuple,ndarray.')
 
     @staticmethod
-    def surface_voxel_structure(id):
+    def surface_voxel_structure(stype):
         '''Get structure element and all its permutations for surface voxels
            as defined by 
            Mullikin, J. C., & Verbeek, P. W. (1993). Surface area estimation of digitized planes. Bioimaging, 1(1), 6-16.
            Also returns weights, where S7-S9 are defined in 
            Windreich, G., Kiryati, N., & Lohmann, G. (2003). Voxel-based surface area estimation: from theory to practice. Pattern Recognition, 36(11), 2531-2541.
         '''
-        if id==1:
+        if stype==1:
             S = np.array([[[ True, True, True], [ True, True, True], [ True, True, True]],[[ True, True, True], [ True, True, True], [ True, True, True]],[[False,False,False], [False,False,False], [False,False,False]]])
             W = 0.894
-        elif id==2:
+        elif stype==2:
             S = np.array([[[ True, True, True], [ True, True, True], [ True, True, True]],[[ True, True,False], [ True, True,False], [ True, True,False]],[[ True,False,False], [ True,False,False], [ True,False,False]]])
             W = 1.3409
-        elif id==3:
+        elif stype==3:
             S = np.array([[[ True, True, True], [ True, True,False], [ True,False,False]],[[ True, True,False], [ True, True,False], [False,False,False]],[[ True,False,False], [False,False,False], [False,False,False]]])
             W = 1.5879
-        elif id==4:
+        elif stype==4:
             S = np.array([[[ True, True, True], [ True, True, True], [ True, True, True]],[[False, True,False], [False, True,False], [False, True,False]],[[False,False,False], [False,False,False], [False,False,False]]])
             W = 2.0
-        elif id==5:
+        elif stype==5:
             S = np.array([[[ True,False,False], [False,False,False], [False,False,False]],[[ True, True,False], [ True, True,False], [False,False,False]],[[ True,False,False], [False,False,False], [False,False,False]]])
             W = 8./3.
-        elif id==6:
+        elif stype==6:
             S = np.array([[[False,False,False], [ True,False,False], [False,False,False]],[[False,False,False], [ True, True,False], [False,False,False]],[[False,False,False], [ True,False,False], [False,False,False]]])
             W = 10./3.
-        elif id==7:
+        elif stype==7:
             S = np.array([[[False,False,False], [False,False,False], [False,False,False]],[[ True, True, True], [ True, True, True], [ True, True, True]],[[False,False,False], [False,False,False], [False,False,False]]])
             W = 1.79
-        elif id==8:
+        elif stype==8:
             S = np.array([[[False,False,False], [False,False,False], [False,False,False]],[[False,False,False], [ True, True, True], [False,False,False]],[[False,False,False], [False,False,False], [False,False,False]]])
             W = 2.68
-        elif id==9:
+        elif stype==9:
             S = np.array([[[False,False,False], [False,False,False], [False,False,False]],[[False,False,False], [False, True,False], [False,False,False]],[[False,False,False], [False,False,False], [False,False,False]]])
             W = 4.08
+        else:
+            raise ValueError('Invalid structure type.')
         return (np.unique(BinaryFieldNd.rotations48(S).reshape(-1,3,3,3),axis=0),W)
     
     @staticmethod