diff --git a/field.py b/field.py
index 7de8769..6decfe7 100644
--- a/field.py
+++ b/field.py
@@ -648,109 +648,82 @@ def gaussian_filter_umean_channel(array,spacing,sigma,truncate=4.0):
 
 
 class Features3d:
-    def __init__(self,input,threshold,origin,spacing,periodicity,invert=False,has_ghost=False):
+    def __init__(self,input,threshold,origin,spacing,periodicity,
+                        invert=False,has_ghost=False,keep_input=False,
+                        contour_method='flying_edges',cellvol_normal_component=2,
+                        report=False):
         assert len(origin)==3,  "'origin' must be of length 3"
         assert len(spacing)==3, "'spacing' must be of length 3"
         assert len(periodicity)==3, "'periodicity' must be of length 3"
         assert isinstance(input,np.ndarray), "'input' must be numpy.ndarray."
         # Assign basic properties to class variables
-        self.origin  = tuple(float(x) for x in origin)
-        self.spacing = tuple(float(x) for x in spacing)
+        self.origin      = np.array(origin,dtype=np.float)
+        self.spacing     = np.array(spacing,dtype=np.float)
+        self.dimensions  = input.shape
         self.periodicity = tuple(bool(x) for x in periodicity)
-        self._invert = invert
-        self._sldata = tuple(slice(0,-1) if x else None for x in periodicity)
         # If regions are supposed to be inverted, i.e. the interior consists of values 
         # smaller than the threshold instead of larger, change the sign of the array.
         sign_invert = -1 if invert else +1
         self._threshold = sign_invert*threshold
-        # '_data' is the array which is to be triangulated. Here one trailing ghost cell is 
+        # '_input' is the array which is to be triangulated. Here one trailing ghost cell is 
         # required to get closed surfaces. The data will always be copied since it probably 
         # needs to be copied for vtk anyway (needs FORTRAN contiguous array) and this way
         # there will never be side effects on the input data.
         if has_ghost:
-            self._data = sign_invert*input
+            self._input = sign_invert*input
         else:
             pw = tuple((0,1) if x else (0,0) for x in periodicity)
-            self._data = np.pad(sign_invert*input,pw,mode='wrap')
-
-
+            self._input = np.pad(sign_invert*input,pw,mode='wrap')
+        # Triangulate
+        self.triangulate(contour_method=contour_method,
+                         cellvol_normal_component=cellvol_normal_component,
+                         report=report)
+        # Drop input if requested: this may free memory in case it was copied/temporary.
+        if not keep_input: self._input = None
+        return
 
     @classmethod
-    def from_field(cls,fld,threshold,periodicity,invert=False,has_ghost=False):
-        return cls(fld.data,threshold,fld.origin,fld.spacing,periodicity,invert=invert,has_ghost=has_ghost)
+    def from_field(cls,fld,threshold,periodicity,invert=False,has_ghost=False,
+                            keep_input=False,contour_method='flying_edges',
+                            cellvol_normal_component=2,report=False):
+        return cls(fld.data,threshold,fld.origin,fld.spacing,periodicity,
+                        invert=invert,has_ghost=has_ghost,keep_input=keep_input,
+                        contour_method=contour_method,
+                        cellvol_normal_component=cellvol_normal_component,
+                        report=report)
 
-    @property
-    def threshold(self): return self._threshold
-
-    @property
-    def points(self): return self._points
-
-    @property
-    def faces(self): return self._faces
-
-    @property
-    def nfeatures(self): return self._nfeatures
-
-    def fill_holes(self):
-        # Check if volume negative -> cell normal direction
-
-        # self.binary.fill_holes()
-        # li = ndimage.binary_erosion(self.binary._data)
-        # self._data[li] = 2*self._threshold # arbitrary, only needs to be above threshold
-        # if self._faces is not None:
-        #     self.triangulate(contour_method=self.__TRI_CONTMETH,
-        #                      cellvol_normal_component=self.__TRI_NORMCOMP)
-        return
-
-    def reduce_noise(self,threshold=1e-5,is_relative=True):
-        '''Removes all objects with smaller (binary-based) volume than a threshold.'''
-        # if is_relative:
-        #     threshold = threshold*self.binary.volume_domain()
-        # vol = self.binary.volumes()
-        # li  = vol<threshold
-        # mask = self.binary.discard_features(li,return_mask=True)
-        # self._data[mask] = np.nan#self._threshold
-        # if self._faces is not None:
-        #     self.triangulate(contour_method=self.__TRI_CONTMETH,
-        #                      cellvol_normal_component=self.__TRI_NORMCOMP)
-        return
-
-    def triangulate(self,contour_method='flying_edges',cellvol_normal_component=2):
+    def triangulate(self,contour_method='flying_edges',cellvol_normal_component=2,report=False):
         import pyvista as pv
         import vtk
         from scipy import ndimage, spatial
-        from time import time
-        # Save arguments in case we need to recreate triangulation later on
-        self.__TRI_CONTMETH = contour_method
-        self.__TRI_NORMCOMP = cellvol_normal_component
-        # Wrap data with pyvista/VTK
+        # Check if '_input' is available: might have been dropped after initialization
+        assert self._input is not None, "'_input' not available. Initialize object with keep_input=True flag."
+        # Wrap data for VTK using pyvista
         datavtk = pv.UniformGrid()
-        datavtk.dimensions = self._data.shape
+        datavtk.dimensions = self._input.shape
         datavtk.origin     = self.origin
         datavtk.spacing    = self.spacing
-        datavtk.point_arrays['data'] = self._data.ravel('F')
+        datavtk.point_arrays['data'] = self._input.ravel('F')
         # Compute full contour: ensure we only have triangles to efficiently extract points
         # later on.
-        t = time()
+        if report: print('[Features3d.triangulate] computing isocontour using {}...'.format(contour_method))
         contour_ = datavtk.contour([self._threshold],method=contour_method,compute_scalars=False)
         assert contour_.is_all_triangles(), "Contouring produced non-triangle cells."
-        print('CONTOUR:',time()-t)
-        # Compute the connectivity of the triangulated surface: first we run a normal 
+        # Compute the connectivity of the triangulated surface: first we run an ordinary 
         # connectivity filter neglecting periodic wrapping.
-        t = time()
+        if report: print('[Features3d.triangulate] computing connectivity...')
         alg = vtk.vtkPolyDataConnectivityFilter() # ~twice as fast as default vtkConnectivityFilter
         alg.SetInputData(contour_)
         alg.SetExtractionModeToAllRegions()
         alg.SetColorRegions(True)
         alg.Update()
         contour_ = pv.filters._get_output(alg)
-        print('connectivity computed in',time()-t,'seconds')
         # Now determine the boundary points, i.e. points on the boundary which have a corresponding
         # vertex on the other side of the domain
-        t = time()
+        if report: print('[Features3d.triangulate] correcting connectivity for periodic boundaries...')
         points = contour_.points
         bd_points_xyz = np.zeros((0,3),dtype=points.dtype)
-        bd_points_idx = np.zeros((0,),dtype=np.int64)
+        bd_points_idx = np.zeros((0,),dtype=np.int)
         for axis in range(3):
             if not self.periodicity[axis]: continue
             # Compute position of boundary, period of domain and set a threshold for "on boundary" condition
@@ -766,18 +739,14 @@ class Features3d:
             li = np.abs(points[:,axis]-bound_pos)<bound_dist
             bd_points_idx = np.append(bd_points_idx,np.nonzero(li)[0],axis=0)
             bd_points_xyz = np.append(bd_points_xyz,points[li,:]-period,axis=0)
-        print('Points selected in',time()-t,'seconds')
         # Construct a KD Tree for efficient neighborhood search
-        t = time()
         kd = spatial.KDTree(bd_points_xyz,leafsize=10,compact_nodes=True,copy_data=False,balanced_tree=True)
-        print('kd tree built in',time()-t,'seconds')
         # Construct a map to get new labels for regions. We search for pairs of points with
         # a very small distance inbetween. Then we check if their labels differ and choose 
         # the lower label as the new joint one.
-        t = time()
         point_labels = contour_.point_arrays['RegionId']
-        nfeatures    = np.max(point_labels)
-        map_         = np.array(range(nfeatures+1),dtype=point_labels.dtype)
+        nfeatures    = np.max(point_labels)+1
+        map_         = np.array(range(nfeatures),dtype=point_labels.dtype)
         overlap_dist = 1e-4*np.sqrt(np.square(datavtk.spacing).sum())
         for (ii,jj) in kd.query_pairs(r=overlap_dist):
             label_ii = point_labels[bd_points_idx[ii]]
@@ -791,29 +760,17 @@ class Features3d:
         # 
         map_         = np.unique(map_,return_inverse=True)[1]
         point_labels = map_[point_labels]
-        nfeatures    = np.max(map_)
-        print('mapped overlaps in',time()-t,'seconds')
-
-        # pl = pv.Plotter()
-        # pl.add_mesh(contour_,scalars=point_labels,opacity=1.0,specular=1.0,interpolate_before_map=True)
-        # pl.show()
-
+        nfeatures    = np.max(map_)+1 # starts with zero
         # Labels are now stored as point data. To efficiently convert it to cell data, the first 
         # point of each cell determines the value for this cell.
-        t = time()
+        if report: print('[Features3d.triangulate] identifying cell based labels...')
         ncells = contour_.n_faces
         faces  = contour_.faces.reshape(ncells,4)[:,:]
-        # cell_labels = np.zeros((ncells,))
-        # for ii,face in enumerate(faces):
-        #     cell_labels[ii] = point_labels[face[1]]
         cell_labels = point_labels[faces[:,1]]
-        print('cell interpolation in',time()-t,'seconds')
-        # While we have the full contour in memory, compute the area and volume per cell. For
-        # 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()
-        # faces  = contour_.faces.reshape(contour_.n_faces,4)[:,:]
-        # points = contour_.points
+        # Compute the volume and area per cell. For 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'.
+        if report: print('[Features3d.triangulate] calculating area and volume per cell...')
         X = points[faces[:,1],:]
         Y = points[faces[:,2],:]
         Z = points[faces[:,3],:]
@@ -821,20 +778,77 @@ class Features3d:
         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 
-        # nlabel tuples of ndarrays.
-        t = time()
-        offset = np.cumsum(np.bincount(cell_labels))[1:-1]
+        # and group them. Therefore we sort the array holding the faces by feature label
+        # and save the offset in another array to easily extract them whenever necessary.
+        if report: print('[Features3d.triangulate] finalizing internal data...')
+        offset = np.zeros((nfeatures+1,),dtype=np.int64)
+        offset[1:] = np.cumsum(np.bincount(cell_labels))
+        print(nfeatures,np.bincount(cell_labels),offset)
         ind    = np.argsort(cell_labels)
         self._offset   = offset
-        self._faces    = faces[ind,1:]
+        self._faces    = faces[ind,:]
         self._points   = points
         self._cell_areas   = area[ind]
         self._cell_volumes = vol[ind]
         self._nfeatures = nfeatures
-        print('FINALIZING:',time()-t)
+        return
+
+    @property
+    def threshold(self): return self._threshold
+
+    @property
+    def points(self): return self._points
+
+    @property
+    def faces(self): return self._faces
+
+    @property
+    def nfeatures(self): return self._nfeatures
+
+    def fill_holes(self,report=False):
+        '''Remove triangulation which is fully enclosed by another. These regions 
+           will have a negative volume due to the direction of their normal vector.'''
+        self.discard_features(np.flatnonzero(self.volumes()<0),report=report)
+        return
+
+    def reduce_noise(self,threshold=1e-5,is_relative=True,report=False):
+        '''Discards all objects with smaller volume than a threshold.'''
+        if is_relative:
+            vol_domain = np.prod(self.spacing*self.dimensions)
+            threshold = threshold*vol_domain
+        self.discard_features(np.flatnonzero(self.volumes()<threshold),report=report)
+        return
+
+    def discard_features(self,labels,clean_points=False,report=False):
+        '''Removes features from triangulated data.'''
+        from collections.abc import Iterable
+        from time import time
+        if not isinstance(labels,Iterable): labels = [labels]
+        # Get index ranges which are to be deleted and also create an array
+        # which determines the size of the block to be deleted
+        idx = []
+        gapsize = np.zeros((self.nfeatures,),dtype=np.int)
+        for label in labels:
+            rng_ = np.arange(self._offset[label],self._offset[label+1])
+            idx.append(rng_)
+            gapsize[label] = len(rng_)
+        idx = np.concatenate(idx,axis=0)
+        # Save former number of faces for report
+        nfaces = self._faces.shape[0]
+        # Delete indexed elements from arrays
+        self._faces        = np.delete(self._faces,idx,axis=0)
+        self._cell_areas   = np.delete(self._cell_areas,idx,axis=0)
+        self._cell_volumes = np.delete(self._cell_volumes,idx,axis=0)
+        # Correct offset
+        self._offset[1:] = self._offset[1:]-np.cumsum(gapsize)
+        self._offset = self._offset[:-len(labels)]
+        if report:
+            print('[Features3d.discard_features]',end=' ')
+            print('discarded {:} features with {:} faces.'.format(
+                            len(labels),nfaces-self._faces.shape[0]))
+        if clean_points:
+            self.clean_points(report=report)
         return
 
     @property
@@ -844,10 +858,10 @@ class Features3d:
     def points(self): return self._points
 
     @property
-    def cell_areas(self): return np.split(self._cell_areas,self._offset)
+    def cell_areas(self): return np.split(self._cell_areas,self._offset[1:-1])
 
     @property
-    def cell_volumes(self): return np.split(self._cell_volumes,self._offset)
+    def cell_volumes(self): return np.split(self._cell_volumes,self._offset[1:-1])
 
     def area(self,feature): 
         '''Returns the surface area of feature. If feature is None, total surface 
@@ -858,36 +872,22 @@ class Features3d:
             return np.sum(self.cell_areas[feature-1])
 
     def areas(self): 
-        '''Returns a tuple with surface areas of all features.'''
-        return tuple(np.sum(x) for x in self.cell_areas)
+        '''Returns an array with surface areas of all features.'''
+        return np.add.reduceat(self._cell_areas,self._offset[:-1])
 
-    def volume(self,feature,method='triangulation'):
+    def volume(self,feature):
         '''Returns volume enclosed by feature. If feature isNone, total volume of 
            all features is returned.'''
-        if method=='triangulation':
-            if self._faces is None: 
-                self.triangulate(contour_method=self.__TRI_CONTMETH,
-                                 cellvol_normal_component=self.__TRI_NORMCOMP)
-            if feature is None: 
-                return np.sum(self._cell_volumes)
-            else: 
-                return np.sum(self.cell_volumes[feature-1])
-        elif method=='binary':
-            return self.binary.volume(feature)
-        else:
-            raise ValueError("Invalid method. Available methods: 'triangulation','binary'")
+        if feature is None: 
+            return np.sum(self._cell_volumes)
+        else: 
+            return np.sum(self.cell_volumes[feature-1])
 
-    def volumes(self,method='triangulation'):
-        if method=='triangulation':
-            self.triangulate(contour_method=self.__TRI_CONTMETH,
-                             cellvol_normal_component=self.__TRI_NORMCOMP)
-            return tuple(np.sum(x) for x in self.cell_volumes)
-        elif method=='binary':
-            return self.binary.volumes()
-        else:
-            raise ValueError("Invalid method. Available methods: 'triangulation','binary'")
+    def volumes(self):
+        '''Returns an array with volumes of all features.'''
+        return np.add.reduceat(self._cell_volumes,self._offset[:-1])
 
-    def which_feature(self,coords,method='binary'):
+    def find_feature(self,coords):
         coords = np.array(coords)
         assert coords.ndim==2 and coords.shape[1]==3, "'coords' need to be provided as 3xN array."
         # if method=='binary':
@@ -895,21 +895,38 @@ class Features3d:
         # elif method=='triangulation':
         return
 
-    def to_vtk(self,label,reduce_points=False):
-        if self._faces is None:
-            raise RuntimeError('Features have not been triangulated yet. Call triangulate() method first.')
+    def clean_points(self,report=False):
+        nfaces_ = self._faces.shape[0]
+        ind,inv = np.unique(self._faces[:,1:],return_inverse=True)
+        if report:
+            print('[Features3d.clean_points]',end=' ') 
+            print('removed {:} orphan points.'.format(self._points.shape[0]-len(ind)))
+        self._faces[:,1:4] = inv.reshape(nfaces_,3)
+        self._points = self._points[ind,:]
+
+    def to_vtk(self,labels,reduce_points=False):
         import pyvista as pv
-        idx = label-1
-        if reduce_points:
-            nfaces_ = self._faces[idx].shape[0]
-            ind,inv = np.unique(self._faces[idx][:,1:4],return_inverse=True)
-            faces_  = np.zeros((nfaces_,4),dtype=self._faces[idx].dtype)
-            faces_[:,0]   = 3
-            faces_[:,1:4] = inv.reshape(nfaces_,3)
-            points_ = self._points[ind,:]
-            return pv.PolyData(points_,faces_)
-        else:
-            return pv.PolyData(self._points,self._faces[idx])
+        from collections.abc import Iterable
+        if labels is None:
+            return pv.PolyData(self._points,self._faces)
+        if not isinstance(labels,Iterable): 
+            labels = [labels]
+        points = self._points
+        faces  = []
+        for label in labels:
+            sl_ = slice(self._offset[label],self._offset[label+1])
+            faces.append(self._faces[sl_,:])
+        faces = np.concatenate(faces,axis=0)
+        # if reduce_points:
+        #     nfaces_ = self._faces[idx].shape[0]
+        #     ind,inv = np.unique(self._faces[:,1:4],return_inverse=True)
+        #     faces_  = np.zeros((nfaces_,4),dtype=self._faces.dtype)
+        #     faces_[:,0]   = 3
+        #     faces_[:,1:4] = inv.reshape(nfaces_,3)
+        #     points_ = self._points[ind,:]
+        #     return pv.PolyData(points_,faces_)
+        # else:
+        return pv.PolyData(points,faces)
 
 class BinaryFieldNd:
     def __init__(self,input,periodicity,connect_diagonals=False,deep=False,has_ghost=False):