diff --git a/field.py b/field.py
index 6755f3d..6ae11a1 100644
--- a/field.py
+++ b/field.py
@@ -578,10 +578,10 @@ class Field3d:
             mesh.point_arrays['data'] = self.data.ravel(order='F') 
         return mesh
 
-    def vtk_contour(self,val,deep=False):
+    def vtk_contour(self,val,deep=False,method='contour'):
         if not isinstance(val,(tuple,list)):
             val = [val]
-        return self.to_vtk(deep=deep).contour(val)
+        return self.to_vtk(deep=deep).contour(val,method=method)
 
     def vtk_slice(self,normal,origin,deep=False):
         assert (normal in ('x','y','z') or (isinstance(normal,(tuple,list)) 
@@ -795,18 +795,24 @@ class BinaryFieldNd:
                 self.wrap[axis] = self.wrap[axis][idx_]
         self.data = self.labels>0
 
-    def isolate_feature(self,selection):
+    def isolate_feature(self,selection,array=None):
         from scipy import ndimage
         if self.labels is None:
             self.label()
         selection = self._select_feature(selection)
-        output = []
+        output1 = []
+        has_array = array is not None
+        if has_array:
+            assert np.all(array.shape==self._dim)
+            output2 = []
         for lab_ in selection:
             # Extract feature of interest
             if lab_==0:
                 data_ = np.logical_not(self.data)
+                if has_array: data2_ = array
             else:
                 data_ = (self.labels[self._feat_slice[lab_-1]]==lab_)
+                if has_array: data2_ = array[self._feat_slice[lab_-1]]
             # If feature is wrapped periodically, duplicate it and extract
             # largest one
             iswrapped = False
@@ -822,20 +828,42 @@ class BinaryFieldNd:
                 il_    = np.argmax(vol_[1:])+1
                 sl_    = ndimage.find_objects(l_==il_)[0]
                 data_  = data_[sl_]
+                if has_array:
+                    data2_ = np.tile(data2_,rep_)[sl_]                    
             # Add to output
-            output.append(data_)
-        return output
+            output1.append(data_)
+            if has_array: output2.append(data2_)
+        if has_array:
+            return (output1,output2)
+        else:
+            return output1
 
-    def triangulate_feature(self,selection,origin=(0,0,0),spacing=(1,1,1)):
+    def triangulate_feature(self,selection,origin=(0,0,0),spacing=(1,1,1),array=None):
         assert self._ndim==3, "Triangulation requires 3D data."
+        from scipy import ndimage
         if self.labels is None:
             self.label()
         selection = self._select_feature(selection)
         output = []
         pw = tuple((1,1) for ii in range(self._ndim))
+        has_array = array is not None
+        if has_array:
+            assert np.all(array.shape==self._dim)
+            output2 = []
         for lab_ in selection:
-            data_ = np.pad(self.isolate_feature(lab_)[0],pw,mode='constant')
-            output.append(Field3d(data_,origin,spacing).vtk_contour(0.5))
+            if has_array:
+                data_,scal_ = self.isolate_feature(lab_,array=array)
+                data_,scal_ = data_[0],scal_[0]
+                # Fill interior in case we filled holes which is not in scal_
+                data_ = ndimage.binary_erosion(data_,structure=self.structure,iterations=2)
+                scal_[data_] = 1.0
+                scal_ = np.pad(scal_,pw,mode='reflect',reflect_type='odd')
+                pd = Field3d(scal_,origin,spacing).vtk_contour(0.0)
+            else:
+                data_ = self.isolate_feature(lab_)[0]
+                data_ = np.pad(data_.astype(float),pw,mode='constant',constant_values=-1.0)
+                pd = Field3d(data_,origin,spacing).vtk_contour(0.5).smooth(1000)
+            output.append(pd)
         return output
 
     def _select_feature(self,selection):