diff --git a/Imaging/Stencil/vtkPolyDataToImageStencil.cxx b/Imaging/Stencil/vtkPolyDataToImageStencil.cxx
index b7cf8c8..99a4ea0 100644
--- a/Imaging/Stencil/vtkPolyDataToImageStencil.cxx
+++ b/Imaging/Stencil/vtkPolyDataToImageStencil.cxx
@@ -53,7 +53,6 @@ POSSIBILITY OF SUCH DAMAGES.
 #include "vtkCellArray.h"
 #include "vtkDoubleArray.h"
 #include "vtkSignedCharArray.h"
-#include "vtkMergePoints.h"
 #include "vtkPoints.h"
 #include "vtkPointData.h"
 #include "vtkCellData.h"
@@ -64,6 +63,8 @@ POSSIBILITY OF SUCH DAMAGES.
 #include "vtkInformationVector.h"
 #include "vtkStreamingDemandDrivenPipeline.h"
 
+#include <map>
+#include <utility>
 #include <math.h>
 
 
@@ -109,10 +110,191 @@ void vtkPolyDataToImageStencil::PrintSelf(ostream& os,
 }
 
 //----------------------------------------------------------------------------
+// A helper class to quickly locate an edge, given the endpoint ids.
+// It uses an stl map rather than a table partitioning scheme, since
+// we have no idea how many entries there will be when we start.  So
+// the performance is approximately log(n).
+namespace {
+
+// A Node in a linked list that contains information about one edge
+class EdgeLocatorNode
+{
+public:
+  EdgeLocatorNode() :
+    ptId0(-1), ptId1(-1), edgeId(-1), next(0) {}
+
+  // Free the list that this node is the head of
+  void FreeList() {
+    EdgeLocatorNode *ptr = this->next;
+    while (ptr)
+      {
+      EdgeLocatorNode *tmp = ptr;
+      ptr = ptr->next;
+      tmp->next = 0;
+      delete tmp;
+      }
+  }
+
+  vtkIdType ptId0;
+  vtkIdType ptId1;
+  vtkIdType edgeId;
+  EdgeLocatorNode *next;
+};
+
+// The EdgeLocator class itself, for keeping track of edges
+class EdgeLocator
+{
+private:
+  typedef std::map<vtkIdType, EdgeLocatorNode> MapType;
+  MapType EdgeMap;
+
+public:
+  EdgeLocator() : EdgeMap() {}
+  ~EdgeLocator() { this->Initialize(); }
+
+  // Description:
+  // Initialize the locator.
+  void Initialize();
+
+  // Description:
+  // If edge (i0, i1) is not in the list, then it will be added and
+  // a pointer for storing the new edgeId will be returned.
+  // If edge (i0, i1) is in the list, then edgeId will be set to the
+  // stored value and a null pointer will be returned.
+  vtkIdType *InsertUniqueEdge(vtkIdType i0, vtkIdType i1, vtkIdType &edgeId);
+
+  // Description:
+  // A helper function for interpolating a new point along an edge.  It
+  // stores the index of the interpolated point in "i", and returns 1 if
+  // a new point was added to the locator.  The values i0, i1, v0, v1 are
+  // the edge endpoints and scalar values, respectively.
+  int InterpolateEdge(
+    vtkPoints *inPoints, vtkPoints *outPoints,
+    vtkIdType i0, vtkIdType i1, double v0, double v1,
+    vtkIdType &i);
+};
+
+void EdgeLocator::Initialize()
+{
+  for (MapType::iterator i = this->EdgeMap.begin();
+       i != this->EdgeMap.end();
+       ++i)
+    {
+    i->second.FreeList();
+    }
+  this->EdgeMap.clear();
+}
+
+vtkIdType *EdgeLocator::InsertUniqueEdge(
+  vtkIdType i0, vtkIdType i1, vtkIdType &edgeId)
+{
+  // Ensure consistent ordering of edge
+  if (i1 < i0)
+    {
+    vtkIdType tmp = i0;
+    i0 = i1;
+    i1 = tmp;
+    }
+
+  // Generate a integer key, try to make it unique
+#ifdef VTK_USE_64BIT_IDS
+  vtkIdType key = ((i1 << 32) ^ i0);
+#else
+  vtkIdType key = ((i1 << 16) ^ i0);
+#endif
+
+  EdgeLocatorNode *node = &this->EdgeMap[key];
+
+  if (node->ptId1 < 0)
+    {
+    // Didn't find key, so add a new edge entry
+    node->ptId0 = i0;
+    node->ptId1 = i1;
+    return &node->edgeId;
+    }
+
+  // Search through the list for i0 and i1
+  if (node->ptId0 == i0 && node->ptId1 == i1)
+    {
+    edgeId = node->edgeId;
+    return 0;
+    }
+
+  int i = 1;
+  while (node->next != 0)
+    {
+    i++;
+    node = node->next;
+
+    if (node->ptId0 == i0 && node->ptId1 == i1)
+      {
+      edgeId = node->edgeId;
+      return 0;
+      }
+    }
+
+  // No entry for i1, so make one and return
+  node->next = new EdgeLocatorNode;
+  node = node->next;
+  node->ptId0 = i0;
+  node->ptId1 = i1;
+  node->edgeId = this->EdgeMap.size()-1;
+  return &node->edgeId;
+}
+
+int EdgeLocator::InterpolateEdge(
+  vtkPoints *points, vtkPoints *outPoints,
+  vtkIdType i0, vtkIdType i1, double v0, double v1,
+  vtkIdType &i)
+{
+  // This swap guarantees that exactly the same point is computed
+  // for both line directions, as long as the endpoints are the same.
+  if (v1 > 0)
+    {
+    vtkIdType tmpi = i0;
+    i0 = i1;
+    i1 = tmpi;
+
+    double tmp = v0;
+    v0 = v1;
+    v1 = tmp;
+    }
+
+  // After the above swap, i0 will be kept, and i1 will be clipped
+
+  // Check to see if this point has already been computed
+  vtkIdType *iptr = this->InsertUniqueEdge(i0, i1, i);
+  if (iptr == 0)
+    {
+    return 0;
+    }
+
+  // Get the edge and interpolate the new point
+  double p0[3], p1[3], p[3];
+  points->GetPoint(i0, p0);
+  points->GetPoint(i1, p1);
+
+  double f = v0/(v0 - v1);
+  double s = 1.0 - f;
+  double t = 1.0 - s;
+
+  p[0] = s*p0[0] + t*p1[0];
+  p[1] = s*p0[1] + t*p1[1];
+  p[2] = s*p0[2] + t*p1[2];
+
+  // Add the point, store the new index in the locator
+  i = outPoints->InsertNextPoint(p);
+  *iptr = i;
+
+  return 1;
+}
+
+} // end anonymous namespace
+
+//----------------------------------------------------------------------------
 // Select contours within slice z
 void vtkPolyDataToImageStencil::PolyDataSelector(
-  vtkPolyData *input, vtkPolyData *output, double z, double thickness,
-  vtkMergePoints *locator)
+  vtkPolyData *input, vtkPolyData *output, double z, double thickness)
 {
   vtkPoints *points = input->GetPoints();
   vtkCellArray *lines = input->GetLines();
@@ -124,11 +306,8 @@ void vtkPolyDataToImageStencil::PolyDataSelector(
   double minz = z - 0.5*thickness;
   double maxz = z + 0.5*thickness;
 
-  double bounds[6];
-  input->GetBounds(bounds);
-  bounds[4] = minz;
-  bounds[5] = maxz;
-  locator->InitPointInsertion(newPoints, bounds);
+  // use a map to avoid adding duplicate points
+  std::map<vtkIdType, vtkIdType> pointLocator;
 
   vtkIdType loc = 0;
   vtkIdType numCells = lines->GetNumberOfCells();
@@ -155,10 +334,21 @@ void vtkPolyDataToImageStencil::PolyDataSelector(
     newLines->InsertNextCell(npts);
     for (i = 0; i < npts; i++)
       {
-      vtkIdType ptId;
-      double point[3];
-      points->GetPoint(ptIds[i], point);
-      locator->InsertUniquePoint(point, ptId);
+      vtkIdType oldId = ptIds[i];
+      std::map<vtkIdType, vtkIdType>::iterator iter =
+        pointLocator.lower_bound(oldId);
+      vtkIdType ptId = 0;
+      if (iter == pointLocator.end() || iter->first != oldId)
+        {
+        double point[3];
+        points->GetPoint(oldId, point);
+        ptId = newPoints->InsertNextPoint(point);
+        pointLocator.insert(iter, std::make_pair(oldId, ptId));
+        }
+      else
+        {
+        ptId = iter->second;
+        }
       newLines->InsertCellPoint(ptId);
       }
     }
@@ -172,75 +362,100 @@ void vtkPolyDataToImageStencil::PolyDataSelector(
 //----------------------------------------------------------------------------
 // This method was taken from vtkCutter and slightly modified
 void vtkPolyDataToImageStencil::PolyDataCutter(
-  vtkPolyData *input, vtkPolyData *output, double z,
-  vtkMergePoints *locator)
+  vtkPolyData *input, vtkPolyData *output, double z)
 {
-  vtkCellData *inCD = input->GetCellData();
-  vtkCellData *outCD = output->GetCellData();
-  vtkDoubleArray *cellScalars = vtkDoubleArray::New();
-
-  // For the new points and lines
+  vtkPoints *points = input->GetPoints();
+  vtkCellArray *inputPolys = input->GetPolys();
+  vtkCellArray *inputStrips = input->GetStrips();
   vtkPoints *newPoints = vtkPoints::New();
   newPoints->Allocate(333);
   vtkCellArray *newLines = vtkCellArray::New();
   newLines->Allocate(1000);
 
-  // No verts or polys are expected
-  vtkCellArray *newVerts = vtkCellArray::New();
-  vtkCellArray *newPolys = vtkCellArray::New();
-
-  // Allocate space for the cell data
-  outCD->CopyAllocate(inCD, 1000);
+  // An edge locator to avoid point duplication while clipping
+  EdgeLocator edgeLocator;
 
-  // locator used to merge potentially duplicate points
-  locator->InitPointInsertion(newPoints, input->GetBounds());
+  // Go through all cells and clip them.
+  vtkIdType numPolys = input->GetNumberOfPolys();
+  vtkIdType numStrips = input->GetNumberOfStrips();
+  vtkIdType numCells = numPolys + numStrips;
 
-  // Compute some information for progress methods
-  vtkGenericCell *cell = vtkGenericCell::New();
-
-  // Loop over all cells; get scalar values for all cell points
-  // and process each cell.
-  vtkIdType numCells = input->GetNumberOfCells();
+  vtkIdType loc = 0;
+  vtkCellArray *cellArray = inputPolys;
   for (vtkIdType cellId = 0; cellId < numCells; cellId++)
     {
-    input->GetCell(cellId, cell);
-    vtkPoints *cellPts = cell->GetPoints();
-    vtkIdList *cellIds = cell->GetPointIds();
+    // switch to strips when polys are done
+    if (cellId == numPolys)
+      {
+      loc = 0;
+      cellArray = inputStrips;
+      }
 
-    if (cell->GetCellDimension() == 2)
+    vtkIdType npts, *ptIds;
+    cellArray->GetCell(loc, npts, ptIds);
+    loc += npts + 1;
+
+    vtkIdType numSubCells = 1;
+    if (cellArray == inputStrips)
+      {
+      numSubCells = npts - 2;
+      npts = 3;
+      }
+
+    for (vtkIdType subId = 0; subId < numSubCells; subId++)
       {
-      vtkIdType numCellPts = cellPts->GetNumberOfPoints();
-      cellScalars->SetNumberOfTuples(numCellPts);
-      for (vtkIdType i = 0; i < numCellPts; i++)
+      vtkIdType i1 = ptIds[npts-1];
+      double point[3];
+      points->GetPoint(i1, point);
+      double v1 = point[2] - z;
+      bool c1 = (v1 > 0);
+      bool odd = ((subId & 1) != 0);
+
+      // To store the ids of the contour line
+      vtkIdType linePts[2];
+      linePts[0] = 0;
+      linePts[1] = 0;
+
+      for (vtkIdType i = 0; i < npts; i++)
+        {
+        // Save previous point info
+        vtkIdType i0 = i1;
+        double v0 = v1;
+        bool c0 = c1;
+
+        // Generate new point info
+        i1 = ptIds[i];
+        points->GetPoint(i1, point);
+        v1 = point[2] - z;
+        c1 = (v1 > 0);
+
+        // If at least one edge end point wasn't clipped
+        if ( (c0 | c1) )
+          {
+          // If only one end was clipped, interpolate new point
+          if ( (c0 ^ c1) )
+            {
+            edgeLocator.InterpolateEdge(
+              points, newPoints, i0, i1, v0, v1, linePts[c0 ^ odd]);
+            }
+          }
+        }
+
+      // Insert the contour line if one was created
+      if (linePts[0] != linePts[1])
         {
-        // scalar value is distance from the specified z plane
-        cellScalars->SetValue(i, input->GetPoint(cellIds->GetId(i))[2]);
+        newLines->InsertNextCell(2, linePts);
         }
 
-      cell->Contour(z, cellScalars, locator,
-                    newVerts, newLines, newPolys, NULL, NULL,
-                    inCD, cellId, outCD);
+      // Increment to get to the next triangle, if cell is a strip
+      ptIds++;
       }
     }
 
-  // Update ourselves.  Because we don't know upfront how many verts, lines,
-  // polys we've created, take care to reclaim memory.
-  cell->Delete();
-  cellScalars->Delete();
-
   output->SetPoints(newPoints);
+  output->SetLines(newLines);
   newPoints->Delete();
-
-  if (newLines->GetNumberOfCells())
-    {
-    output->SetLines(newLines);
-    }
   newLines->Delete();
-  newVerts->Delete();
-  newPolys->Delete();
-
-  //release any extra memory
-  locator->Initialize();
 }
 
 //----------------------------------------------------------------------------
@@ -279,9 +494,6 @@ void vtkPolyDataToImageStencil::ThreadedExecute(
   // get the input data
   vtkPolyData *input = this->GetInput();
 
-  // the locator to use with the data
-  vtkMergePoints *locator = vtkMergePoints::New();
-
   // the output produced by cutting the polydata with the Z plane
   vtkPolyData *slice = vtkPolyData::New();
 
@@ -312,12 +524,12 @@ void vtkPolyDataToImageStencil::ThreadedExecute(
     // Step 1: Cut the data into slices
     if (input->GetNumberOfPolys() > 0 || input->GetNumberOfStrips() > 0)
       {
-      this->PolyDataCutter(input, slice, z, locator);
+      this->PolyDataCutter(input, slice, z);
       }
     else
       {
       // if no polys, select polylines instead
-      this->PolyDataSelector(input, slice, z, spacing[2], locator);
+      this->PolyDataSelector(input, slice, z, spacing[2]);
       }
 
     if (!slice->GetNumberOfLines())
@@ -360,7 +572,7 @@ void vtkPolyDataToImageStencil::ThreadedExecute(
       lines->InitTraversal();
       vtkIdType npts;
       vtkIdType *pointIds;
-      while( lines->GetNextCell(npts, pointIds) )
+      while (lines->GetNextCell(npts, pointIds))
         {
         int isClosed = 0;
         if (pointIds[0] == pointIds[npts-1])
@@ -395,82 +607,114 @@ void vtkPolyDataToImageStencil::ThreadedExecute(
       if (numberOfNeighbors == 1)
         {
         // store the loose end
-        looseEndIdList->InsertNextId( i );
-        looseEndNeighborList->InsertNextId( neighborId );
+        looseEndIdList->InsertNextId(i);
+        looseEndNeighborList->InsertNextId(neighborId);
         }
       // mark inflection points
-      inflectionPointList->InsertNextValue( bottomPoint | topPoint );
+      inflectionPointList->InsertNextValue(bottomPoint | topPoint);
       }
 
+    // join any loose ends
     while (looseEndIdList->GetNumberOfIds() >= 2)
       {
-      // first loose end point in the list
-      vtkIdType firstLooseEndId = looseEndIdList->GetId(0);
-      vtkIdType neighborId = looseEndNeighborList->GetId(0);
-      double firstLooseEnd[3];
-      slice->GetPoint( firstLooseEndId, firstLooseEnd );
-      double neighbor[3];
-      slice->GetPoint( neighborId, neighbor);
-
-      // second loose end in the list
-      vtkIdType secondLooseEndId = looseEndIdList->GetId(1);
-      double secondLooseEnd[3];
-      slice->GetPoint( secondLooseEndId, secondLooseEnd );
+      vtkIdType n = looseEndIdList->GetNumberOfIds();
 
-      // search for the loose end closest to the first one
+      // search for the two closest loose ends
       double maxval = -VTK_FLOAT_MAX;
+      vtkIdType firstIndex = 0;
+      vtkIdType secondIndex = 1;
+      bool isCoincident = false;
 
-      for(vtkIdType j = 1; j < looseEndIdList->GetNumberOfIds(); j++)
+      for (vtkIdType i = 0; i < n && !isCoincident; i++)
         {
-        vtkIdType currentLooseEndId = looseEndIdList->GetId( j );
-        if (currentLooseEndId != neighborId)
+        // first loose end
+        vtkIdType firstLooseEndId = looseEndIdList->GetId(i);
+        vtkIdType neighborId = looseEndNeighborList->GetId(i);
+        double firstLooseEnd[3];
+        slice->GetPoint(firstLooseEndId, firstLooseEnd);
+        double neighbor[3];
+        slice->GetPoint(neighborId, neighbor);
+
+        for (vtkIdType j = i+1; j < n && !isCoincident; j++)
           {
-          double currentLooseEnd[3];
-          slice->GetPoint( currentLooseEndId, currentLooseEnd );
-
-          // When connecting loose ends, use dot product to favor
-          // continuing in same direction as the line already
-          // connected to the loose end, but also favour short
-          // distances by dividing dotprod by square of distance.
-          double v1[2], v2[2];
-          v1[0] = firstLooseEnd[0] - neighbor[0];
-          v1[1] = firstLooseEnd[1] - neighbor[1];
-          v2[0] = currentLooseEnd[0] - firstLooseEnd[0];
-          v2[1] = currentLooseEnd[1] - firstLooseEnd[1];
-          double dotprod = v1[0]*v2[0] + v1[1]*v2[1];
-          double distance2 = v2[0]*v2[0] + v2[1]*v2[1];
-
-          if (dotprod > maxval*distance2 && distance2 > 0.0)
+          vtkIdType secondLooseEndId = looseEndIdList->GetId(j);
+          if (secondLooseEndId != neighborId)
             {
-            maxval = dotprod/distance2;
-            secondLooseEndId = currentLooseEndId;
+            double currentLooseEnd[3];
+            slice->GetPoint(secondLooseEndId, currentLooseEnd);
+
+            // When connecting loose ends, use dot product to favor
+            // continuing in same direction as the line already
+            // connected to the loose end, but also favour short
+            // distances by dividing dotprod by square of distance.
+            double v1[2], v2[2];
+            v1[0] = firstLooseEnd[0] - neighbor[0];
+            v1[1] = firstLooseEnd[1] - neighbor[1];
+            v2[0] = currentLooseEnd[0] - firstLooseEnd[0];
+            v2[1] = currentLooseEnd[1] - firstLooseEnd[1];
+            double dotprod = v1[0]*v2[0] + v1[1]*v2[1];
+            double distance2 = v2[0]*v2[0] + v2[1]*v2[1];
+
+            if (distance2 == 0)
+              {
+              firstIndex = i;
+              secondIndex = j;
+              isCoincident = true;
+              }
+            else if (dotprod > maxval*distance2)
+              {
+              firstIndex = i;
+              secondIndex = j;
+              maxval = dotprod/distance2;
+              }
             }
           }
         }
 
-      // create a new line segment by connecting these two points
-      looseEndIdList->DeleteId( firstLooseEndId );
-      looseEndIdList->DeleteId( secondLooseEndId );
-      looseEndNeighborList->DeleteId( firstLooseEndId );
-      looseEndNeighborList->DeleteId( secondLooseEndId );
-
-      lines->InsertNextCell( 2 );
-      lines->InsertCellPoint( firstLooseEndId );
-      lines->InsertCellPoint( secondLooseEndId );
+      // get info about the two loose ends and their neighbors
+      vtkIdType firstLooseEndId = looseEndIdList->GetId(firstIndex);
+      vtkIdType neighborId = looseEndNeighborList->GetId(firstIndex);
+      double firstLooseEnd[3];
+      slice->GetPoint(firstLooseEndId, firstLooseEnd);
+      double neighbor[3];
+      slice->GetPoint(neighborId, neighbor);
 
-      // check if the new vertices are vertical inflection points
-      slice->GetPoint( secondLooseEndId, secondLooseEnd );
-      vtkIdType secondNeighborId = looseEndNeighborList->GetId(0);
+      vtkIdType secondLooseEndId = looseEndIdList->GetId(secondIndex);
+      vtkIdType secondNeighborId = looseEndNeighborList->GetId(secondIndex);
+      double secondLooseEnd[3];
+      slice->GetPoint(secondLooseEndId, secondLooseEnd);
       double secondNeighbor[3];
-      slice->GetPoint( secondNeighborId, secondNeighbor);
+      slice->GetPoint(secondNeighborId, secondNeighbor);
 
-      inflectionPointList->SetValue(firstLooseEndId,
-        ((firstLooseEnd[1] - neighbor[1])*
-        (secondLooseEnd[1] - firstLooseEnd[1]) <= 0));
+      // remove these loose ends from the list
+      looseEndIdList->DeleteId(firstLooseEndId);
+      looseEndIdList->DeleteId(secondLooseEndId);
+      looseEndNeighborList->DeleteId(firstLooseEndId);
+      looseEndNeighborList->DeleteId(secondLooseEndId);
 
-      inflectionPointList->SetValue(secondLooseEndId,
-        ((secondLooseEnd[1] - firstLooseEnd[1])*
-         (secondNeighbor[1] - secondLooseEnd[1]) <= 0));
+      if (isCoincident)
+        {
+        // check if the new vertex is a vertical inflection points
+        inflectionPointList->SetValue(firstLooseEndId,
+          ((firstLooseEnd[1] - neighbor[1])*
+          (secondNeighbor[1] - firstLooseEnd[1]) <= 0));
+        }
+      else
+        {
+        // create a new line segment by connecting these two points
+        lines->InsertNextCell(2);
+        lines->InsertCellPoint(firstLooseEndId);
+        lines->InsertCellPoint(secondLooseEndId);
+
+        // check if the new vertices are vertical inflection points
+        inflectionPointList->SetValue(firstLooseEndId,
+          ((firstLooseEnd[1] - neighbor[1])*
+          (secondLooseEnd[1] - firstLooseEnd[1]) <= 0));
+
+        inflectionPointList->SetValue(secondLooseEndId,
+          ((secondLooseEnd[1] - firstLooseEnd[1])*
+           (secondNeighbor[1] - secondLooseEnd[1]) <= 0));
+        }
       }
 
     // Step 3: Go through all the line segments for this slice,
@@ -480,7 +724,7 @@ void vtkPolyDataToImageStencil::ThreadedExecute(
     vtkIdType *pts = 0;
     vtkIdType npts = 0;
 
-    while ( lines->GetNextCell(npts, pts) )
+    while (lines->GetNextCell(npts, pts))
       {
       for (vtkIdType j = 1; j < npts; j++)
         {
@@ -507,7 +751,6 @@ void vtkPolyDataToImageStencil::ThreadedExecute(
     }
 
   slice->Delete();
-  locator->Delete();
 }
 
 //----------------------------------------------------------------------------
diff --git a/Imaging/Stencil/vtkPolyDataToImageStencil.h b/Imaging/Stencil/vtkPolyDataToImageStencil.h
index 1b79281..751130f 100644
--- a/Imaging/Stencil/vtkPolyDataToImageStencil.h
+++ b/Imaging/Stencil/vtkPolyDataToImageStencil.h
@@ -94,11 +94,10 @@ protected:
                        int extent[6], int threadId);
 
   static void PolyDataCutter(vtkPolyData *input, vtkPolyData *output,
-                             double z, vtkMergePoints *locator);
+                             double z);
 
   static void PolyDataSelector(vtkPolyData *input, vtkPolyData *output,
-                               double z, double thickness,
-                               vtkMergePoints *locator);
+                               double z, double thickness);
 
   virtual int RequestData(vtkInformation *, vtkInformationVector **,
                           vtkInformationVector *);