diff --git a/Common/DataModel/CMakeLists.txt b/Common/DataModel/CMakeLists.txt
index 37ac41beb36cca8a6d70a3068daa4fe2fd3e2d6f..a092506071e63631362d5b7ed76befbe86cb5d03 100644
--- a/Common/DataModel/CMakeLists.txt
+++ b/Common/DataModel/CMakeLists.txt
@@ -9,6 +9,7 @@ set(Module_SRCS
   vtkAnnotation.cxx
   vtkAnnotationLayers.cxx
   vtkArrayData.cxx
+  vtkArrayListTemplate.txx
   vtkAttributesErrorMetric.cxx
   vtkBiQuadraticQuad.cxx
   vtkBiQuadraticQuadraticHexahedron.cxx
@@ -213,6 +214,7 @@ set(Module_SRCS
   )
 
 set(${vtk-module}_HDRS
+  vtkArrayListTemplate.h
   vtkCellType.h
   vtkMappedUnstructuredGrid.h
   vtkMappedUnstructuredGridCellIterator.h
@@ -260,6 +262,7 @@ set_source_files_properties(
 
 set_source_files_properties(
   vtkAMRBox
+  vtkArrayListTemplate.txx
   vtkAtom
   vtkBond
   vtkBoundingBox
diff --git a/Filters/Points/vtkArrayListTemplate.h b/Common/DataModel/vtkArrayListTemplate.h
similarity index 67%
rename from Filters/Points/vtkArrayListTemplate.h
rename to Common/DataModel/vtkArrayListTemplate.h
index 00e494cd90edd2f9353bd805fe0b129b04c07426..b926a2c1967af940cc5ef411b2938e0087321bd4 100644
--- a/Filters/Points/vtkArrayListTemplate.h
+++ b/Common/DataModel/vtkArrayListTemplate.h
@@ -26,10 +26,10 @@
 // vtkDataSetAttributes::CopyInterpolate() or InterpolateAllocate() which
 // performs the initial magic of constructing input and output arrays. Then
 // the input attributes, and output attributes, are passed to initialize the
-// internal structures. Internally, pairs of typed arrays are created; the
-// operations (e.g., interpolate) occur on these typed arrays using a
-// typeless, virtual dispatch base class.
-
+// internal structures. Essentially these internal structures are pairs of
+// arrays of the same type, which can be efficently accessed and
+// assigned. The operations on these array pairs (e.g., interpolation) occur
+// using a typeless, virtual dispatch base class.
 
 // .SECTION See Also
 // vtkFieldData vtkDataSetAttributes vtkPointData vtkCellData
@@ -46,8 +46,10 @@ struct BaseArrayPair
 {
   vtkIdType Num;
   int NumComp;
+  vtkDataArray *OutputArray;
 
-  BaseArrayPair(vtkIdType num, int numComp) : Num(num), NumComp(numComp)
+  BaseArrayPair(vtkIdType num, int numComp, vtkDataArray *outArray) :
+    Num(num), NumComp(numComp), OutputArray(outArray)
     {
     }
   virtual ~BaseArrayPair()
@@ -56,8 +58,11 @@ struct BaseArrayPair
 
   virtual void Copy(vtkIdType inId, vtkIdType outId) = 0;
   virtual void Interpolate(int numWeights, const vtkIdType *ids,
-                           const double *weights, vtkIdType outPtId) = 0;
+                           const double *weights, vtkIdType outId) = 0;
+  virtual void InterpolateEdge(vtkIdType v0, vtkIdType v1,
+                               double t, vtkIdType outId) = 0;
   virtual void AssignNullValue(vtkIdType outId) = 0;
+  virtual void Realloc(vtkIdType sze) = 0;
 };
 
 // Type specific interpolation on a matched pair of data arrays
@@ -68,8 +73,8 @@ struct ArrayPair : public BaseArrayPair
   T *Output;
   T  NullValue;
 
-  ArrayPair(T *in, T *out, vtkIdType num, int numComp, T null) :
-    BaseArrayPair(num,numComp), Input(in), Output(out), NullValue(null)
+  ArrayPair(T *in, T *out, vtkIdType num, int numComp, vtkDataArray *outArray, T null) :
+    BaseArrayPair(num,numComp,outArray), Input(in), Output(out), NullValue(null)
     {
     }
   virtual ~ArrayPair()  //calm down some finicky compilers
@@ -85,7 +90,7 @@ struct ArrayPair : public BaseArrayPair
     }
 
   virtual void Interpolate(int numWeights, const vtkIdType *ids,
-                           const double *weights, vtkIdType outPtId)
+                           const double *weights, vtkIdType outId)
     {
     for (int j=0; j < this->NumComp; ++j)
       {
@@ -94,7 +99,19 @@ struct ArrayPair : public BaseArrayPair
         {
         v += weights[i] * static_cast<double>(this->Input[ids[i]*this->NumComp+j]);
         }
-      this->Output[outPtId*this->NumComp+j] = static_cast<T>(v);
+      this->Output[outId*this->NumComp+j] = static_cast<T>(v);
+      }
+    }
+
+  virtual void InterpolateEdge(vtkIdType v0, vtkIdType v1, double t, vtkIdType outId)
+    {
+    double v;
+    vtkIdType numComp=this->NumComp;
+    for (int j=0; j < numComp; ++j)
+      {
+      v = this->Input[v0*numComp+j] +
+        t * (this->Input[v1*numComp+j] - this->Input[v0*numComp+j]);
+      this->Output[outId*numComp+j] = static_cast<T>(v);
       }
     }
 
@@ -106,6 +123,12 @@ struct ArrayPair : public BaseArrayPair
       }
     }
 
+  virtual void Realloc(vtkIdType sze)
+    {
+      this->OutputArray->WriteVoidPointer(0,sze*this->NumComp);
+      this->Output = static_cast<T*>(this->OutputArray->GetVoidPointer(0));
+    }
+
 };
 
 // Forward declarations. This makes working with vtkTemplateMacro easier.
@@ -137,16 +160,26 @@ struct ArrayList
     }
 
   // Loop over the arrays and have them interpolate themselves
-  void Interpolate(int numWeights, const vtkIdType *ids, const double *weights, vtkIdType outPtId)
+  void Interpolate(int numWeights, const vtkIdType *ids, const double *weights, vtkIdType outId)
     {
       for (std::vector<BaseArrayPair*>::iterator it = Arrays.begin();
            it != Arrays.end(); ++it)
         {
-        (*it)->Interpolate(numWeights, ids, weights, outPtId);
+        (*it)->Interpolate(numWeights, ids, weights, outId);
         }
     }
 
-  // Loop over the arrays and have them interpolate themselves
+  // Loop over the arrays perform edge interpolation
+  void InterpolateEdge(vtkIdType v0, vtkIdType v1, double t, vtkIdType outId)
+    {
+      for (std::vector<BaseArrayPair*>::iterator it = Arrays.begin();
+           it != Arrays.end(); ++it)
+        {
+        (*it)->InterpolateEdge(v0, v1, t, outId);
+        }
+    }
+
+  // Loop over the arrays and assign the null value
   void AssignNullValue(vtkIdType outId)
     {
       for (std::vector<BaseArrayPair*>::iterator it = Arrays.begin();
@@ -156,6 +189,16 @@ struct ArrayList
         }
     }
 
+  // Extend (realloc) the arrays
+  void Realloc(vtkIdType sze)
+    {
+      for (std::vector<BaseArrayPair*>::iterator it = Arrays.begin();
+           it != Arrays.end(); ++it)
+        {
+        (*it)->Realloc(sze);
+        }
+    }
+
   // Only you can prevent memory leaks!
   ~ArrayList()
     {
diff --git a/Filters/Points/vtkArrayListTemplate.txx b/Common/DataModel/vtkArrayListTemplate.txx
similarity index 92%
rename from Filters/Points/vtkArrayListTemplate.txx
rename to Common/DataModel/vtkArrayListTemplate.txx
index c83f01b9f122b154115e3c5a3ad98be32cb6757b..a47ebe8d93c1465923d55e586ab2264f22a2ce7e 100644
--- a/Filters/Points/vtkArrayListTemplate.txx
+++ b/Common/DataModel/vtkArrayListTemplate.txx
@@ -22,9 +22,9 @@
 // Sort of a little object factory (in conjunction w/ vtkTemplateMacro())
 template <typename T>
 void CreateArrayPair(ArrayList *list, T *inData, T *outData,
-                     vtkIdType numPts, int numComp, T nullValue)
+                     vtkIdType numPts, int numComp, vtkDataArray *outArray, T nullValue)
 {
-  ArrayPair<T> *pair = new ArrayPair<T>(inData,outData,numPts,numComp, nullValue);
+  ArrayPair<T> *pair = new ArrayPair<T>(inData,outData,numPts,numComp,outArray,nullValue);
   list->Arrays.push_back(pair);
 }
 
@@ -64,7 +64,7 @@ AddArrays(vtkIdType numOutPts, vtkDataSetAttributes *inPD, vtkDataSetAttributes
             {
             vtkTemplateMacro(CreateArrayPair(this, static_cast<VTK_TT *>(iD),
                              static_cast<VTK_TT *>(oD),numOutPts,oNumComp,
-                                             static_cast<VTK_TT>(nullValue)));
+                             oArray,static_cast<VTK_TT>(nullValue)));
             }//over all VTK types
           }//if matching types
         }//if matching input array
diff --git a/Filters/Core/Testing/Data/Baseline/TestFlyingEdges3DWithInterpolation.png.md5 b/Filters/Core/Testing/Data/Baseline/TestFlyingEdges3DWithInterpolation.png.md5
new file mode 100644
index 0000000000000000000000000000000000000000..2f1e9ed59b639f4a652e18e4402837e8e3fe06b9
--- /dev/null
+++ b/Filters/Core/Testing/Data/Baseline/TestFlyingEdges3DWithInterpolation.png.md5
@@ -0,0 +1 @@
+11b30d52f07417d89c088a3d69bb3d7a
diff --git a/Filters/Core/Testing/Python/CMakeLists.txt b/Filters/Core/Testing/Python/CMakeLists.txt
index 9eb9c78e453fa9f92c925aceb0b02400afbbbe6b..a34a5feb55851ded0c36551edcc3af2a1deb6b5f 100644
--- a/Filters/Core/Testing/Python/CMakeLists.txt
+++ b/Filters/Core/Testing/Python/CMakeLists.txt
@@ -18,6 +18,7 @@ vtk_add_test_python(
   TestElevationFilter.py
   TestFlyingEdges2D.py
   TestFlyingEdges3D.py
+  TestFlyingEdges3DWithInterpolation.py
   TestFlyingEdgesExtents.py
   TestFlyingEdgesPlaneCutter.py
   TestGridSynchronizedTemplates3D.py
diff --git a/Filters/Core/Testing/Python/TestFlyingEdges3D.py b/Filters/Core/Testing/Python/TestFlyingEdges3D.py
index c6da8d826841c4b83ab13707baa2d1de3fb7e8cb..f19b37e48ba54b0ba35f4a042d38d2430b41b545 100644
--- a/Filters/Core/Testing/Python/TestFlyingEdges3D.py
+++ b/Filters/Core/Testing/Python/TestFlyingEdges3D.py
@@ -31,6 +31,8 @@ iso.SetInputConnection(sample.GetOutputPort())
 iso.SetValue(0,0.25)
 iso.ComputeNormalsOn()
 iso.ComputeGradientsOn()
+iso.ComputeScalarsOn()
+iso.InterpolateAttributesOff()
 
 isoMapper = vtk.vtkPolyDataMapper()
 isoMapper.SetInputConnection(iso.GetOutputPort())
diff --git a/Filters/Core/Testing/Python/TestFlyingEdges3DWithInterpolation.py b/Filters/Core/Testing/Python/TestFlyingEdges3DWithInterpolation.py
new file mode 100755
index 0000000000000000000000000000000000000000..ab5d79f231ff5b531254952ed6212b227a2b439c
--- /dev/null
+++ b/Filters/Core/Testing/Python/TestFlyingEdges3DWithInterpolation.py
@@ -0,0 +1,86 @@
+#!/usr/bin/env python
+import vtk
+from vtk.test import Testing
+from vtk.util.misc import vtkGetDataRoot
+VTK_DATA_ROOT = vtkGetDataRoot()
+
+# Create the RenderWindow, Renderer and both Actors
+#
+ren1 = vtk.vtkRenderer()
+renWin = vtk.vtkRenderWindow()
+renWin.SetMultiSamples(0)
+renWin.AddRenderer(ren1)
+iren = vtk.vtkRenderWindowInteractor()
+iren.SetRenderWindow(renWin)
+
+# Create a synthetic source
+sphere = vtk.vtkSphere()
+sphere.SetCenter( 0.0,0.0,0.0)
+sphere.SetRadius(0.25)
+
+# Iso-surface to create geometry. Demonstrate the ability to
+# interpolate data attributes.
+sample = vtk.vtkSampleFunction()
+sample.SetImplicitFunction(sphere)
+sample.SetModelBounds(-0.5,0.5, -0.5,0.5, -0.5,0.5)
+sample.SetSampleDimensions(100,100,100)
+
+# Now create some new attributes
+cyl = vtk.vtkCylinder()
+cyl.SetRadius(0.1)
+cyl.SetAxis(1,1,1)
+
+attr = vtk.vtkSampleImplicitFunctionFilter()
+attr.SetInputConnection(sample.GetOutputPort())
+attr.SetImplicitFunction(cyl)
+attr.ComputeGradientsOn()
+attr.Update()
+
+iso = vtk.vtkFlyingEdges3D()
+iso.SetInputConnection(attr.GetOutputPort())
+iso.SetInputArrayToProcess(0, 0, 0, vtk.vtkDataObject.FIELD_ASSOCIATION_POINTS, "scalars")
+iso.SetValue(0,0.25)
+iso.ComputeNormalsOn()
+iso.ComputeGradientsOn()
+iso.ComputeScalarsOn()
+iso.InterpolateAttributesOn()
+
+# Time execution
+timer = vtk.vtkTimerLog()
+timer.StartTimer()
+iso.Update()
+timer.StopTimer()
+time = timer.GetElapsedTime()
+print("Flying edges with attributes: {0}".format(time))
+
+isoMapper = vtk.vtkPolyDataMapper()
+isoMapper.SetInputConnection(iso.GetOutputPort())
+isoMapper.ScalarVisibilityOn()
+isoMapper.SetScalarModeToUsePointFieldData()
+isoMapper.SelectColorArray("Implicit scalars")
+isoMapper.SetScalarRange(0,.3)
+
+isoActor = vtk.vtkActor()
+isoActor.SetMapper(isoMapper)
+isoActor.GetProperty().SetColor(1,1,1)
+isoActor.GetProperty().SetOpacity(1)
+
+outline = vtk.vtkOutlineFilter()
+outline.SetInputConnection(sample.GetOutputPort())
+outlineMapper = vtk.vtkPolyDataMapper()
+outlineMapper.SetInputConnection(outline.GetOutputPort())
+outlineActor = vtk.vtkActor()
+outlineActor.SetMapper(outlineMapper)
+outlineProp = outlineActor.GetProperty()
+
+# Add the actors to the renderer, set the background and size
+#
+ren1.AddActor(outlineActor)
+ren1.AddActor(isoActor)
+ren1.SetBackground(0,0,0)
+renWin.SetSize(300,300)
+ren1.ResetCamera()
+iren.Initialize()
+
+renWin.Render()
+# --- end of script --
diff --git a/Filters/Core/vtkFlyingEdges3D.cxx b/Filters/Core/vtkFlyingEdges3D.cxx
index 3ad4df0a3eaacb5ca609080edf83c2332f3a94aa..cf353f2185ebf7c06149456448cedc74389c23da 100644
--- a/Filters/Core/vtkFlyingEdges3D.cxx
+++ b/Filters/Core/vtkFlyingEdges3D.cxx
@@ -26,6 +26,7 @@
 #include "vtkFloatArray.h"
 #include "vtkStreamingDemandDrivenPipeline.h"
 #include "vtkMarchingCubesTriangleCases.h"
+#include "vtkArrayListTemplate.h"
 #include "vtkSMPTools.h"
 
 #include <cmath>
@@ -33,6 +34,9 @@
 vtkStandardNewMacro(vtkFlyingEdges3D);
 
 //----------------------------------------------------------------------------
+namespace {
+#include "vtkArrayListTemplate.h" // For processing attribute data
+}
 
 // This templated class implements the heart of the algorithm.
 // vtkFlyingEdges3D populates the information in this class and
@@ -121,6 +125,8 @@ public:
   float     *NewGradients;
   float     *NewNormals;
   bool       NeedGradients;
+  bool       InterpolateAttributes;
+  ArrayList  Arrays;
 
   // Setup algorithm
   vtkFlyingEdges3DAlgorithm();
@@ -212,7 +218,8 @@ public:
                            const int incs[3],
                            float x1[3],
                            vtkIdType vId,
-                           vtkIdType ijk[3],
+                           vtkIdType ijk0[3],
+                           vtkIdType ijk1[3],
                            float g0[3])
     {
 
@@ -224,7 +231,7 @@ public:
       if ( this->NeedGradients )
         {
         float gTmp[3], g1[3];
-        this->ComputeGradient(loc,ijk,
+        this->ComputeGradient(loc,ijk1,
                               s + incs[0], s - incs[0],
                               s + incs[1], s - incs[1],
                               s + incs[2], s - incs[2],
@@ -244,6 +251,13 @@ public:
           vtkMath::Normalize(n);
           }
         }//if normals or gradients required
+
+      if ( this->InterpolateAttributes )
+        {
+        vtkIdType v0=ijk0[0] + ijk0[1]*incs[1] + ijk0[2]*incs[2];
+        vtkIdType v1=ijk1[0] + ijk1[1]*incs[1] + ijk1[2]*incs[2];;
+        this->Arrays.InterpolateEdge(v0,v1,t,vId);
+        }
     }
 
   // Compute the gradient on a point which may be on the boundary of the volume.
@@ -380,7 +394,7 @@ public:
   // Interface between VTK and templated functions
   static void Contour(vtkFlyingEdges3D *self, vtkImageData *input,
                       int extent[6], vtkIdType *incs, T *scalars,
-                      vtkPoints *newPts, vtkCellArray *newTris,
+                      vtkPolyData *output, vtkPoints *newPts, vtkCellArray *newTris,
                       vtkDataArray *newScalars,vtkFloatArray *newNormals,
                       vtkFloatArray *newGradients);
 };
@@ -683,6 +697,13 @@ InterpolateEdge(double value, vtkIdType ijk[3],
       vtkMath::Normalize(n);
       }
     }//if normals or gradients required
+
+  if ( this->InterpolateAttributes )
+    {
+    vtkIdType v0=ijk0[0] + ijk0[1]*incs[1] + ijk0[2]*incs[2];
+    vtkIdType v1=ijk1[0] + ijk1[1]*incs[1] + ijk1[2]*incs[2];;
+    this->Arrays.InterpolateEdge(v0,v1,t,vId);
+    }
 }
 
 //----------------------------------------------------------------------------
@@ -719,7 +740,7 @@ GeneratePoints(double value, unsigned char loc, vtkIdType ijk[3],
 
       T const * const sPtr2 = (sPtr+incs[i]);
       double t = (value - *sPtr) / (*sPtr2 - *sPtr);
-      this->InterpolateAxesEdge(t, loc, x, sPtr2, incs, x1, eIds[i*4], ijk1, g0);
+      this->InterpolateAxesEdge(t, loc, x, sPtr2, incs, x1, eIds[i*4], ijk, ijk1, g0);
       }
     }
 
@@ -794,15 +815,14 @@ ProcessXEdge(double value, T const* const inPtr, vtkIdType row, vtkIdType slice)
   vtkIdType nxcells=this->Dims[0]-1;
   vtkIdType minInt=nxcells, maxInt = 0;
   vtkIdType *edgeMetaData;
-  unsigned char *ePtr = this->XCases + slice*this->SliceOffset + row*nxcells;
+  unsigned char edgeCase, *ePtr=this->XCases+slice*this->SliceOffset+row*nxcells;
   double s0, s1 = static_cast<double>(*inPtr);
+  vtkIdType sum = 0;
 
   //run along the entire x-edge computing edge cases
   edgeMetaData = this->EdgeMetaData + (slice*this->Dims[1] + row)*6;
   std::fill_n(edgeMetaData, 6, 0);
 
-  vtkIdType sum = 0;
-
   //pull this out help reduce false sharing
   vtkIdType inc0 = this->Inc0;
 
@@ -811,11 +831,14 @@ ProcessXEdge(double value, T const* const inPtr, vtkIdType row, vtkIdType slice)
     s0 = s1;
     s1 = static_cast<double>(*(inPtr + (i+1)*inc0));
 
-    unsigned char edgeCase = vtkFlyingEdges3DAlgorithm::Below;
     if (s0 >= value)
       {
       edgeCase = vtkFlyingEdges3DAlgorithm::LeftAbove;
       }
+    else
+      {
+      edgeCase = vtkFlyingEdges3DAlgorithm::Below;
+      }
     if( s1 >= value)
       {
       edgeCase |= vtkFlyingEdges3DAlgorithm::RightAbove;
@@ -828,7 +851,10 @@ ProcessXEdge(double value, T const* const inPtr, vtkIdType row, vtkIdType slice)
          edgeCase == vtkFlyingEdges3DAlgorithm::RightAbove )
       {
       ++sum; //increment number of intersections along x-edge
-      minInt = ( i < minInt ? i : minInt);
+      if ( i < minInt )
+        {
+        minInt = i;
+        }
       maxInt = i + 1;
       }//if contour interacts with this x-edge
     }//for all x-cell edges along this x-edge
@@ -936,6 +962,7 @@ ProcessYZEdges(vtkIdType row, vtkIdType slice)
   // voxel axes. Also check the number of primitives generated.
   unsigned char *edgeUses, eCase, numTris;
   ePtr[0] += xL; ePtr[1] += xL; ePtr[2] += xL; ePtr[3] += xL;
+  const vtkIdType dim0Wall = this->Dims[0]-2;
   for (i=xL; i < xR; ++i) //run along the trimmed x-voxels
     {
     eCase = this->GetEdgeCase(ePtr);
@@ -950,7 +977,7 @@ ProcessYZEdges(vtkIdType row, vtkIdType slice)
       edgeUses = this->GetEdgeUses(eCase);
       eMD[0][1] += edgeUses[4]; //y-voxel axes edge always counted
       eMD[0][2] += edgeUses[8]; //z-voxel axes edge always counted
-      loc = yzLoc | (i >= (this->Dims[0]-2) ? MaxBoundary : Interior);
+      loc = yzLoc | (i >= dim0Wall ? MaxBoundary : Interior);
       if ( loc != 0 )
         {
         this->CountBoundaryYZInts(loc,edgeUses,eMD);
@@ -1030,6 +1057,8 @@ GenerateOutput(double value, T* rowPtr, vtkIdType row, vtkIdType slice)
   //load the inc0/inc1/inc2 into local memory
   const int incs[3] = { this->Inc0, this->Inc1, this->Inc2 };
   const T* sPtr = rowPtr + xL*incs[0];
+  const double xSpace = this->Spacing[0];
+  const vtkIdType dim0Wall = this->Dims[0]-2;
 
   for (i=xL; i < xR; ++i)
     {
@@ -1042,7 +1071,7 @@ GenerateOutput(double value, T* rowPtr, vtkIdType row, vtkIdType slice)
       // Now generate point(s) along voxel axes if needed. Remember to take
       // boundary into account.
       loc = yzLoc | (i < 1 ? MinBoundary :
-          (i >= (this->Dims[0]-2) ? MaxBoundary : Interior));
+                     (i >= dim0Wall ? MaxBoundary : Interior));
       if ( this->CaseIncludesAxes(eCase) || loc != Interior )
         {
         unsigned char const * const edgeUses = this->GetEdgeUses(eCase);
@@ -1059,7 +1088,7 @@ GenerateOutput(double value, T* rowPtr, vtkIdType row, vtkIdType slice)
 
     ++ijk[0];
     sPtr += incs[0];
-    x[0]+= this->Spacing[0];
+    x[0] += xSpace;
     } //for all non-trimmed cells along this x-edge
 }
 
@@ -1069,9 +1098,9 @@ GenerateOutput(double value, T* rowPtr, vtkIdType row, vtkIdType slice)
 // class. It also invokes the three passes of the Flying Edges algorithm.
 template <class T> void vtkFlyingEdges3DAlgorithm<T>::
 Contour(vtkFlyingEdges3D *self, vtkImageData *input, int extent[6],
-        vtkIdType *incs, T *scalars, vtkPoints *newPts, vtkCellArray *newTris,
-        vtkDataArray *newScalars, vtkFloatArray *newNormals,
-        vtkFloatArray *newGradients)
+        vtkIdType *incs, T *scalars, vtkPolyData *output, vtkPoints *newPts,
+        vtkCellArray *newTris, vtkDataArray *newScalars,
+        vtkFloatArray *newNormals, vtkFloatArray *newGradients)
 {
   double value, *values = self->GetValues();
   int numContours = self->GetNumberOfContours();
@@ -1114,6 +1143,12 @@ Contour(vtkFlyingEdges3D *self, vtkImageData *input, int extent[6],
   // for computational trimming).
   algo.EdgeMetaData = new vtkIdType [algo.NumberOfEdges*6];
 
+  // Interpolating attributes and other stuff. Interpolate extra attributes only if they
+  // exist and the user requests it.
+  algo.NeedGradients = (newGradients || newNormals);
+  algo.InterpolateAttributes = (self->GetInterpolateAttributes() &&
+                                 input->GetPointData()->GetNumberOfArrays() > 1) ? true : false;
+
   // Loop across each contour value. This encompasses all three passes.
   for (vidx = 0; vidx < numContours; vidx++)
     {
@@ -1195,7 +1230,19 @@ Contour(vtkFlyingEdges3D *self, vtkImageData *input, int extent[6],
         newNormals->WriteVoidPointer(0,3*totalPts);
         algo.NewNormals = static_cast<float*>(newNormals->GetVoidPointer(0));
         }
-      algo.NeedGradients = (algo.NewGradients || algo.NewNormals);
+      if ( algo.InterpolateAttributes )
+        {
+        if ( vidx == 0 ) //first contour
+          {
+//          output->GetPointData()->CopyScalarsOff();
+          output->GetPointData()->InterpolateAllocate(input->GetPointData(),totalPts);
+          algo.Arrays.AddArrays(totalPts,input->GetPointData(),output->GetPointData());
+          }
+        else
+          {
+          algo.Arrays.Realloc(totalPts);
+          }
+        }
 
       // PASS 4: Fourth and final pass: Process voxel rows and generate output.
       // Note that we are simultaneously generating triangles and interpolating
@@ -1226,6 +1273,7 @@ vtkFlyingEdges3D::vtkFlyingEdges3D()
   this->ComputeNormals = 1;
   this->ComputeGradients = 0;
   this->ComputeScalars = 1;
+  this->InterpolateAttributes = 0;
   this->ArrayComponent = 0;
 
   // by default process active point scalars
@@ -1365,7 +1413,7 @@ int vtkFlyingEdges3D::RequestData(
     {
     vtkTemplateMacro(vtkFlyingEdges3DAlgorithm<VTK_TT>::
                      Contour(this, input, exExt, incs, (VTK_TT *)ptr,
-                             newPts, newTris, newScalars, newNormals,
+                             output, newPts, newTris, newScalars, newNormals,
                              newGradients));
     }
 
@@ -1422,5 +1470,6 @@ void vtkFlyingEdges3D::PrintSelf(ostream& os, vtkIndent indent)
   os << indent << "Compute Normals: " << (this->ComputeNormals ? "On\n" : "Off\n");
   os << indent << "Compute Gradients: " << (this->ComputeGradients ? "On\n" : "Off\n");
   os << indent << "Compute Scalars: " << (this->ComputeScalars ? "On\n" : "Off\n");
+  os << indent << "Interpolate Attributes: " << (this->InterpolateAttributes ? "On\n" : "Off\n");
   os << indent << "ArrayComponent: " << this->ArrayComponent << endl;
 }
diff --git a/Filters/Core/vtkFlyingEdges3D.h b/Filters/Core/vtkFlyingEdges3D.h
index b3ddafa858c70858a37253ac7ffa8b258fbf1fd2..a0cd740d8ad5fed59598697cb0e4256eab0bee39 100644
--- a/Filters/Core/vtkFlyingEdges3D.h
+++ b/Filters/Core/vtkFlyingEdges3D.h
@@ -100,6 +100,15 @@ public:
   vtkGetMacro(ComputeScalars,int);
   vtkBooleanMacro(ComputeScalars,int);
 
+  // Description:
+  // Indicate whether to interpolate other attribute data. That is, as the
+  // isosurface is generated, interpolate all point attribute data across
+  // the edge. This is independent of scalar interpolation, which is
+  // controlled by the ComputeScalars flag.
+  vtkSetMacro(InterpolateAttributes,int);
+  vtkGetMacro(InterpolateAttributes,int);
+  vtkBooleanMacro(InterpolateAttributes,int);
+
   // Description:
   // Set a particular contour value at contour number i. The index i ranges
   // between 0<=i<NumberOfContours.
@@ -157,6 +166,7 @@ protected:
   int ComputeNormals;
   int ComputeGradients;
   int ComputeScalars;
+  int InterpolateAttributes;
   int ArrayComponent;
   vtkContourValues *ContourValues;
 
diff --git a/Filters/Points/CMakeLists.txt b/Filters/Points/CMakeLists.txt
index 4b3e6508dbed3334a92e750cc49f9a19af2336a1..c2fc0c25e4870950713e83c6ad3f68c174c5122b 100644
--- a/Filters/Points/CMakeLists.txt
+++ b/Filters/Points/CMakeLists.txt
@@ -1,5 +1,4 @@
 set(Module_SRCS
-  vtkArrayListTemplate.txx
   vtkEllipsoidalGaussianKernel.cxx
   vtkGaussianKernel.cxx
   vtkInterpolationKernel.cxx
@@ -10,17 +9,9 @@ set(Module_SRCS
   vtkVoronoiKernel.cxx
   )
 
-set(${vtk-module}_HDRS
-  vtkArrayListTemplate.h
-  )
-
 set_source_files_properties(
   vtkInterpolationKernel
   ABSTRACT
   )
 
-set_source_files_properties(
-  vtkArrayListTemplate.txx
-  WRAP_EXCLUDE
-)
 vtk_module_library(vtkFiltersPoints ${Module_SRCS})