diff --git a/Code/Source/sv/Geometry/sv_sys_geom.cxx b/Code/Source/sv/Geometry/sv_sys_geom.cxx
index 5a720ebe0..a7569345a 100644
--- a/Code/Source/sv/Geometry/sv_sys_geom.cxx
+++ b/Code/Source/sv/Geometry/sv_sys_geom.cxx
@@ -4126,19 +4126,23 @@ int sys_geom_set_ids_for_caps( cvPolyData *pd,cvPolyData **outpd,int **doublecap
   for (int i = 0; i < facemax; i++)
   {
     double facecaprange[2];
-    vtkNew(vtkThreshold,threshold1);
+    auto threshold_surface = VtkUtils_ThresholdSurface(i+1.0, i+1.0, "ModelFaceID", geom);
+    /* dp
+    vtkNew(vtkThreshold, threshold1);
     threshold1->SetInputData(geom);
     threshold1->SetInputArrayToProcess(0,0,0,1,"ModelFaceID");
-    //dp threshold1->ThresholdBetween(i+1,i+1);
+    threshold1->ThresholdBetween(i+1,i+1);
     threshold1->Update();
 
     vtkNew(vtkDataSetSurfaceFilter,surfacer);
     surfacer->SetInputData(threshold1->GetOutput());
     surfacer->Update();
+    */
 
-    vtkNew(vtkIntArray,modelfacecaps);
+    vtkNew(vtkIntArray, modelfacecaps);
 
-    if (VtkUtils_PDCheckArrayName(surfacer->GetOutput(),1,"CapID") != SV_OK)
+    if (VtkUtils_PDCheckArrayName(threshold_surface, 1, "CapID") != SV_OK)
+    //dp if (VtkUtils_PDCheckArrayName(surfacer->GetOutput(),1,"CapID") != SV_OK)
     {
       fprintf(stderr,"Second\n");
       fprintf(stderr,"CapID Array is not on the surface\n");
@@ -4146,12 +4150,16 @@ int sys_geom_set_ids_for_caps( cvPolyData *pd,cvPolyData **outpd,int **doublecap
       delete [] captwo;
       return SV_ERROR;
     }
-    modelfacecaps = vtkIntArray::SafeDownCast(surfacer->GetOutput()->
-	GetCellData()->GetArray("CapID"));
+
+    modelfacecaps = vtkIntArray::SafeDownCast(threshold_surface->GetCellData()->GetArray("CapID"));
+    //dp modelfacecaps = vtkIntArray::SafeDownCast(surfacer->GetOutput()->GetCellData()->GetArray("CapID"));
+
     capone[i] = 0;
     captwo[i] = 0;
     (*doublecaps)[i] = 0;
-    for (int j = 0; j < surfacer->GetOutput()->GetNumberOfCells();j++)
+
+    for (int j = 0; j < threshold_surface->GetNumberOfCells();j++)
+    //dp for (int j = 0; j < surfacer->GetOutput()->GetNumberOfCells();j++)
     {
       if (modelfacecaps->GetValue(j) == 1)
         capone[i] = 1;
diff --git a/Code/Source/sv/Mesh/VMTKUtils/sv_vmtk_utils.cxx b/Code/Source/sv/Mesh/VMTKUtils/sv_vmtk_utils.cxx
index fb368c784..f61ee9c30 100644
--- a/Code/Source/sv/Mesh/VMTKUtils/sv_vmtk_utils.cxx
+++ b/Code/Source/sv/Mesh/VMTKUtils/sv_vmtk_utils.cxx
@@ -40,6 +40,7 @@
  */
 
 #include "sv_vmtk_utils.h"
+#include "sv_vtk_utils.h"
 
 #include "SimVascular.h"
 #include "sv_tetgenmesh_utils.h"
@@ -1612,10 +1613,11 @@ int VMTKUtils_CreateBoundaryLayerSurfaceAndCaps(vtkUnstructuredGrid* boundaryMes
 
   // Create boundary layer surface mesh from cells with 'isSurface' array value 1. 
   //
+  /* dp
   auto thresholder = vtkSmartPointer<vtkThreshold>::New();
   thresholder->SetInputData(boundaryMesh);
   thresholder->SetInputArrayToProcess(0, 0, 0, 1, "isSurface");
-  //dp thresholder->ThresholdBetween(1, 1);
+  thresholder->ThresholdBetween(1, 1);
   thresholder->Update();
 
   auto surfacer = vtkSmartPointer<vtkDataSetSurfaceFilter>::New();
@@ -1623,27 +1625,35 @@ int VMTKUtils_CreateBoundaryLayerSurfaceAndCaps(vtkUnstructuredGrid* boundaryMes
   surfacer->Update();
 
   boundaryMeshSurface->DeepCopy(surfacer->GetOutput());
+  */
+  auto threshold_surface = VtkUtils_ThresholdSurface(1.0, 1.0, "isSurface", boundaryMesh);
+  boundaryMeshSurface->DeepCopy(threshold_surface);
 
   // Create boundary layer volume mesh from cells with 'isSurface' array value 0. 
   //
+  /* dp
   thresholder->SetInputData(boundaryMesh);
   thresholder->SetInputArrayToProcess(0, 0, 0, 1, "isSurface");
-  //dp thresholder->ThresholdBetween(0, 0);
+  thresholder->ThresholdBetween(0, 0);
   thresholder->Update();
-
   boundaryMeshVolume->DeepCopy(thresholder->GetOutput());
+  */
+  auto threshold_volume = VtkUtils_ThresholdUgrid(0.0, 0.0, "isSurface", boundaryMesh);
+  boundaryMeshVolume->DeepCopy(threshold_volume);
 
   // Create boundary layer mesh caps from cells with 'WallID' array value 0. 
   //
+  /* dp
   thresholder->SetInputData(surfaceWithSize);
   thresholder->SetInputArrayToProcess(0,0,0,1,"WallID");
-  //dp thresholder->ThresholdBetween(0,0);
+  dp thresholder->ThresholdBetween(0,0);
   thresholder->Update();
-
   surfacer->SetInputData(thresholder->GetOutput());
   surfacer->Update();
-
   surfaceMeshCaps->DeepCopy(surfacer->GetOutput());
+  */
+  auto threshold_caps = VtkUtils_ThresholdSurface(0.0, 0.0, "WallID", surfaceWithSize);
+  surfaceMeshCaps->DeepCopy(threshold_caps);
 
   // Set the values of the 'ModelFaceID' array for the caps to 9999(?). 
   //
diff --git a/Code/Source/sv/Model/OCCTSolidModel/sv_OCCTSolidModel.cxx b/Code/Source/sv/Model/OCCTSolidModel/sv_OCCTSolidModel.cxx
index 780ac868b..63f3310d5 100644
--- a/Code/Source/sv/Model/OCCTSolidModel/sv_OCCTSolidModel.cxx
+++ b/Code/Source/sv/Model/OCCTSolidModel/sv_OCCTSolidModel.cxx
@@ -1707,33 +1707,36 @@ int cvOCCTSolidModel::GetOnlyPD(vtkPolyData *pd,double &max_dist) const
     cleaner->SetInputData(pd);
     cleaner->PointMergingOn();
     cleaner->Update();
+
+    auto threshold_surface = VtkUtils_ThresholdSurface(7.0, 7.0, "MESH_TYPES", cleaner->GetOutput());
+
+    #ifdef dp_old 
     vtkSmartPointer<vtkThreshold> thresholder =
       vtkSmartPointer<vtkThreshold>::New();
     thresholder->SetInputData(cleaner->GetOutput());
+
     //Set Input Array to 0 port,0 connection,1 for Cell Data, and MESh_TYPES is the type name
-    thresholder->SetInputArrayToProcess(0,0,0,1,"MESH_TYPES");
     //Source polydata is on MESH_TYPE 7
-    //dp thresholder->ThresholdBetween(7,7);
+    thresholder->SetInputArrayToProcess(0,0,0,1,"MESH_TYPES");
+    thresholder->ThresholdBetween(7,7);
     thresholder->Update();
     //Extract surface
     vtkSmartPointer<vtkDataSetSurfaceFilter> surfacer =
       vtkSmartPointer<vtkDataSetSurfaceFilter>::New();
     surfacer->SetInputData(thresholder->GetOutput());
     surfacer->Update();
-    //For polydata of just edges, MESH_TYPE is not 7
-    //Only want if not edges
-    //fprintf(stderr,"Num Points bef! %d\n",surfacer->GetOutput()->GetNumberOfPoints());
-    //  vtkSmartPointer<vtkQuadricDecimation> decimator =
-    //    vtkSmartPointer<vtkQuadricDecimation>::New();
-    //  decimator->SetInputData(surfacer->GetOutput());
-    //  decimator->SetTargetReduction(0.2);
-    //  decimator->GetOutput();
-    //  pd->DeepCopy(decimator->GetOutput());
+    #endif
+
+    if (threshold_surface->GetNumberOfPoints() != 0)
+    {
+      pd->DeepCopy(threshold_surface);
+    }
+    /* dp
     if (surfacer->GetOutput()->GetNumberOfPoints() != 0)
     {
       pd->DeepCopy(surfacer->GetOutput());
     }
-    //fprintf(stderr,"Num Points now! %d\n",pd->GetNumberOfPoints());
+    */
   }
 
   return SV_OK;
diff --git a/Code/Source/sv/Utils/sv_vtk_utils.cxx b/Code/Source/sv/Utils/sv_vtk_utils.cxx
index 774cc7300..f7eb59028 100644
--- a/Code/Source/sv/Utils/sv_vtk_utils.cxx
+++ b/Code/Source/sv/Utils/sv_vtk_utils.cxx
@@ -38,13 +38,16 @@
 #include "sv_misc_utils.h"
 #include "sv_cgeom.h"
 
+#include <vtkDataSetSurfaceFilter.h>
+#include <vtkSmartPointer.h>
+#include <vtkThreshold.h>
+#include <vtkUnstructuredGrid.h>
 
 #define START  -1
 #define INTERMED  0
 #define DEADEND   1
 #define CLOSED    2
 
-
 typedef struct {
   double x, y, z;
   int postId;
@@ -78,6 +81,47 @@ static int GetStartPt( int *state, int stateSz );
 static int UpdateCells( int *cells, int numCells,
 			vtkIdType **newCells, int *numNewCells );
 
+//-------------------------
+// VtkUtils_ThresholdUgrid
+//-------------------------
+// Create a vtkUnstructuredGrid object from a threshold of a cellular 
+// data array contained in a vtkDataObject object.
+//
+vtkSmartPointer<vtkUnstructuredGrid>
+VtkUtils_ThresholdUgrid(const double lower, const double upper, const std::string& data_name, 
+    vtkDataObject* vtk_data)
+{
+  int idx = 0;
+  int port = 0;
+  int connection = 0;
+  auto fieldAssociation = vtkDataObject::FieldAssociations::FIELD_ASSOCIATION_CELLS;
+
+  auto threshold = vtkSmartPointer<vtkThreshold>::New();
+  threshold->SetLowerThreshold(lower);
+  threshold->SetUpperThreshold(upper);
+  threshold->SetInputData(vtk_data);
+  threshold->SetInputArrayToProcess(idx, port, connection, fieldAssociation, data_name.c_str());
+  threshold->Update();
+  return threshold->GetOutput();
+}
+
+//---------------------------
+// VtkUtils_ThresholdSurface
+//---------------------------
+// Create a vtkPolyData object from a threshold of a cellular 
+// data array contained in a vtkDataObject object.
+//
+vtkSmartPointer<vtkPolyData>
+VtkUtils_ThresholdSurface(const double lower, const double upper, const std::string& data_name, 
+    vtkDataObject* vtk_data)
+{
+  auto threshold_ugrid = VtkUtils_ThresholdUgrid(lower, upper, data_name, vtk_data);
+
+  auto surfacer = vtkSmartPointer<vtkDataSetSurfaceFilter>::New();
+  surfacer->SetInputData(vtk_data);
+  surfacer->Update();
+  return surfacer->GetOutput();
+}
 
 // -----------------------
 // VtkUtils_NewVtkPolyData
diff --git a/Code/Source/sv/Utils/sv_vtk_utils.h b/Code/Source/sv/Utils/sv_vtk_utils.h
index ab56b4a2f..7d1b9bf73 100644
--- a/Code/Source/sv/Utils/sv_vtk_utils.h
+++ b/Code/Source/sv/Utils/sv_vtk_utils.h
@@ -59,6 +59,14 @@ int SV_EXPORT_UTILS VtkUtils_NewVtkPolyDataLines( vtkPolyData **pd, int numPts,
 // same (potentially redundant) point set, and just change the id's
 // that get used by Lines and Polys.
 
+vtkSmartPointer<vtkUnstructuredGrid>
+VtkUtils_ThresholdUgrid(const double lower, const double upper, const std::string& data_name, 
+    vtkDataObject* vtk_data);
+
+vtkSmartPointer<vtkPolyData>
+VtkUtils_ThresholdSurface(const double lower, const double upper, const std::string& data_name, 
+    vtkDataObject* vtk_data);
+
 int SV_EXPORT_UTILS VtkUtils_FixTopology( vtkPolyData *pd, double tol );
 
 int SV_EXPORT_UTILS VtkUtils_GetPoints( vtkPolyData *pd, double **pts, int *numPts );
diff --git a/Code/Source/sv2/ImageProcessing/sv2_img_threshold.cxx b/Code/Source/sv2/ImageProcessing/sv2_img_threshold.cxx
index 7633aa512..8686b4958 100644
--- a/Code/Source/sv2/ImageProcessing/sv2_img_threshold.cxx
+++ b/Code/Source/sv2/ImageProcessing/sv2_img_threshold.cxx
@@ -33,6 +33,8 @@
 #include "sv_VTK.h"
 #include "vtkThresholdPoints.h"
 
+// [TODO:DaveP] This does not seem to be used anywhere.
+//
 int img_threshold (vtkStructuredPoints *image, vtkFloatingPointType thrMin, vtkFloatingPointType thrMax,
                    int max_num_pts, cvPolyData **result) {
 
diff --git a/Code/Source/sv3/Path/sv3_PathUtils.cxx b/Code/Source/sv3/Path/sv3_PathUtils.cxx
index 0150cb709..fb668e7ec 100644
--- a/Code/Source/sv3/Path/sv3_PathUtils.cxx
+++ b/Code/Source/sv3/Path/sv3_PathUtils.cxx
@@ -31,6 +31,7 @@
 
 #include "SimVascular.h"
 #include "sv3_PathUtils.h"
+#include "sv_vtk_utils.h"
 
 #include <vtkCellData.h>
 #include <vtkDataSetSurfaceFilter.h>
@@ -72,33 +73,40 @@ PathUtils::ExtractCenterlinesSections(vtkSmartPointer<vtkPolyData>& centerlines)
   // Extract sections based on cells with group IDs for each centerline ID.
   //
   for (int cid = minId; cid <= maxId; cid++) {
+      auto centerlinesCidThreshold = VtkUtils_ThresholdSurface(cid, cid, CenterlineIdsArrayName, centerlines);
+      /*dp
       auto threshold = vtkSmartPointer<vtkThreshold>::New();
       threshold->SetInputData(centerlines);
       threshold->SetInputArrayToProcess(0, 0, 0, "vtkDataObject::FIELD_ASSOCIATION_CELLS", CenterlineIdsArrayName.c_str());
-      //dp threshold->ThresholdBetween(cid, cid);
+      threshold->ThresholdBetween(cid, cid);
       threshold->Update();
-
+      auto threshold_mesh = VtkUtils_ThresholdPolyData(CenterlineIdsArrayName, centerlines); 
       auto surfacer = vtkSmartPointer<vtkDataSetSurfaceFilter>::New();
-      surfacer->SetInputData(threshold->GetOutput());
+      surfacer->SetInputData(threshold_mesh);
       surfacer->Update();
       auto centerlinesCidThreshold = surfacer->GetOutput();
+      */
 
       auto groupData = centerlinesCidThreshold->GetCellData()->GetArray(GroupIdsArrayName.c_str());
       double lowerValue = groupData->GetRange()[0];
       double upperValue = groupData->GetRange()[1];
 
+      auto group_threshold = VtkUtils_ThresholdSurface(lowerValue, upperValue, GroupIdsArrayName, centerlinesCidThreshold);
+      /* dp
       auto groupThreshold = vtkSmartPointer<vtkThreshold>::New();
       groupThreshold->SetInputData(centerlinesCidThreshold);
       groupThreshold->SetInputArrayToProcess(0, 0, 0, "vtkDataObject::FIELD_ASSOCIATION_CELLS", GroupIdsArrayName.c_str());
-      //dp groupThreshold->ThresholdBetween(lowerValue, upperValue);
+      groupThreshold->ThresholdBetween(lowerValue, upperValue);
       groupThreshold->Update();
 
       auto groupSurfacer = vtkSmartPointer<vtkDataSetSurfaceFilter>::New();
       groupSurfacer->SetInputData(groupThreshold->GetOutput());
       groupSurfacer->Update();
+      */
 
       auto groupCenterlines = vtkSmartPointer<vtkPolyData>::New();
-      groupCenterlines->DeepCopy(groupSurfacer->GetOutput());
+      groupCenterlines->DeepCopy(group_threshold);
+      //dp groupCenterlines->DeepCopy(groupSurfacer->GetOutput());
       pathsGeometry.push_back(groupCenterlines);
   }
 
diff --git a/Code/Source/sv4gui/Modules/Mesh/Common/sv4gui_MeshLegacyIO.cxx b/Code/Source/sv4gui/Modules/Mesh/Common/sv4gui_MeshLegacyIO.cxx
index 5190fa8df..96ea6427d 100644
--- a/Code/Source/sv4gui/Modules/Mesh/Common/sv4gui_MeshLegacyIO.cxx
+++ b/Code/Source/sv4gui/Modules/Mesh/Common/sv4gui_MeshLegacyIO.cxx
@@ -33,6 +33,7 @@
 #include "sv4gui_MitkMeshIO.h"
 
 #include "sv_polydatasolid_utils.h"
+#include "sv_vtk_utils.h"
 
 #include <QDir>
 
@@ -185,10 +186,15 @@ bool sv4guiMeshLegacyIO::WriteFiles(mitk::DataNode::Pointer meshNode, sv4guiMode
 
       // Extract surface mesh.
       //
+      auto threshold_surface = VtkUtils_ThresholdSurface(i, i, "ModelRegionID", surfaceMesh);
+      if (threshold_surface->GetNumberOfCells() == 0) { 
+        continue;
+      }
+      /* dp
       auto surfThresholder = vtkSmartPointer<vtkThreshold>::New();
       surfThresholder->SetInputData(surfaceMesh);
       surfThresholder->SetInputArrayToProcess(0,0,0,1,"ModelRegionID");
-      //dp surfThresholder->ThresholdBetween(i,i);
+      surfThresholder->ThresholdBetween(i,i);
       surfThresholder->Update();
       if (surfThresholder->GetOutput()->GetNumberOfCells() == 0) {
         continue;
@@ -200,23 +206,30 @@ bool sv4guiMeshLegacyIO::WriteFiles(mitk::DataNode::Pointer meshNode, sv4guiMode
       if (surfacer->GetOutput()->GetNumberOfCells() == 0) {
         continue;
       }
+      */
 
       // Extract volume mesh.
+      auto threshold_volume = VtkUtils_ThresholdUgrid(i, i, "ModelRegionID", volumeMesh);
+      if (threshold_volume->GetNumberOfCells() == 0) {
+        continue;
+      }
+      /* dp
       auto volumeThresholder = vtkSmartPointer<vtkThreshold>::New();
       volumeThresholder->SetInputData(volumeMesh);
       volumeThresholder->SetInputArrayToProcess(0,0,0,1,"ModelRegionID");
-      //dp volumeThresholder->ThresholdBetween(i,i);
+      volumeThresholder->ThresholdBetween(i,i);
       volumeThresholder->Update();
       if (volumeThresholder->GetOutput()->GetNumberOfCells() == 0) {
         continue;
       }
+      */
 
       // Write meshes.
       //
       fprintf(stdout, "[sv4guiMeshLegacyIO::WriteFiles] Writing domain %d\n", i);
       QString newDir = meshDir+"_domain-" + QString::number(i);
       QDir().mkpath(newDir);
-      if (WriteFiles(surfacer->GetOutput(), volumeThresholder->GetOutput(), modelElement, newDir) == false) {
+      if (WriteFiles(threshold_surface, threshold_volume, modelElement, newDir) == false) {
         return false;
       }
     }
@@ -322,10 +335,14 @@ bool sv4guiMeshLegacyIO::WriteFiles(vtkSmartPointer<vtkPolyData> surfaceMesh, vt
       if (connectFilter->GetNumberOfExtractedRegions() > 1) {
 
         for (int j = 0; j < connectFilter->GetNumberOfExtractedRegions(); j++) {
+          auto region_surface = VtkUtils_ThresholdSurface(j, j, "RegionId", connectFilter->GetOutput());
+          region_surface->GetCellData()->RemoveArray("RegionId");
+          region_surface->GetPointData()->RemoveArray("RegionId");
+          /* dp
           auto thresholder = vtkSmartPointer<vtkThreshold>::New();
           thresholder->SetInputData(connectFilter->GetOutput());
           thresholder->SetInputArrayToProcess(0,0,0,1,"RegionId");
-          //dp thresholder->ThresholdBetween(j, j);
+          thresholder->ThresholdBetween(j, j);
           thresholder->Update();
 
           auto surfacer = vtkSmartPointer<vtkDataSetSurfaceFilter>::New();
@@ -334,6 +351,7 @@ bool sv4guiMeshLegacyIO::WriteFiles(vtkSmartPointer<vtkPolyData> surfaceMesh, vt
           surfacer->GetOutput()->GetCellData()->RemoveArray("RegionId");
           surfacer->GetOutput()->GetPointData()->RemoveArray("RegionId");
           auto region_surface = surfacer->GetOutput(); 
+          */
 
           vtpFilePath = meshDir + "/walls_combined_connected_region_" + QString::number(j) + ".vtp";
           vtpFilePath = QDir::toNativeSeparators(vtpFilePath);
diff --git a/Code/Source/sv4gui/Modules/Model/Common/sv4gui_ModelUtils.cxx b/Code/Source/sv4gui/Modules/Model/Common/sv4gui_ModelUtils.cxx
index 328b9e414..77cfee743 100644
--- a/Code/Source/sv4gui/Modules/Model/Common/sv4gui_ModelUtils.cxx
+++ b/Code/Source/sv4gui/Modules/Model/Common/sv4gui_ModelUtils.cxx
@@ -42,6 +42,7 @@
 #include "sv_vmtk_utils.h"
 #include "sv_PolyData.h"
 #include "sv_polydatasolid_utils.h"
+#include "sv_vtk_utils.h"
 
 #include <vtkCellType.h>
 #include <vtkFillHolesFilter.h>
@@ -1288,12 +1289,14 @@ std::vector<sv4guiPathElement::sv4guiPathPoint> sv4guiModelUtils::ConvertToPathP
     return pathPoints;
 }
 
-vtkSmartPointer<vtkPolyData> sv4guiModelUtils::GetThresholdRegion(vtkSmartPointer<vtkPolyData> pd, vtkDataObject::FieldAssociations dataType, std::string arrayName, double minValue, double maxValue )
+vtkSmartPointer<vtkPolyData> sv4guiModelUtils::GetThresholdRegion(vtkSmartPointer<vtkPolyData> pd, 
+    vtkDataObject::FieldAssociations dataType, std::string arrayName, double minValue, double maxValue )
 {
-    vtkSmartPointer<vtkThreshold> thresholder=vtkSmartPointer<vtkThreshold>::New();
+    auto thresholder = vtkSmartPointer<vtkThreshold>::New();
     thresholder->SetInputData(pd);
     thresholder->SetInputArrayToProcess(0, 0, 0, dataType, arrayName.c_str());
-    //dp thresholder->ThresholdBetween(minValue, maxValue);
+    thresholder->SetLowerThreshold(minValue);
+    thresholder->SetUpperThreshold(maxValue);
     thresholder->Update();
 
     vtkSmartPointer<vtkDataSetSurfaceFilter> surfacer=vtkSmartPointer<vtkDataSetSurfaceFilter>::New();
diff --git a/Code/Source/vtkSV/Modules/Boolean/vtkSVLoopBooleanPolyDataFilter.cxx b/Code/Source/vtkSV/Modules/Boolean/vtkSVLoopBooleanPolyDataFilter.cxx
index e4570775c..8f2364ff8 100644
--- a/Code/Source/vtkSV/Modules/Boolean/vtkSVLoopBooleanPolyDataFilter.cxx
+++ b/Code/Source/vtkSV/Modules/Boolean/vtkSVLoopBooleanPolyDataFilter.cxx
@@ -29,6 +29,8 @@
  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  */
 
+#include "sv_vtk_utils.h"
+
 #include "vtkSVLoopBooleanPolyDataFilter.h"
 
 #include "vtkAppendPolyData.h"
@@ -1507,53 +1509,59 @@ void vtkSVLoopBooleanPolyDataFilter::Impl::PerformBoolean(
   //  {
   //  surfaces[i]->Delete();
   //  }
-  vtkSmartPointer<vtkThreshold> thresholder =
-    vtkSmartPointer<vtkThreshold>::New();
-  vtkSmartPointer<vtkDataSetSurfaceFilter> surfacer =
-    vtkSmartPointer<vtkDataSetSurfaceFilter>::New();
 
-  vtkSmartPointer<vtkPolyData> surface1_A =
-	  vtkSmartPointer<vtkPolyData>::New();
+  auto threshold_surface_a = VtkUtils_ThresholdSurface(-1.0, -1.0, "BooleanRegion", this->Mesh[0]); 
+  vtkSmartPointer<vtkPolyData> surface1_A = vtkSmartPointer<vtkPolyData>::New();
+  surface1_A->DeepCopy(threshold_surface_a);
+  /* dp
+  vtkSmartPointer<vtkThreshold> thresholder = vtkSmartPointer<vtkThreshold>::New();
   thresholder->SetInputData(this->Mesh[0]);
   thresholder->SetInputArrayToProcess(0,0,0,1,"BooleanRegion");
-  //dp thresholder->ThresholdBetween(-1,-1);
+  thresholder->ThresholdBetween(-1,-1);
   thresholder->Update();
+
+  vtkSmartPointer<vtkDataSetSurfaceFilter> surfacer = vtkSmartPointer<vtkDataSetSurfaceFilter>::New();
   surfacer->SetInputData(thresholder->GetOutput());
   surfacer->Update();
-  surface1_A->DeepCopy(surfacer->GetOutput());
+  */
 
-  vtkSmartPointer<vtkPolyData> surface1_B =
-	  vtkSmartPointer<vtkPolyData>::New();
+  auto threshold_surface_b = VtkUtils_ThresholdSurface(1.0, 1.0, "BooleanRegion", this->Mesh[0]); 
+  auto surface1_B = vtkSmartPointer<vtkPolyData>::New();
+  surface1_B->DeepCopy(threshold_surface_b);
+   /* dp
   thresholder->SetInputData(this->Mesh[0]);
   thresholder->SetInputArrayToProcess(0,0,0,1,"BooleanRegion");
-  //dp thresholder->ThresholdBetween(1,1);
+  thresholder->ThresholdBetween(1,1);
   thresholder->Update();
   surfacer->SetInputData(thresholder->GetOutput());
   surfacer->Update();
-  surface1_B->DeepCopy(surfacer->GetOutput());
+  */
 
-  vtkSmartPointer<vtkPolyData> surface2_A =
-	  vtkSmartPointer<vtkPolyData>::New();
+  auto threshold_surface2_a = VtkUtils_ThresholdSurface(1.0, 1.0, "BooleanRegion", this->Mesh[1]); 
+  auto surface2_A = vtkSmartPointer<vtkPolyData>::New();
+  surface2_A->DeepCopy(threshold_surface2_a);
+  /* dp
   thresholder->SetInputData(this->Mesh[1]);
   thresholder->SetInputArrayToProcess(0,0,0,1,"BooleanRegion");
-  //dp thresholder->ThresholdBetween(1,1);
+  thresholder->ThresholdBetween(1,1);
   thresholder->Update();
   surfacer->SetInputData(thresholder->GetOutput());
   surfacer->Update();
-  surface2_A->DeepCopy(surfacer->GetOutput());
+  */
 
-  vtkSmartPointer<vtkPolyData> surface2_B =
-	  vtkSmartPointer<vtkPolyData>::New();
+  auto threshold_surface2_b = VtkUtils_ThresholdSurface(-1.0, -1.0, "BooleanRegion", this->Mesh[1]); 
+  auto surface2_B = vtkSmartPointer<vtkPolyData>::New();
+  surface2_B->DeepCopy(threshold_surface2_b);
+  /* dp
   thresholder->SetInputData(this->Mesh[1]);
   thresholder->SetInputArrayToProcess(0,0,0,1,"BooleanRegion");
-  //dp thresholder->ThresholdBetween(-1,-1);
+  thresholder->ThresholdBetween(-1,-1);
   thresholder->Update();
   surfacer->SetInputData(thresholder->GetOutput());
   surfacer->Update();
-  surface2_B->DeepCopy(surfacer->GetOutput());
+  */
 
-  vtkSmartPointer<vtkAppendPolyData> appender =
-    vtkSmartPointer<vtkAppendPolyData>::New();
+  auto appender = vtkSmartPointer<vtkAppendPolyData>::New();
 
   if (booleanOperation == 0)
   {
diff --git a/Code/Source/vtkSV/Modules/Misc/vtkSVGetSphereRegions.cxx b/Code/Source/vtkSV/Modules/Misc/vtkSVGetSphereRegions.cxx
index a9af57a62..c03228bb9 100644
--- a/Code/Source/vtkSV/Modules/Misc/vtkSVGetSphereRegions.cxx
+++ b/Code/Source/vtkSV/Modules/Misc/vtkSVGetSphereRegions.cxx
@@ -380,11 +380,15 @@ int vtkSVGetSphereRegions::SetSphereRegions(vtkPolyData *pd, vtkPolyData *lines,
   radius = new double[numLoops];
 
   // Set up loop thresholder
+
   thresholder->SetInputData(lines);
   thresholder->SetInputArrayToProcess(0,0,0,1,"LoopId");
+
   for (int i=0;i<numLoops;i++)
   {
     // Threshold the loop
+    thresholder->SetLowerThreshold(i);
+    thresholder->SetUpperThreshold(i);
     //dp thresholder->ThresholdBetween(i,i);
     thresholder->Update();
 
diff --git a/Code/Source/vtkSV/Modules/Segmentation/vtkSVCenterlineGraph.cxx b/Code/Source/vtkSV/Modules/Segmentation/vtkSVCenterlineGraph.cxx
index 3cc5256bd..7b5b03fce 100644
--- a/Code/Source/vtkSV/Modules/Segmentation/vtkSVCenterlineGraph.cxx
+++ b/Code/Source/vtkSV/Modules/Segmentation/vtkSVCenterlineGraph.cxx
@@ -29,6 +29,8 @@
  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  */
 
+#include "sv_vtk_utils.h"
+
 #include "vtkSVCenterlineGraph.h"
 
 #include "vtkCellArray.h"
@@ -397,26 +399,31 @@ int vtkSVCenterlineGraph::GetConnectingLineGroups(const int groupId, std::vector
  */
 int vtkSVCenterlineGraph::ComputeGlobalReferenceVectors(vtkSVCenterlineGCell *parent)
 {
+  auto threshold = VtkUtils_ThresholdUgrid(parent->GroupId, parent->GroupId,  this->GroupIdsArrayName, this->Lines); 
+  /* dp
   vtkNew(vtkThreshold, thresholder);
   thresholder->SetInputData(this->Lines);
   thresholder->SetInputArrayToProcess(0, 0, 0, 1, this->GroupIdsArrayName);
   //dp thresholder->ThresholdBetween(parent->GroupId, parent->GroupId);
   thresholder->Update();
   int numPts = thresholder->GetOutput()->GetNumberOfPoints();
+  */
+  int numPts = threshold->GetNumberOfPoints();
 
   double endPt0[3], endPt1[3];
   if (parent->Parent == nullptr && parent->Children.size() == 0)
   {
-    thresholder->GetOutput()->GetPoint(1, endPt0);
-    thresholder->GetOutput()->GetPoint(0, endPt1);
+    threshold->GetPoint(1, endPt0);
+    threshold->GetPoint(0, endPt1);
   }
   else
   {
-    thresholder->GetOutput()->GetPoint(numPts-2, endPt0);
-    thresholder->GetOutput()->GetPoint(numPts-1, endPt1);
+    threshold->GetPoint(numPts-2, endPt0);
+    threshold->GetPoint(numPts-1, endPt1);
   }
-  fprintf(stdout,"End pt 0 is: %.4f %.4f %.4f\n", endPt0[0], endPt0[1], endPt0[2]);
-  fprintf(stdout,"End pt 1 is: %.4f %.4f %.4f\n", endPt1[0], endPt1[1], endPt1[2]);
+
+  //fprintf(stdout,"End pt 0 is: %.4f %.4f %.4f\n", endPt0[0], endPt0[1], endPt0[2]);
+  //fprintf(stdout,"End pt 1 is: %.4f %.4f %.4f\n", endPt1[0], endPt1[1], endPt1[2]);
   vtkMath::Subtract(endPt0, endPt1, this->ReferenceVecs[0]);
   vtkMath::Normalize(this->ReferenceVecs[0]);
 
@@ -425,14 +432,21 @@ int vtkSVCenterlineGraph::ComputeGlobalReferenceVectors(vtkSVCenterlineGCell *pa
   {
     for (int i=0; i<numChildren; i++)
     {
+      double lower = parent->Children[i]->GroupId; 
+      double upper = lower;
+      auto threshold = VtkUtils_ThresholdSurface(lower, upper, this->GroupIdsArrayName, this->Lines); 
+      /*dp
       thresholder->SetInputData(this->Lines);
       thresholder->SetInputArrayToProcess(0, 0, 0, 1, this->GroupIdsArrayName);
-      //dp thresholder->ThresholdBetween(parent->Children[i]->GroupId, parent->Children[i]->GroupId);
+      thresholder->ThresholdBetween(parent->Children[i]->GroupId, parent->Children[i]->GroupId);
       thresholder->Update();
+      */
 
       double startPt[3], secondPt[3];
-      thresholder->GetOutput()->GetPoint(0, startPt);
-      thresholder->GetOutput()->GetPoint(1, secondPt);
+      threshold->GetPoint(0, startPt);
+      threshold->GetPoint(1, secondPt);
+      //dp thresholder->GetOutput()->GetPoint(0, startPt);
+      //dp thresholder->GetOutput()->GetPoint(1, secondPt);
 
       vtkMath::Subtract(secondPt, startPt, parent->Children[i]->BranchVec);
       vtkMath::Normalize(parent->Children[i]->BranchVec);
@@ -507,16 +521,21 @@ int vtkSVCenterlineGraph::ComputeBranchReferenceVectors(vtkSVCenterlineGCell *pa
 {
   //fprintf(stdout,"Child %d of parent %d, dir: %d\n", parent->Children[0]->GroupId, parent->GroupId, parent->Children[0]->BranchDir);
   //fprintf(stdout,"Child %d of parent %d, dir: %d\n", parent->Children[1]->GroupId, parent->GroupId, parent->Children[1]->BranchDir);
+
+  auto threshold = VtkUtils_ThresholdSurface(parent->GroupId, parent->GroupId,  this->GroupIdsArrayName, this->Lines); 
+  /*dp
   vtkNew(vtkThreshold, thresholder);
   thresholder->SetInputData(this->Lines);
   thresholder->SetInputArrayToProcess(0, 0, 0, 1, this->GroupIdsArrayName);
-  //dp thresholder->ThresholdBetween(parent->GroupId, parent->GroupId);
+  thresholder->ThresholdBetween(parent->GroupId, parent->GroupId);
   thresholder->Update();
   int numPts = thresholder->GetOutput()->GetNumberOfPoints();
+  */
+  int numPts = threshold->GetNumberOfPoints();
 
   double endPt0[3], endPt1[3];
-  thresholder->GetOutput()->GetPoint(numPts-2, endPt0);
-  thresholder->GetOutput()->GetPoint(numPts-1, endPt1);
+  threshold->GetPoint(numPts-2, endPt0);
+  threshold->GetPoint(numPts-1, endPt1);
 
   double vec0[3];
   vtkMath::Subtract(endPt0, endPt1, vec0);
@@ -529,14 +548,19 @@ int vtkSVCenterlineGraph::ComputeBranchReferenceVectors(vtkSVCenterlineGCell *pa
   int numChildren = parent->Children.size();
   for (int i=0; i<numChildren; i++)
   {
+    double lower = parent->Children[i]->GroupId; 
+    double upper = lower;
+    auto threshold = VtkUtils_ThresholdSurface(lower, upper, this->GroupIdsArrayName, this->Lines); 
+    /* dp
     thresholder->SetInputData(this->Lines);
     thresholder->SetInputArrayToProcess(0, 0, 0, 1, this->GroupIdsArrayName);
-    //dp thresholder->ThresholdBetween(parent->Children[i]->GroupId, parent->Children[i]->GroupId);
+    thresholder->ThresholdBetween(parent->Children[i]->GroupId, parent->Children[i]->GroupId);
     thresholder->Update();
+    */
 
     double startPt[3], secondPt[3];
-    thresholder->GetOutput()->GetPoint(0, startPt);
-    thresholder->GetOutput()->GetPoint(1, secondPt);
+    threshold->GetPoint(0, startPt);
+    threshold->GetPoint(1, secondPt);
 
     vtkMath::Subtract(secondPt, startPt, parent->Children[i]->BranchVec);
     vtkMath::Normalize(parent->Children[i]->BranchVec);
diff --git a/Code/Source/vtkSV/Modules/Segmentation/vtkSVCenterlines.cxx b/Code/Source/vtkSV/Modules/Segmentation/vtkSVCenterlines.cxx
index d98255d8b..8d7a1b21a 100644
--- a/Code/Source/vtkSV/Modules/Segmentation/vtkSVCenterlines.cxx
+++ b/Code/Source/vtkSV/Modules/Segmentation/vtkSVCenterlines.cxx
@@ -29,6 +29,8 @@
  * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  */
 
+#include "sv_vtk_utils.h"
+
 #include "vtkSVCenterlines.h"
 
 #include "vtkArrayCalculator.h"
@@ -454,32 +456,40 @@ int vtkSVCenterlines::RequestData(
   int mAbsThr = this->AbsoluteThreshold; // TODO: FIGURE OUT GOOD VALUE FOR THIS!!!
   double mAbsRange[2];
   newEdgePd->GetCellData()->GetArray("MAbs")->GetRange(mAbsRange);
+  auto mAbs_threshold = VtkUtils_ThresholdUgrid(mAbsThr, mAbsRange[1], "MAbs", newEdgePd);
+  /*dp 
   vtkNew(vtkThreshold, mAbsThresholder);
   mAbsThresholder->SetInputData(newEdgePd);
   mAbsThresholder->SetInputArrayToProcess(0, 0, 0, 1, "MAbs");
-  //dp mAbsThresholder->ThresholdBetween(mAbsThr, mAbsRange[1]);
+  mAbsThresholder->ThresholdBetween(mAbsThr, mAbsRange[1]);
   mAbsThresholder->Update();
   vtkDebugMacro("Thresholded MAbs: " << mAbsThresholder->GetOutput()->GetNumberOfCells());
-  // ------------------------------------------------------------------------
+  */
 
   // ------------------------------------------------------------------------
   // Threshold based on relative retention
   double mRelThr = this->RelativeThreshold;
   double mRelRange[2];
   newEdgePd->GetCellData()->GetArray("MRel")->GetRange(mRelRange);
+
+  auto MRel_threshold = VtkUtils_ThresholdUgrid(mAbsThr, mAbsRange[1], "MRel", mAbs_threshold);
+  /*dp
   vtkNew(vtkThreshold, mRelThresholder);
-  mRelThresholder->SetInputData(mAbsThresholder->GetOutput());
+  mRelThresholder->SetInputData(mAbs_threshold);
   mRelThresholder->SetInputArrayToProcess(0, 0, 0, 1, "MRel");
-  //dp mRelThresholder->ThresholdBetween(mRelThr, mRelRange[1]);
+  mRelThresholder->ThresholdBetween(mRelThr, mRelRange[1]);
   mRelThresholder->Update();
   vtkDebugMacro("Thresholded MRel: " << mRelThresholder->GetOutput()->GetNumberOfCells());
+  */
+
   // ------------------------------------------------------------------------
 
   // ------------------------------------------------------------------------
   // Get the medial edges
   // Get all separated regions from the thresholded polydata
   vtkNew(vtkConnectivityFilter, connector);
-  connector->SetInputData(mRelThresholder->GetOutput());
+  connector->SetInputData(MRel_threshold);
+  //dp connector->SetInputData(mRelThresholder->GetOutput());
   connector->SetExtractionModeToAllRegions();
   connector->ColorRegionsOn();
   connector->Update();