diff --git a/field.py b/field.py
index 9d02f31..453817e 100644
--- a/field.py
+++ b/field.py
@@ -587,18 +587,56 @@ def gaussian_filter_umean_channel(array,spacing,sigma,truncate=4.0):
     array = ndimage.gaussian_filter1d(array,sigma_img,axis=1,truncate=truncate,mode='mirror')
     return array
 
+class VoxelThreshold:
+    def __init__(self,data,threshold,invert=False):
+        assert isinstance(data,np.ndarray),\
+               "'data' must be a numpy array."
+        self._dim  = data.shape
+        self._ndim = data.ndim
+        if invert:
+            self.data = data<threshold
+        else:
+            self.data = data>=threshold
+
+    @classmethod
+    def from_field(cls,fld3d,threshold,invert=False):
+        return cls(fld3d.data,threshold,invert=invert)
+
+    def fill_holes(self,periodicity=(False,False,False)):
+        '''Fills topological holes in threshold regions.'''
+        assert all([isinstance(x,(bool,int)) for x in periodicity]),\
+                "'periodicity' requires bool values."
+        from scipy import ndimage
+        binarr = ndimage.binary_fill_holes(self.data)
+        for axis in range(self._ndim):
+            if periodicity[axis]:
+                n = binarr.shape[axis]
+                binarr = np.roll(binarr,n//2,axis=axis)
+                binarr = ndimage.binary_fill_holes(binarr)
+                binarr = np.roll(binarr,-n//2,axis=axis)
+        self.data = binarr
+        return
+
+    def probe(self,idx):
+        '''Returns whether or not point at index is inside threshold region or not.'''
+        return self.data[tuple(idx)]
+
+    def volume(self):
+        '''Returns volume of region above threshold.'''
+        return np.sum(self.data)
+
 class ConnectedRegions:
-    def __init__(self,boolarr,periodicity,connect_diagonals=False,bytes_label=32):
-        assert isinstance(boolarr,np.ndarray) and boolarr.dtype==np.dtype('bool'),\
-               "'boolarr' must be a numpy array of dtype('bool')."
+    def __init__(self,binarr,periodicity,connect_diagonals=False,fill_holes=False,bytes_label=32):
+        assert isinstance(binarr,np.ndarray) and binarr.dtype==np.dtype('bool'),\
+               "'binarr' must be a numpy array of dtype('bool')."
         assert all([isinstance(x,(bool,int)) for x in periodicity]),\
                 "'periodicity' requires bool values."
         assert bytes_label in (8,16,32,64),\
                "'bytes_label' must be one of {8,16,32,64}."
-        self._dim  = boolarr.shape
-        self._ndim = boolarr.ndim
+        self._dim  = binarr.shape
+        self._ndim = binarr.ndim
         assert self._ndim in (2,3),\
-               "'boolarr' must be either two or three dimensional."
+               "'binarr' must be either two or three dimensional."
         assert len(periodicity)==self._ndim,\
                "Length of 'periodicity' must match number of dimensions of data."
         from scipy import ndimage
@@ -625,7 +663,7 @@ class ConnectedRegions:
         # this does not take into account periodic wrapping
         dtype_label = np.dtype('uint'+str(bytes_label))
         self.label = np.empty(self._dim,dtype=dtype_label)
-        ndimage.label(boolarr,structure=connectivity,output=self.label)
+        ndimage.label(binarr,structure=connectivity,output=self.label)
         self.count = np.max(self.label)
         # Merge labels if there are periodic overlaps
         map_tgt = np.array(range(0,self.count+1),dtype=dtype_label)
@@ -635,9 +673,9 @@ class ConnectedRegions:
             # Merge the first and last plane and compute connectivity
             sl = self._ndim*[slice(None)]
             sl[axis] = (-1,0)
-            boolarr_ = boolarr[tuple(sl)]
-            label_   = np.empty(boolarr_.shape,dtype=dtype_label)
-            ndimage.label(boolarr_,structure=connectivity,output=label_)
+            binarr_  = binarr[tuple(sl)]
+            label_   = np.empty(binarr_.shape,dtype=dtype_label)
+            ndimage.label(binarr_,structure=connectivity,output=label_)
             for val_ in np.unique(label_):
                 # Get all global labels which are associated to a region 
                 # connected over the boundary
@@ -659,13 +697,18 @@ class ConnectedRegions:
         self.count = np.max(map_tgt)
 
     @classmethod
-    def from_field(cls,fld3d,val,periodicity,connect_diagonals=False,bytes_label=32,invert_threshold=False):
-        if invert_threshold:
-            return cls(fld3d.data<val,periodicity,
-                       connect_diagonals=connect_diagonals,bytes_label=bytes_label)
-        else:
-            return cls(fld3d.data>=val,periodicity,
-                       connect_diagonals=connect_diagonals,bytes_label=bytes_label)
+    def from_field(cls,fld3d,threshold,periodicity,connect_diagonals=False,bytes_label=32,invert=False):
+        voxthr = VoxelThreshold.from_field(fld3d,threshold,invert=invert)
+        return cls.from_voxelthresh(voxthr,periodicity,
+                            connect_diagonals=connect_diagonals,
+                            bytes_label=bytes_label)
+
+    @classmethod
+    def from_voxelthresh(cls,voxthr,periodicity,connect_diagonals=False,bytes_label=32):
+        return cls(voxthr.data,periodicity,
+                       connect_diagonals=connect_diagonals,
+                       bytes_label=bytes_label)
+
 
     def volume(self,label=None):
         '''Returns volume of labeled regions. If 'label' is None all volumes
@@ -702,6 +745,10 @@ class ConnectedRegions:
         self.label = map_tgt[self.label]
         self.count = np.max(map_tgt)
 
+    def probe(self,idx):
+        '''Returns label for given index.'''
+        return self.label[tuple(idx)]
+
     def vtk_contour(self,fld3,val,selection):
         '''Computes contours of a Field3d only within selected structures.'''
         assert isinstance(fld3,Field3d), "'fld3' must be a Field3d instance."