diff --git a/src/vtkh/filters/CMakeLists.txt b/src/vtkh/filters/CMakeLists.txt
index f7f85fc7..15058c64 100644
--- a/src/vtkh/filters/CMakeLists.txt
+++ b/src/vtkh/filters/CMakeLists.txt
@@ -13,6 +13,7 @@ set(vtkh_filters_headers
   GhostStripper.hpp
   HistSampling.hpp
   Histogram.hpp
+  ValGradHistogram.hpp
   Log.hpp
   IsoVolume.hpp
   NoOp.hpp
@@ -54,6 +55,7 @@ set(vtkh_filters_sources
   GhostStripper.cpp
   HistSampling.cpp
   Histogram.cpp
+  ValGradHistogram.cpp
   Log.cpp
   IsoVolume.cpp
   NoOp.cpp
diff --git a/src/vtkh/filters/HistSampling.cpp b/src/vtkh/filters/HistSampling.cpp
index 86590709..b9483a93 100644
--- a/src/vtkh/filters/HistSampling.cpp
+++ b/src/vtkh/filters/HistSampling.cpp
@@ -2,8 +2,13 @@
 #include <vtkh/filters/GhostStripper.hpp>
 #include <vtkh/filters/Threshold.hpp>
 #include <vtkh/filters/Histogram.hpp>
+#include <vtkh/filters/ValGradHistogram.hpp>
 #include <vtkh/Error.hpp>
 
+#include <vtkh/filters/Gradient.hpp>
+#include <vtkh/filters/VectorMagnitude.hpp>
+
+
 #include <vtkm/worklet/DispatcherMapField.h>
 #include <vtkm/worklet/WorkletMapField.h>
 
@@ -15,7 +20,15 @@
 #include <iostream>
 #include <algorithm>
 #include <vtkm/worklet/WorkletMapField.h>
-#include <iostream>
+
+#include <vtkh/filters/CleanGrid.hpp>
+
+#include <vtkm/worklet/Keys.h>
+#include <vtkm/worklet/ExtractPoints.h>
+#include <vtkm/filter/CleanGrid.h>
+#include <vtkh/filters/Recenter.hpp>
+
+#include <vtkh/vtkm_filters/vtkmGradient.hpp>
 
 namespace vtkh
 {
@@ -57,7 +70,7 @@ class RandomGenerate : public vtkm::worklet::WorkletMapField
 
 
 vtkm::cont::ArrayHandle<vtkm::Float32>
-calculate_pdf(const vtkm::Int32 tot_points,
+calculate_1d_pdf(const vtkm::Int32 tot_points,
               const vtkm::Int32 num_bins,
               const vtkm::Float32 sample_percent,
               vtkm::cont::ArrayHandle<vtkm::Id> mybins)
@@ -130,12 +143,112 @@ calculate_pdf(const vtkm::Int32 tot_points,
   return acceptanceProbsVec;
 }
 
+vtkm::cont::ArrayHandle<vtkm::Float32>
+calculate_2d_pdf(const vtkm::Int32 tot_points,
+              const vtkm::Float32 sample_percent,
+              ValGradHistogram::NonSparseHistogramResult bivar_histogram)
+{
+
+  vtkm::Int32 num_bins = bivar_histogram.m_nob;
+
+  Histogram::HistogramResult valMarginal = bivar_histogram.getValDist();
+
+  vtkm::cont::ArrayHandleIndex indexArray(num_bins);
+  vtkm::cont::ArrayHandle<vtkm::Id> indices;
+  vtkm::cont::Algorithm::Copy(indexArray,indices);
+
+  vtkm::cont::ArrayHandle<vtkm::Id> val_counts_copy;
+  vtkm::cont::Algorithm::Copy(valMarginal.m_bins,val_counts_copy);
+
+  vtkm::cont::ArrayHandleZip <vtkm::cont:: ArrayHandle <vtkm::Id >, 
+                              vtkm::cont:: ArrayHandle <vtkm::Id >> 
+                                zipArray(val_counts_copy , indices );
+
+  //sorts the histogram bins based on their counts
+  vtkm::cont:: Algorithm ::Sort(zipArray );
+  vtkm::cont::ArrayHandleZip <vtkm::cont::ArrayHandle <vtkm::Id >, 
+                              vtkm::cont:: ArrayHandle <vtkm::Id >>
+                                ::PortalConstControl binPortal = zipArray.GetPortalConstControl();
+
+  
+
+  vtkm::Float32 remainingSamples = sample_percent*tot_points;
+  vtkm::Float32 remainingBins = num_bins;
+  std::vector<vtkm::Float32> targetSamples(num_bins, 0.0);
+
+  for (int i = 0; i < num_bins; ++i)
+  {
+    vtkm::Float32 targetNeededSamples = remainingSamples/(1.0*remainingBins);
+    vtkm::Float32 curCount = (vtkm::Float32)binPortal.Get(i).first;
+    vtkm::Float32 samplesTaken;
+    
+    if (curCount<targetNeededSamples)
+    {
+      samplesTaken = curCount;
+    }
+    else
+    {
+      samplesTaken = targetNeededSamples;
+    }
+    targetSamples[binPortal.Get(i).second] = samplesTaken;
+    remainingBins = remainingBins-1;
+    remainingSamples = remainingSamples - samplesTaken; 
+  }
+
+
+  // Start calculating the 2d acceptance probability histogram (flattened)
+  vtkm::cont::ArrayHandle<vtkm::Float32> acceptanceProbsVec;
+  acceptanceProbsVec.Allocate(num_bins*num_bins);
+  auto acceptance_portal = acceptanceProbsVec.WritePortal();
+  for(int i = 0; i < num_bins*num_bins; ++i)
+  {
+    acceptance_portal.Set(i, -1.f);
+  }
+  
+
+  // for each bin in the val histogram, identify how many samples to pick from each 
+  // bin of the corresponding marginal gradient histogram bin
+
+  vtkm::cont::ArrayHandle<vtkm::Id> biVarFreq;
+  vtkm::cont::Algorithm::Copy(bivar_histogram.m_freqs,biVarFreq);
+  auto biVarFreqPortal = biVarFreq.ReadPortal();
+
+  for(int i=0; i<num_bins; i++)
+  {
+    vtkm::Float32 current_bin_target_samples = targetSamples[i];
+
+    for(int j=num_bins-1; j>=0 ; j--)
+    {
+      vtkm::Id curGradFreq = biVarFreqPortal.Get(i*num_bins + j);
+      if(current_bin_target_samples <= curGradFreq)
+      {
+        acceptance_portal.Set(i*num_bins + j, current_bin_target_samples);
+        current_bin_target_samples = 0;
+      }
+      else{
+        acceptance_portal.Set(i*num_bins + j, curGradFreq);
+        current_bin_target_samples -= curGradFreq;
+      }
+      if(curGradFreq > 0)
+      {
+        vtkm::Float32 temp_val = (acceptance_portal.Get(i*num_bins + j)*1.0)/curGradFreq;
+        acceptance_portal.Set(i*num_bins + j, temp_val);
+      }
+      else
+        acceptance_portal.Set(i*num_bins + j, 0.0);
+    }
+  }
+
+  return acceptanceProbsVec;
+}
+
 
 }
 
 HistSampling::HistSampling()
   : m_sample_percent(0.1f),
-    m_num_bins(128)
+    m_num_bins(128),
+    m_use_gradient(false)
 {
 
 }
@@ -165,6 +278,12 @@ HistSampling::SetNumBins(const int num_bins)
   m_num_bins = num_bins;
 }
 
+void
+HistSampling::SetGradientSampling(const bool use_gradient)
+{
+  m_use_gradient = use_gradient;
+}
+
 void
 HistSampling::SetField(const std::string &field_name)
 {
@@ -194,14 +313,14 @@ HistSampling::GetField() const
   return m_field_name;
 }
 
-struct LookupWorklet : public vtkm::worklet::WorkletMapField
+struct Lookup_1d_Worklet : public vtkm::worklet::WorkletMapField
 {
 protected:
   vtkm::Id m_num_bins;
   vtkm::Float64 m_min;
   vtkm::Float64 m_bin_delta;
 public:
-  LookupWorklet(const vtkm::Id num_bins,
+  Lookup_1d_Worklet(const vtkm::Id num_bins,
                 const vtkm::Float64 min_value,
                 const vtkm::Float64 bin_delta)
     : m_num_bins(num_bins),
@@ -231,6 +350,186 @@ struct LookupWorklet : public vtkm::worklet::WorkletMapField
   }
 };
 
+struct LookupID_worklet : public vtkm::worklet::WorkletMapField
+{
+public:
+  using ControlSignature = void(FieldIn, FieldOut);
+  using ExecutionSignature = _2(_1, InputIndex);
+
+  
+  VTKM_EXEC vtkm::Id operator()(const vtkm::Id &field_value,
+                                   vtkm::Id inIndex) const
+  {
+    vtkm::Id boolVal = static_cast<vtkm::Id>(field_value);
+    
+
+    if(boolVal)
+      return inIndex;
+    else
+      return 0;
+  }
+};
+
+struct SamplingWorklet_1d : public vtkm::worklet::WorkletMapField
+{
+protected:
+  vtkm::Id m_num_bins;
+  vtkm::Float64 m_min;
+  vtkm::Float64 m_bin_delta;
+public:
+  SamplingWorklet_1d(const vtkm::Id num_bins,
+                const vtkm::Float64 min_value,
+                const vtkm::Float64 bin_delta)
+    : m_num_bins(num_bins),
+      m_min(min_value),
+      m_bin_delta(bin_delta)
+  {}
+
+  using ControlSignature = void(FieldIn, FieldOut, WholeArrayIn, FieldIn);
+  using ExecutionSignature = _2(_1, _3, _4, InputIndex);
+
+  template <typename TablePortal>
+  VTKM_EXEC vtkm::Id operator()(const vtkm::Float64 &field_value,
+                                   TablePortal table,
+                                   const vtkm::Float32 &random,
+                                   vtkm::Id inIndex) const
+  {
+    vtkm::Id bin = static_cast<vtkm::Id>((field_value - m_min) / m_bin_delta);
+    if(bin < 0)
+    {
+      bin = 0;
+    }
+    if(bin >= m_num_bins)
+    {
+      bin = m_num_bins - 1;
+    }
+
+    if(random < table.Get(bin))
+      return inIndex;
+    else
+      return 0;
+  }
+};
+
+struct Lookup_2d_Worklet : public vtkm::worklet::WorkletMapField
+{
+protected:
+  vtkm::Id m_num_bins;
+  vtkm::Float64 m_min_val;
+  vtkm::Float64 m_bin_delta_val;
+  vtkm::Float64 m_min_grad;
+  vtkm::Float64 m_bin_delta_grad;
+public:
+  Lookup_2d_Worklet(const vtkm::Id num_bins,
+                const vtkm::Float64 min_val,
+                const vtkm::Float64 delta_val,
+                const vtkm::Float64 min_grad,
+                const vtkm::Float64 delta_grad)
+    : m_num_bins(num_bins),
+      m_min_val(min_val),
+      m_bin_delta_val(delta_val),
+      m_min_grad(min_grad),
+      m_bin_delta_grad(delta_grad)
+  {}
+
+  using ControlSignature = void(FieldIn, FieldIn, FieldOut, WholeArrayIn, FieldIn);
+  using ExecutionSignature = _3(_1, _2, _4, _5);
+
+  template <typename TablePortal>
+  VTKM_EXEC vtkm::UInt8 operator()(const vtkm::Float64 &field_val,
+                                   const vtkm::Float64 &field_grad,
+                                   TablePortal table,
+                                   const vtkm::Float32 &random) const
+  {
+    vtkm::Id bin_val = static_cast<vtkm::Id>((field_val - m_min_val) / m_bin_delta_val);
+    if(bin_val < 0)
+    {
+      bin_val = 0;
+    }
+    if(bin_val >= m_num_bins)
+    {
+      bin_val = m_num_bins - 1;
+    }
+
+    vtkm::Id bin_grad = static_cast<vtkm::Id>((field_grad - m_min_grad) / m_bin_delta_grad);
+    if(bin_grad < 0)
+    {
+      bin_grad = 0;
+    }
+    if(bin_grad >= m_num_bins)
+    {
+      bin_grad = m_num_bins - 1;
+    }
+
+    //flattened bin index
+    vtkm::Id bin = bin_val*m_num_bins + bin_grad;
+
+    return random < table.Get(bin);
+  }
+};
+
+struct SamplingWorklet_2d : public vtkm::worklet::WorkletMapField
+{
+protected:
+  vtkm::Id m_num_bins;
+  vtkm::Float64 m_min_val;
+  vtkm::Float64 m_bin_delta_val;
+  vtkm::Float64 m_min_grad;
+  vtkm::Float64 m_bin_delta_grad;
+public:
+  SamplingWorklet_2d(const vtkm::Id num_bins,
+                const vtkm::Float64 min_val,
+                const vtkm::Float64 delta_val,
+                const vtkm::Float64 min_grad,
+                const vtkm::Float64 delta_grad)
+    : m_num_bins(num_bins),
+      m_min_val(min_val),
+      m_bin_delta_val(delta_val),
+      m_min_grad(min_grad),
+      m_bin_delta_grad(delta_grad)
+  {}
+
+  using ControlSignature = void(FieldIn, FieldIn, FieldOut, WholeArrayIn, FieldIn);
+  using ExecutionSignature = _3(_1, _2, _4, _5, InputIndex);
+
+  template <typename TablePortal>
+  VTKM_EXEC vtkm::Id operator()(const vtkm::Float64 &field_val,
+                                   const vtkm::Float64 &field_grad,
+                                   TablePortal table,
+                                   const vtkm::Float32 &random,
+                                   vtkm::Id inIndex) const
+  {
+    vtkm::Id bin_val = static_cast<vtkm::Id>((field_val - m_min_val) / m_bin_delta_val);
+    if(bin_val < 0)
+    {
+      bin_val = 0;
+    }
+    if(bin_val >= m_num_bins)
+    {
+      bin_val = m_num_bins - 1;
+    }
+
+    vtkm::Id bin_grad = static_cast<vtkm::Id>((field_grad - m_min_grad) / m_bin_delta_grad);
+    if(bin_grad < 0)
+    {
+      bin_grad = 0;
+    }
+    if(bin_grad >= m_num_bins)
+    {
+      bin_grad = m_num_bins - 1;
+    }
+
+    //flattened bin index
+    vtkm::Id bin = bin_val*m_num_bins + bin_grad;
+
+    if(random < table.Get(bin))
+      return inIndex;
+    else
+      return 0;
+
+  }
+};
+
 
 void PrintStatInfo(vtkm::worklet::FieldStatistics<vtkm::Float64>::StatInfo statinfo)
 {
@@ -254,7 +553,6 @@ void PrintStatInfo(vtkm::worklet::FieldStatistics<vtkm::Float64>::StatInfo stati
 
 void HistSampling::DoExecute()
 {
-
   vtkh::DataSet *input = this->m_input;
   bool has_ghosts = m_ghost_field != "";
 
@@ -270,100 +568,349 @@ void HistSampling::DoExecute()
     input = stripper.GetOutput();
   }
 
+  // make sure we have a node-centered field
+  bool valid_field = false;
+  bool is_cell_assoc = m_input->GetFieldAssociation(m_field_name, valid_field) ==
+                       vtkm::cont::Field::Association::CELL_SET;
+  //before recenter
+  if(valid_field && is_cell_assoc)
+  {
+    std::cout << "****THERE!****\n";
+    Recenter recenter;
+    recenter.SetInput(m_input);
+    recenter.SetField(m_field_name);
+    recenter.SetResultAssoc(vtkm::cont::Field::Association::POINTS);
+    recenter.Update();
+    m_input = recenter.GetOutput();
+  }
+  //new dataset
+  DataSet data_set;
+
+   
+  vtkm::cont::ArrayHandle <vtkm::Float32 > probArray;
+  vtkm::Range val_range, grad_range;
+  vtkm::Float64 val_delta, grad_delta;
   const int num_domains = input->GetNumberOfDomains();
+  vtkh::DataSet *mag_output;
 
-  Histogram histogrammer;
-  histogrammer.SetNumBins(m_num_bins);
-  Histogram::HistogramResult histogram = histogrammer.Run(*input,m_field_name);
-  //histogram.Print(std::cout);
+  vtkh::DataSet *ps_mag_output;
 
-  vtkm::Id numberOfBins = histogram.m_bins.GetNumberOfValues();
 
-  bool valid_field;
-  vtkm::cont::Field::Association assoc = input->GetFieldAssociation(m_field_name,
-                                                                    valid_field);
+  
 
-  for(int i = 0; i < num_domains; ++i)
+  if(!m_use_gradient)
   {
-    vtkm::Range range;
-    vtkm::Float64 delta;
-    vtkm::cont::DataSet &dom = input->GetDomain(i);
+    Histogram histogrammer;
+    histogrammer.SetNumBins(m_num_bins);
 
-    if(!dom.HasField(m_field_name))
+    Histogram::HistogramResult histogram = histogrammer.Run(*input, m_field_name);
+    // histogram.Print(std::cout);
+
+    vtkm::Id global_num_values = histogram.totalCount();
+    vtkm::cont:: ArrayHandle <vtkm::Id > globCounts = histogram.m_bins;
+
+    std::cout << "global_num_values = " << global_num_values << "\n";
+
+    val_range = histogram.m_range;
+    val_delta = histogram.m_bin_delta;
+
+    // calculate the acceptance probability histogram 
+    probArray = detail::calculate_1d_pdf(global_num_values, m_num_bins, m_sample_percent, globCounts);
+  }
+  else
+  {
+    //Calculate the gradient magnitude field
+    vtkh::Gradient grad;
+    grad.SetInput(input);
+    grad.SetField(m_field_name);
+
+    vtkh::GradientParameters params;
+    params.output_name = "grad";
+    params.use_point_gradient = false;
+    grad.SetParameters(params);
+
+    grad.Update();
+
+    vtkh::DataSet *grad_output = grad.GetOutput();
+
+    vtkh::VectorMagnitude mag;
+    mag.SetInput(grad_output);
+    mag.SetField("grad");
+
+    mag.SetResultName("mag");
+    mag.Update();
+
+    mag_output = mag.GetOutput();
+
+    //-----------------------------------------------------
+    //Calculate the gradient magnitude field for point dataset
+    vtkh::Gradient ps_grad;
+    ps_grad.SetInput(input);
+    ps_grad.SetField(m_field_name);
+
+    vtkh::GradientParameters ps_params;
+    ps_params.output_name = "grad";
+    ps_params.use_point_gradient = true;
+    ps_grad.SetParameters(ps_params);
+
+    ps_grad.Update();
+
+    vtkh::DataSet *ps_grad_output = ps_grad.GetOutput();
+
+    vtkh::VectorMagnitude ps_mag;
+    ps_mag.SetInput(ps_grad_output);
+    ps_mag.SetField("grad");
+
+    ps_mag.SetResultName("mag");
+    ps_mag.Update();
+
+    ps_mag_output = ps_mag.GetOutput();
+    std::cout << "ps_mag_output dom = " << ps_mag_output->GetNumberOfDomains() << " \n";
+    //-----------------------------------------------------------
+
+
+    const int grad_num_domains = mag_output->GetNumberOfDomains();
+    //make sure both the fields have same number of domains
+    if( num_domains != grad_num_domains)
     {
-      // We have already check to see if the field exists globally,
-      // so just skip if this particular domain doesn't have the field
-      continue;
+      throw Error("Number of domains does not match for the scalar value field and the gradient magnitude field");
     }
 
-    vtkm::cont::ArrayHandle<vtkm::Float64> data;
-    dom.GetField(m_field_name).GetData().CopyTo(data);
+    ValGradHistogram vg_histogrammer;
+    vg_histogrammer.SetNumBins(m_num_bins);
 
 
-    //vtkm::worklet::FieldStatistics<vtkm::Float64>::StatInfo statinfo;
-    //vtkm::worklet::FieldStatistics<vtkm::Float64>().Run(data, statinfo);
 
-    //std::cout << "Statistics for CELL data:" << std::endl;
-    //PrintStatInfo(statinfo);
+    ValGradHistogram::NonSparseHistogramResult bivar_histogram = vg_histogrammer.Run(*input, m_field_name, *mag_output);
+    // bivar_histogram.Print(std::cout);
 
-    vtkm::cont:: ArrayHandle <vtkm::Id > globCounts = histogram.m_bins;
+    vtkm::Id global_num_values = bivar_histogram.totalCount();
 
-    // start doing sampling
+    std::cout << "global_num_values = " << global_num_values << "\n";
 
-    vtkm::Int32 tot_points = data.GetNumberOfValues();
+    val_range = bivar_histogram.m_range_val;
+    val_delta = bivar_histogram.m_bin_delta_val;
+    grad_range = bivar_histogram.m_range_grad;
+    grad_delta = bivar_histogram.m_bin_delta_grad;
 
-    // use the acceptance probabilities to create a stencil buffer
-    vtkm::cont::ArrayHandle<vtkm::Float32> randArray;
+    // calculate the acceptance probability histogram 
+    probArray = detail::calculate_2d_pdf(global_num_values, m_sample_percent, bivar_histogram);
+  }
 
-    randArray.Allocate(tot_points);
+  //bool valid_field;
+  vtkm::cont::Field::Association assoc = input->GetFieldAssociation(m_field_name, valid_field);
+  
+  for(int i = 0; i < num_domains; ++i)
+  {
+    vtkm::Id domain_id;
+    //vtkm::cont::DataSet &dom = input->GetDomain(i);
+    vtkm::cont::DataSet dom;
+    this->m_input->GetDomain(i, dom, domain_id);
 
-    const vtkm::Int32 seed = 0;
+    if(!dom.HasField(m_field_name))
+    {
+      // We have already check to see if the field exists globally,
+      // so just skip if this particular domain doesn't have the field
+      continue;
+    }
 
-    vtkm::worklet::DispatcherMapField<detail::RandomGenerate>(seed).Invoke(randArray);
+    vtkm::cont::ArrayHandle<vtkm::Float64> val_data;
+    dom.GetField(m_field_name).GetData().CopyTo(val_data);
 
-    vtkm::cont::ArrayHandle <vtkm::Float32 > probArray;
-    probArray = detail::calculate_pdf(tot_points, numberOfBins, m_sample_percent, globCounts);
+    vtkm::Int32 tot_points = val_data.GetNumberOfValues();
 
-    vtkm::cont::ArrayHandle <vtkm::UInt8> stencilBool;
-    vtkm::worklet::DispatcherMapField<LookupWorklet>(LookupWorklet{numberOfBins,
-                                                     histogram.m_range.Min,
-                                                     histogram.m_bin_delta}).Invoke(data,
-                                                                                    stencilBool,
-                                                                                    probArray,
-                                                                                    randArray);
+    // use the acceptance probabilities to create a stencil buffer
+    vtkm::cont::ArrayHandle<vtkm::Float32> randArray;
+    randArray.Allocate(tot_points);
+    const vtkm::Int32 seed = 0;
+    vtkm::worklet::DispatcherMapField<detail::RandomGenerate>(seed).Invoke(randArray);
 
+    vtkm::cont::ArrayHandle <vtkm::Id> stencilBool;
 
     vtkm::cont::ArrayHandle <vtkm::Float32> output;
-    vtkm::cont::Algorithm ::Copy(stencilBool , output );
 
-    vtkm::cont:: DataSetFieldAdd dataSetFieldAdd;
+    vtkm::cont::ArrayHandle <vtkm::Id> stencilIds; 
 
-    if(assoc == vtkm::cont::Field::Association::POINTS)
+    if(!m_use_gradient)
     {
-      dataSetFieldAdd.AddPointField(dom , "valSampled", output );
+      vtkm::worklet::DispatcherMapField<Lookup_1d_Worklet>(Lookup_1d_Worklet{m_num_bins,
+                                                       val_range.Min,
+                                                       val_delta}).Invoke(val_data,
+                                                                          stencilBool,
+                                                                          probArray,
+                                                                          randArray);
+      vtkm::cont::Algorithm ::Copy(stencilBool, output);
+
+      std::cout << "test:::" << "tot_points = " << tot_points << "\n";
+
+      vtkm::worklet::DispatcherMapField<SamplingWorklet_1d>(SamplingWorklet_1d{m_num_bins,
+                                                     val_range.Min,
+                                                     val_delta}).Invoke(val_data,
+                                                                          stencilIds,
+                                                                          probArray,
+                                                                          randArray);
+
+
     }
     else
     {
-      dataSetFieldAdd.AddCellField(dom , "valSampled", output );
+
+      vtkm::cont::DataSet &grad_dom = ps_mag_output->GetDomain(i);
+      
+
+      if(!grad_dom.HasField("mag"))
+      {
+        // We have already check to see if the field exists globally,
+        // so just skip if this particular domain doesn't have the field
+        continue;
+      }
+
+      vtkm::cont::ArrayHandle<vtkm::Float64> grad_data;
+      grad_dom.GetField("mag").GetData().CopyTo(grad_data);
+
+      if(tot_points != grad_data.GetNumberOfValues())
+      {
+        //Sanity check: the number of points in the scalar value field 
+        // and the gradient magnitude field must match.
+        std::cout << "test:::" << "tot_points = " << tot_points << "  grad_numval = " << grad_data.GetNumberOfValues() << "\n";
+        throw Error("Mismatch in the number of values for scalar values and gradient magnitude");
+      }
+
+
+      
+
+      
+      vtkm::worklet::DispatcherMapField<Lookup_2d_Worklet>(Lookup_2d_Worklet{m_num_bins,
+                                                         val_range.Min,
+                                                         val_delta,
+                                                         grad_range.Min,
+                                                         grad_delta}).Invoke(val_data,
+                                                                            grad_data,
+                                                                            stencilBool,
+                                                                            probArray,
+                                                                            randArray);
+
+      vtkm::cont::Algorithm ::Copy(stencilBool, output);
+
+      vtkm::worklet::DispatcherMapField<SamplingWorklet_2d>(SamplingWorklet_2d{m_num_bins,
+                                                         val_range.Min,
+                                                         val_delta,
+                                                         grad_range.Min,
+                                                         grad_delta}).Invoke(val_data,
+                                                                            grad_data,
+                                                                            stencilIds,
+                                                                            probArray,
+                                                                            randArray);
+    }
+
+    
+    // vtkm::worklet::DispatcherMapField<LookupID_worklet>(LookupID_worklet.Invoke(stencilBool,
+    //                                                                       stencilIds);
+    
+
+    // // Test code: Verify the number of eventually selected samples by summing the stencilBool
+    vtkm::Float32 test_sum = 0.0;
+    for(int j=0; j<tot_points; j++)
+    {
+      test_sum += output.ReadPortal().Get(j);
     }
+    std::cout << "Domain::"<< i
+              << ": Samples taken=" << test_sum 
+              << ". Total points=" << tot_points
+              << std::endl;
+
+    // vtkm::cont:: DataSetFieldAdd dataSetFieldAdd;
+
+    // if(assoc == vtkm::cont::Field::Association::POINTS)
+    // {
+    //   dataSetFieldAdd.AddPointField(dom , "valSampled", output );
+    // }
+    // else
+    // {
+    //   dataSetFieldAdd.AddCellField(dom , "valSampled", output );
+    // }
+
+    //////////////Extracting Particle data [Subhashis: new code begin]//////////////////
+    vtkm::cont::ArrayHandle<vtkm::Id> sortedIds;
+    vtkm::cont::ArrayCopy(stencilIds, sortedIds);
+    vtkm::worklet::Keys<vtkm::Id> keys(sortedIds);
+    vtkm::cont::ArrayHandle<vtkm::Id> sampleIds;
+    vtkm::cont::ArrayCopy(keys.GetUniqueKeys(), sampleIds);  //get the uniques keys to have a list of only the sampled indices
+
+    vtkm::Id numUniqueIds = sampleIds.GetNumberOfValues();
+
+    std::cout << "numUniqueIds = " << stencilIds.GetNumberOfValues() << " \n";
+
+    // vtkm::Float32 test_sum1 = 0.0;
+    // for(int j=0; j<stencilIds.GetNumberOfValues(); j++)
+    // {
+    //   test_sum1 = stencilIds.ReadPortal().Get(j);
+    //   std::cout << test_sum1 << "\t ";
+    // }
+    // std::cout << std::endl;
+
+    vtkm::worklet::ExtractPoints extractPoints;
+    vtkm::cont::CellSetSingleType<> outCellSet = extractPoints.Run(dom.GetCellSet(), sampleIds);
+
+    vtkm::cont::DataSet outDataSet;
+    outDataSet.SetCellSet(outCellSet);
+    outDataSet.AddCoordinateSystem(dom.GetCoordinateSystem(0));
+
+    // Associate the relevant field to the extracted cellset
+    vtkm::cont::Field in_field = dom.GetField(m_field_name);
+    outDataSet.AddField(in_field);
+
+    data_set.AddDomain(outDataSet, domain_id);
+
+
   }
 
-  vtkh::Threshold thresher;
-  thresher.SetInput(input);
-  thresher.SetField("valSampled");
+  CleanGrid cleaner;
+  cleaner.SetInput(&data_set);
+  cleaner.Update();
+  this->m_output = cleaner.GetOutput();
+
+  //verification of output dataset
+  vtkm::Id num_output_domain = this->m_output->GetNumberOfDomains();
+  vtkm::Id num_output_cells = this->m_output->GetNumberOfCells();
+
 
-  double upper_bound = 1.;
-  double lower_bound = 1.;
+  std::cout << "outDataSet: num_domains=" << num_output_domain
+              << " num_cells=" << num_output_cells
+              << std::endl;
 
-  thresher.SetUpperThreshold(upper_bound);
-  thresher.SetLowerThreshold(lower_bound);
-  thresher.Update();
-  this->m_output = thresher.GetOutput();
 
+  for(int i=0; i<num_output_domain; i++)
+  {
+    vtkm::cont::DataSet &out_dom = this->m_output->GetDomain(i);
+    std::cout << "Domain " << i << " of the output :";
+    std::cout << " num_coord=" << out_dom.GetNumberOfCoordinateSystems() 
+              << " num_fields=" << out_dom.GetNumberOfFields()
+              << " num_values=" << out_dom.GetField(m_field_name).GetData().GetNumberOfValues()
+              << std::endl << std::endl;
+
+  }
+
+  // vtkh::Threshold thresher;
+  // thresher.SetInput(input);
+  // thresher.SetField("valSampled");
+
+  // double upper_bound = 1.;
+  // double lower_bound = 1.;
+
+  // thresher.SetUpperThreshold(upper_bound);
+  // thresher.SetLowerThreshold(lower_bound);
+  // thresher.Update();
+  // this->m_output = thresher.GetOutput();
+  
   if(has_ghosts)
   {
     delete input;
   }
+  if(m_use_gradient){
+    delete mag_output;
+    delete ps_mag_output;
+  }
 }
 
 std::string
diff --git a/src/vtkh/filters/HistSampling.hpp b/src/vtkh/filters/HistSampling.hpp
index 83b4813d..4920c0a6 100644
--- a/src/vtkh/filters/HistSampling.hpp
+++ b/src/vtkh/filters/HistSampling.hpp
@@ -20,6 +20,7 @@ class VTKH_API HistSampling : public Filter
   void SetGhostField(const std::string &field_name);
   std::string GetField() const;
   void SetSamplingPercent(const float percent);
+  void SetGradientSampling(const bool use_gradient);
 protected:
   void PreExecute() override;
   void PostExecute() override;
@@ -29,6 +30,7 @@ class VTKH_API HistSampling : public Filter
   std::string m_ghost_field;
   float m_sample_percent;
   int m_num_bins;
+  bool m_use_gradient;
 };
 
 } //namespace vtkh
diff --git a/src/vtkh/filters/Histogram.cpp b/src/vtkh/filters/Histogram.cpp
index a3b80254..4d49448d 100644
--- a/src/vtkh/filters/Histogram.cpp
+++ b/src/vtkh/filters/Histogram.cpp
@@ -194,4 +194,17 @@ Histogram::HistogramResult::Print(std::ostream &out)
   out<<"total points: "<<sum<<"\n";
 }
 
+vtkm::Id
+Histogram::HistogramResult::totalCount()
+{
+  auto binPortal = m_bins.ReadPortal();
+  const int num_bins = m_bins.GetNumberOfValues();
+  vtkm::Id sum = 0;
+  for (vtkm::Id i = 0; i < num_bins; i++)
+  {
+    sum += binPortal.Get(i);
+  }
+  return sum;
+}
+
 } //  namespace vtkh
diff --git a/src/vtkh/filters/Histogram.hpp b/src/vtkh/filters/Histogram.hpp
index 70317421..17562cc5 100644
--- a/src/vtkh/filters/Histogram.hpp
+++ b/src/vtkh/filters/Histogram.hpp
@@ -24,6 +24,7 @@ class VTKH_API Histogram
     vtkm::Range m_range;
     vtkm::Float64 m_bin_delta;
     void Print(std::ostream &out);
+    vtkm::Id totalCount();
   };
 
   HistogramResult Run(vtkh::DataSet &data_set, const std::string &field_name);
diff --git a/src/vtkh/filters/ValGradHistogram.cpp b/src/vtkh/filters/ValGradHistogram.cpp
new file mode 100644
index 00000000..14e1266a
--- /dev/null
+++ b/src/vtkh/filters/ValGradHistogram.cpp
@@ -0,0 +1,316 @@
+#include <vtkh/filters/ValGradHistogram.hpp>
+#include <vtkh/filters/Histogram.hpp>
+
+#include <vtkh/Error.hpp>
+#include <vtkh/Logger.hpp>
+#include <vtkh/utils/vtkm_array_utils.hpp>
+
+#include <vtkm/worklet/FieldHistogram.h>
+#include <vtkm/worklet/NDimsHistogram.h>
+
+#ifdef VTKH_PARALLEL
+#include <mpi.h>
+#endif
+
+namespace vtkh
+{
+
+namespace valgrad_detail
+{
+
+ValGradHistogram::NonSparseHistogramResult
+merge_histograms(std::vector<ValGradHistogram::SparseHistogramResult> &histograms)
+{
+  const int size = histograms.size();
+  if(size < 1)
+  {
+    throw Error("2d histogram size is 0");
+  }
+  ValGradHistogram::NonSparseHistogramResult res;
+
+  //initialize the frequencies
+  std::vector<vtkm::Id> freq_vec(histograms[0].m_nob*histograms[0].m_nob, 0);
+
+  for(int h=0; h<size; ++h)
+  {
+    vtkm::Id nonSparseBins = histograms[h].m_binIds[0].ReadPortal().GetNumberOfValues();
+    for (int i = 0; i < nonSparseBins; i++)
+    {
+      vtkm::Id idx0 = histograms[h].m_binIds[0].ReadPortal().Get(i);
+      vtkm::Id idx1 = histograms[h].m_binIds[1].ReadPortal().Get(i);
+      vtkm::Id f = histograms[h].m_freqs.ReadPortal().Get(i);
+
+      freq_vec[idx0*histograms[0].m_nob + idx1] += f;
+    }
+  }
+
+  res.m_freqs = vtkm::cont::make_ArrayHandle(freq_vec, vtkm::CopyFlag::On);
+
+  res.m_nob = histograms[0].m_nob;
+  res.m_range_val = histograms[0].m_range_val;
+  res.m_range_grad = histograms[0].m_range_grad;
+  res.m_bin_delta_val = histograms[0].m_bin_delta_val;
+  res.m_bin_delta_grad = histograms[0].m_bin_delta_grad;  
+
+  return res;
+}
+
+template<typename T>
+void reduce(T *array, int size);
+
+template<>
+void reduce<vtkm::Int32>(vtkm::Int32 *array, int size)
+{
+#ifdef VTKH_PARALLEL
+  MPI_Comm mpi_comm = MPI_Comm_f2c(vtkh::GetMPICommHandle());
+  MPI_Allreduce(MPI_IN_PLACE,array,size, MPI_INT,MPI_SUM,mpi_comm);
+#else
+  (void) array;
+  (void) size;
+#endif
+}
+
+template<>
+void reduce<vtkm::Int64>(vtkm::Int64 *array, int size)
+{
+#ifdef VTKH_PARALLEL
+  MPI_Comm mpi_comm = MPI_Comm_f2c(vtkh::GetMPICommHandle());
+  MPI_Allreduce(MPI_IN_PLACE,array,size, MPI_LONG_LONG,MPI_SUM,mpi_comm);
+#else
+  (void) array;
+  (void) size;
+#endif
+}
+
+
+} // namespace valgrad_detail
+
+ValGradHistogram::ValGradHistogram()
+  : m_num_bins(256)
+{
+
+}
+
+ValGradHistogram::~ValGradHistogram()
+{
+
+}
+
+void
+ValGradHistogram::SetNumBins(const int num_bins)
+{
+  m_num_bins = num_bins;
+}
+
+ValGradHistogram::NonSparseHistogramResult
+ValGradHistogram::Run(vtkh::DataSet &data_set, const std::string &field_name, vtkh::DataSet &grad_mag_ds)
+{
+  VTKH_DATA_OPEN("histogram");
+  VTKH_DATA_ADD("device", GetCurrentDevice());
+  VTKH_DATA_ADD("bins", m_num_bins);
+  VTKH_DATA_ADD("input_cells", data_set.GetNumberOfCells());
+  VTKH_DATA_ADD("input_domains", data_set.GetNumberOfDomains());
+
+  if(!data_set.GlobalFieldExists(field_name))
+  {
+    throw Error("Histogram: field '"+field_name+"' does not exist");
+  }
+
+  vtkm::Range val_range;
+  vtkm::cont::ArrayHandle<vtkm::Range> val_ranges = data_set.GetGlobalRange(field_name);
+  if(val_ranges.GetNumberOfValues() != 1)
+  {
+    throw Error("Histogram: field must have a single component");
+  }
+  val_range = val_ranges.ReadPortal().Get(0);
+
+  vtkm::Range grad_range;
+  vtkm::cont::ArrayHandle<vtkm::Range> grad_ranges = grad_mag_ds.GetGlobalRange("mag");
+  if(grad_ranges.GetNumberOfValues() != 1)
+  {
+    throw Error("Grad Mag Histogram: field must have a single component");
+  }
+  grad_range = grad_ranges.ReadPortal().Get(0);
+
+  // std::cout << "[new]scalar val global range: " << val_range << std::endl;
+  // std::cout << "[new]gradient mag global range: " << grad_range << std::endl;
+
+  const int num_domains = data_set.GetNumberOfDomains(); 
+
+
+  std::vector<ValGradHistogram::SparseHistogramResult> local_sparse_histograms;
+  for(int i = 0; i < num_domains; ++i)
+  {
+    vtkm::Id val_domain_id, grad_domain_id;
+    vtkm::cont::DataSet val_dom, grad_dom;
+    
+    data_set.GetDomain(i, val_dom, val_domain_id);
+    grad_mag_ds.GetDomain(i, grad_dom, grad_domain_id);
+
+    if(!val_dom.HasField(field_name)) continue;
+    vtkm::cont::Field val_field = val_dom.GetField(field_name);
+    
+    if(!grad_dom.HasField("mag")) continue;
+    vtkm::cont::Field grad_field = grad_dom.GetField("mag");
+    
+    
+    std::cout << "VGhist test:::" << "val_field nov = " << val_field.GetNumberOfValues() << "  grad_field nov = " << grad_field.GetNumberOfValues() << "\n";
+    if(val_field.GetNumberOfValues() != grad_field.GetNumberOfValues() )
+    {
+      throw Error("value and gradient magintude should have same number of datapoints");
+    }
+
+    vtkm::worklet::NDimsHistogram ndHistogram;
+    // Set the number of data points
+    ndHistogram.SetNumOfDataPoints(val_field.GetNumberOfValues());
+    vtkm::Float64 delta_val;
+    ndHistogram.AddField(val_field.GetData(), m_num_bins, val_range, delta_val, true);
+    vtkm::Float64 delta_grad;
+    ndHistogram.AddField(grad_field.GetData(), m_num_bins, grad_range, delta_grad, true);
+
+
+    std::vector<vtkm::cont::ArrayHandle<vtkm::Id>> binIds;
+    vtkm::cont::ArrayHandle<vtkm::Id> freqs;
+    ndHistogram.Run(binIds, freqs);
+
+    ValGradHistogram::SparseHistogramResult dom_2d_hist;
+    dom_2d_hist.m_binIds = binIds;
+    dom_2d_hist.m_freqs = freqs;
+    dom_2d_hist.m_nob = m_num_bins;
+    dom_2d_hist.m_range_val = val_range;
+    dom_2d_hist.m_range_grad = grad_range;
+    dom_2d_hist.m_bin_delta_val = delta_val;
+    dom_2d_hist.m_bin_delta_grad = delta_grad;
+
+    //dom_2d_hist.Print(std::cout);
+    local_sparse_histograms.push_back(dom_2d_hist);
+
+
+
+  }
+
+  ValGradHistogram::NonSparseHistogramResult local = valgrad_detail::merge_histograms(local_sparse_histograms);
+
+  
+  vtkm::Id * bin_ptr = GetVTKMPointer(local.m_freqs);
+  valgrad_detail::reduce(bin_ptr, m_num_bins*m_num_bins);
+
+  VTKH_DATA_CLOSE();
+  return local;
+}
+
+void
+ValGradHistogram::SparseHistogramResult::Print(std::ostream &out)
+{
+  vtkm::Id nonSparseBins = m_binIds[0].ReadPortal().GetNumberOfValues();
+
+  vtkm::Id f_sum = 0;
+  for (int i = 0; i < nonSparseBins; i++)
+  {
+    vtkm::Id idx0 = m_binIds[0].ReadPortal().Get(i);
+    vtkm::Id idx1 = m_binIds[1].ReadPortal().Get(i);
+    vtkm::Id f = m_freqs.ReadPortal().Get(i);
+    f_sum += f;
+    out << "[" << idx0 << ", " << idx1 << "] --> " << f << "\n";
+  }
+  out << "Total num of points in 2d histogram = " << f_sum << "\n"; 
+}
+
+void
+ValGradHistogram::NonSparseHistogramResult::Print(std::ostream &out)
+{
+
+  auto freqPortal = m_freqs.ReadPortal();  
+  
+  vtkm::Id f_sum = 0;
+  for(int i=0; i<m_nob; ++i)
+  {
+    for(int j=0; j<m_nob; ++j)
+    {
+      out << freqPortal.Get(i*m_nob + j) << "  ";
+      f_sum += freqPortal.Get(i*m_nob + j);
+    }
+    out << "\n";
+  }
+
+  out << "Total num of points in 2d histogram = " << f_sum << "\n"; 
+}
+
+vtkm::Id
+ValGradHistogram::NonSparseHistogramResult::totalCount()
+{
+  auto freqPortal = m_freqs.ReadPortal();  
+  vtkm::Id f_sum = 0;
+  for(int i=0; i<m_nob*m_nob; ++i)
+  {
+   f_sum += freqPortal.Get(i);
+  }
+  return f_sum;
+}
+
+Histogram::HistogramResult
+ValGradHistogram::NonSparseHistogramResult::getValDist()
+{
+  auto freqPortal = m_freqs.ReadPortal();
+  std::vector<vtkm::Id> marginal_freqs(m_nob, 0);
+  for(int i=0; i<m_nob; ++i)
+  {
+    for(int j=0; j<m_nob; ++j)
+    {
+      marginal_freqs[i] += freqPortal.Get(i*m_nob + j);
+    }
+  }
+
+  Histogram::HistogramResult marginal_hist;
+  marginal_hist.m_bins = vtkm::cont::make_ArrayHandle(marginal_freqs, vtkm::CopyFlag::On); 
+  marginal_hist.m_range = m_range_val;
+  marginal_hist.m_bin_delta = m_bin_delta_val;
+
+  return marginal_hist;
+}
+
+Histogram::HistogramResult
+ValGradHistogram::NonSparseHistogramResult::getGradDist()
+{
+  auto freqPortal = m_freqs.ReadPortal();
+  std::vector<vtkm::Id> marginal_freqs(m_nob, 0);
+  for(int i=0; i<m_nob; ++i)
+  {
+    for(int j=0; j<m_nob; ++j)
+    {
+      marginal_freqs[j] += freqPortal.Get(i*m_nob + j);
+    }
+  }
+
+  Histogram::HistogramResult marginal_hist;
+  marginal_hist.m_bins = vtkm::cont::make_ArrayHandle(marginal_freqs, vtkm::CopyFlag::On); 
+  marginal_hist.m_range = m_range_grad;
+  marginal_hist.m_bin_delta = m_bin_delta_grad;
+
+  return marginal_hist;
+}
+
+Histogram::HistogramResult
+ValGradHistogram::NonSparseHistogramResult::getConditionalGradDist(int val_bin_id)
+{
+  if(val_bin_id < 0 || val_bin_id >= m_nob)
+  {
+    throw Error("Conditional Grad Histogram: value bin id out of range");
+  }
+
+  auto freqPortal = m_freqs.ReadPortal();
+  std::vector<vtkm::Id> conditional_freqs(m_nob, 0);
+  for(int j=0; j<m_nob; ++j)
+  {
+    conditional_freqs[j] = freqPortal.Get(val_bin_id*m_nob + j);
+  }
+
+  Histogram::HistogramResult conditional_hist;
+  conditional_hist.m_bins = vtkm::cont::make_ArrayHandle(conditional_freqs, vtkm::CopyFlag::On); 
+  conditional_hist.m_range = m_range_grad;
+  conditional_hist.m_bin_delta = m_bin_delta_grad;
+
+  return conditional_hist;
+}
+
+} //  namespace vtkh
diff --git a/src/vtkh/filters/ValGradHistogram.hpp b/src/vtkh/filters/ValGradHistogram.hpp
new file mode 100644
index 00000000..d60a2cf9
--- /dev/null
+++ b/src/vtkh/filters/ValGradHistogram.hpp
@@ -0,0 +1,55 @@
+#ifndef VTK_H_VALGRADHISTOGRAM_HPP
+#define VTK_H_VALGRADHISTOGRAM_HPP
+
+#include <vtkh/vtkh.hpp>
+#include <vtkh/vtkh_exports.h>
+#include <vtkh/filters/Filter.hpp>
+#include <vtkh/DataSet.hpp>
+#include <vtkh/filters/Histogram.hpp>
+
+#include <vector>
+#include <iostream>
+
+namespace vtkh
+{
+
+class VTKH_API ValGradHistogram
+{
+public:
+  ValGradHistogram();
+  virtual ~ValGradHistogram();
+
+  struct SparseHistogramResult
+  {
+    std::vector<vtkm::cont::ArrayHandle<vtkm::Id>> m_binIds;
+    vtkm::cont::ArrayHandle<vtkm::Id> m_freqs;
+    int m_nob;
+    
+    vtkm::Range m_range_val, m_range_grad;
+    vtkm::Float64 m_bin_delta_val, m_bin_delta_grad;
+    void Print(std::ostream &out);
+  };
+
+  struct NonSparseHistogramResult
+  {
+    //this is a flattened 2d histogram
+    vtkm::cont::ArrayHandle<vtkm::Id> m_freqs;
+    int m_nob;
+    
+    vtkm::Range m_range_val, m_range_grad;
+    vtkm::Float64 m_bin_delta_val, m_bin_delta_grad;
+    void Print(std::ostream &out);
+    Histogram::HistogramResult getValDist();
+    Histogram::HistogramResult getGradDist();
+    Histogram::HistogramResult getConditionalGradDist(int val_bin_id);
+    vtkm::Id totalCount();
+  };
+
+  NonSparseHistogramResult Run(vtkh::DataSet &data_set, const std::string &field_name, vtkh::DataSet &grad_mag_ds);
+  void SetNumBins(const int num_bins);
+protected:
+  int m_num_bins;
+};
+
+} //namespace vtkh
+#endif