add histsampling filter

docs/changelog/hist-sampling.md

@@ -0,0 +1,4 @@
+
+# Added a `HistSampling` filter
+
+This filter assumes the field data are point clouds. It samples the field data according to its importance level. The importance level (sampling rate) is computed based on the histogram. The rarer values can provide more importance. More details can be found in the following paper. “In Situ Data-Driven Adaptive Sampling for Large-scale Simulation Data Summarization”, Ayan Biswas, Soumya Dutta, Jesus Pulido, and James Ahrens, In Situ Infrastructures for Enabling Extreme-scale Analysis and Visualization (ISAV 2018), co-located with Supercomputing 2018.

vtkm/filter/resampling/CMakeLists.txt

@@ -9,10 +9,13 @@
 ##============================================================================
 set(resampling_headers
   Probe.h
+  HistSampling.h
   )
 
 set(resampling_sources
-  Probe.cxx)
+  Probe.cxx
+  HistSampling.cxx
+  )
 
 vtkm_library(
   NAME vtkm_filter_resampling

vtkm/filter/resampling/HistSampling.cxx

@@ -0,0 +1,113 @@
+//============================================================================
+//  Copyright (c) Kitware, Inc.
+//  All rights reserved.
+//  See LICENSE.txt for details.
+//
+//  This software is distributed WITHOUT ANY WARRANTY; without even
+//  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+//  PURPOSE.  See the above copyright notice for more information.
+//============================================================================
+#include <vtkm/cont/Algorithm.h>
+#include <vtkm/cont/ArrayCopy.h>
+#include <vtkm/cont/ArrayHandleRandomUniformReal.h>
+#include <vtkm/filter/MapFieldPermutation.h>
+#include <vtkm/filter/density_estimate/Histogram.h>
+#include <vtkm/filter/entity_extraction/ThresholdPoints.h>
+#include <vtkm/filter/resampling/HistSampling.h>
+#include <vtkm/filter/resampling/worklet/HistSampling.h>
+#include <vtkm/worklet/DispatcherMapField.h>
+#include <vtkm/worklet/WorkletMapField.h>
+
+namespace vtkm
+{
+namespace filter
+{
+namespace resampling
+{
+namespace
+{
+vtkm::cont::ArrayHandle<vtkm::FloatDefault> CalculatPdf(vtkm::Id TotalPoints,
+                                                        vtkm::FloatDefault SamplePercent,
+                                                        vtkm::cont::ArrayHandle<vtkm::Id> BinCount)
+{
+  vtkm::Id NumBins = BinCount.GetNumberOfValues();
+  vtkm::cont::ArrayHandleIndex indexArray(NumBins);
+  vtkm::cont::ArrayHandle<vtkm::Id> BinIndices;
+  vtkm::cont::Algorithm::Copy(indexArray, BinIndices);
+  vtkm::cont::Algorithm::SortByKey(BinCount, BinIndices);
+
+  vtkm::FloatDefault remainingSamples = SamplePercent * TotalPoints;
+  vtkm::FloatDefault remainingBins = static_cast<vtkm::FloatDefault>(NumBins);
+  vtkm::cont::ArrayHandle<vtkm::FloatDefault> targetSamples;
+  targetSamples.Allocate(NumBins);
+
+  auto binCountPortal = BinCount.ReadPortal();
+  auto targetWritePortal = targetSamples.WritePortal();
+
+  for (int i = 0; i < NumBins; ++i)
+  {
+    vtkm::FloatDefault targetNeededSamples = remainingSamples / remainingBins;
+    vtkm::FloatDefault curCount = static_cast<vtkm::FloatDefault>(binCountPortal.Get(i));
+    vtkm::FloatDefault samplesTaken = vtkm::Min(curCount, targetNeededSamples);
+    targetWritePortal.Set(i, samplesTaken);
+    remainingBins = remainingBins - 1;
+    remainingSamples = remainingSamples - samplesTaken;
+  }
+
+  vtkm::cont::ArrayHandle<vtkm::FloatDefault> acceptanceProbsVec;
+  acceptanceProbsVec.AllocateAndFill(NumBins, -1.f);
+
+  vtkm::cont::Invoker invoker;
+  invoker(vtkm::worklet::AcceptanceProbsWorklet{},
+          targetSamples,
+          BinIndices,
+          BinCount,
+          acceptanceProbsVec);
+  return acceptanceProbsVec;
+}
+
+} // anonymous namespace
+
+vtkm::cont::DataSet HistSampling::DoExecute(const vtkm::cont::DataSet& input)
+{
+  //computing histogram based on input
+  vtkm::filter::density_estimate::Histogram histogram;
+  histogram.SetNumberOfBins(this->NumberOfBins);
+  histogram.SetActiveField(this->GetActiveFieldName());
+  auto histogramOutput = histogram.Execute(input);
+  vtkm::cont::ArrayHandle<vtkm::Id> binCountArray;
+  vtkm::cont::ArrayCopyShallowIfPossible(
+    histogramOutput.GetField(histogram.GetOutputFieldName()).GetData(), binCountArray);
+  vtkm::Id TotalPoints = input.GetNumberOfPoints();
+  //computing pdf
+  vtkm::cont::ArrayHandle<vtkm::FloatDefault> probArray;
+  probArray = CalculatPdf(TotalPoints, this->SamplePercent, binCountArray);
+  // use the acceptance probabilities and random array to create 0-1 array
+  // generating random array between 0 to 1
+  vtkm::cont::ArrayHandle<vtkm::Int8> outputArray;
+  auto resolveType = [&](const auto& concrete) {
+    vtkm::Id NumFieldValues = concrete.GetNumberOfValues();
+    auto randArray = vtkm::cont::ArrayHandleRandomUniformReal<vtkm::FloatDefault>(
+      NumFieldValues, { this->GetSeed() });
+    vtkm::worklet::DispatcherMapField<vtkm::worklet::LookupWorklet>(
+      vtkm::worklet::LookupWorklet{
+        this->NumberOfBins, histogram.GetComputedRange().Min, histogram.GetBinDelta() })
+      .Invoke(concrete, outputArray, probArray, randArray);
+  };
+  const auto& inField = this->GetFieldFromDataSet(input);
+  this->CastAndCallScalarField(inField, resolveType);
+
+  vtkm::cont::DataSet sampledDataSet =
+    this->CreateResultField(input, "ifsampling", inField.GetAssociation(), outputArray);
+  vtkm::filter::entity_extraction::ThresholdPoints threshold;
+  threshold.SetActiveField("ifsampling");
+  threshold.SetCompactPoints(true);
+  threshold.SetThresholdAbove(0.5);
+  // filter out the results with zero in it
+  vtkm::cont::DataSet thresholdDataSet = threshold.Execute(sampledDataSet);
+  return thresholdDataSet;
+}
+
+} // namespace resampling
+} // namespace filter
+} // namespace vtkm

vtkm/filter/resampling/HistSampling.h

@@ -0,0 +1,51 @@
+//============================================================================
+//  Copyright (c) Kitware, Inc.
+//  All rights reserved.
+//  See LICENSE.txt for details.
+//
+//  This software is distributed WITHOUT ANY WARRANTY; without even
+//  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+//  PURPOSE.  See the above copyright notice for more information.
+//============================================================================
+#ifndef vtk_m_filter_resampling_HistSampling_h
+#define vtk_m_filter_resampling_HistSampling_h
+
+#include <vtkm/filter/FilterField.h>
+#include <vtkm/filter/resampling/vtkm_filter_resampling_export.h>
+
+namespace vtkm
+{
+namespace filter
+{
+namespace resampling
+{
+// This filter can sample particles according to its importance level
+// The source code of this filter comes from
+// https://github.com/Alpine-DAV/ascent/blob/develop/src/libs/vtkh/filters/HistSampling.cpp
+// More details can be found in the following paper:
+// "In Situ Data-Driven Adaptive Sampling for Large-scale Simulation Data Summarization",
+// Ayan Biswas, Soumya Dutta, Jesus Pulido, and James Ahrens, In Situ Infrastructures for Enabling Extreme-scale Analysis and Visualization (ISAV 2018), co-located with Supercomputing 2018
+class VTKM_FILTER_RESAMPLING_EXPORT HistSampling : public vtkm::filter::FilterField
+{
+public:
+  VTKM_CONT void SetNumberOfBins(vtkm::Id numberOfBins) { this->NumberOfBins = numberOfBins; }
+  VTKM_CONT vtkm::Id GetNumberOfBins() { return this->NumberOfBins; }
+  VTKM_CONT void SetSamplePercent(vtkm::FloatDefault samplePercent)
+  {
+    this->SamplePercent = samplePercent;
+  }
+  VTKM_CONT vtkm::FloatDefault GetSamplePercent() { return this->SamplePercent; }
+  VTKM_CONT vtkm::UInt32 GetSeed() { return this->Seed; }
+  VTKM_CONT void SetSeed(vtkm::UInt32 seed) { this->Seed = seed; }
+
+private:
+  VTKM_CONT vtkm::cont::DataSet DoExecute(const vtkm::cont::DataSet& input) override;
+  vtkm::Id NumberOfBins;
+  vtkm::FloatDefault SamplePercent = static_cast<vtkm::FloatDefault>(0.1);
+  vtkm::UInt32 Seed = 0;
+};
+} // namespace resampling
+} // namespace filter
+} // namespace vtkm
+
+#endif // vtk_m_filter_resampling_HistSampling_h

vtkm/filter/resampling/testing/CMakeLists.txt

@@ -7,18 +7,20 @@
 ##  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
 ##  PURPOSE.  See the above copyright notice for more information.
 ##============================================================================
-
 set(unit_tests
   UnitTestProbe.cxx
   )
+set(unit_tests_device
+  UnitTestHistSampling.cxx
+  )
 
 set(libraries
-  vtkm_filter_clean_grid
   vtkm_filter_resampling
   )
 
 vtkm_unit_tests(
   SOURCES ${unit_tests}
+  DEVICE_SOURCES ${unit_tests_device}
   LIBRARIES ${libraries}
   USE_VTKM_JOB_POOL
 )

vtkm/filter/resampling/testing/UnitTestHistSampling.cxx

@@ -0,0 +1,92 @@
+//============================================================================
+//  Copyright (c) Kitware, Inc.
+//  All rights reserved.
+//  See LICENSE.txt for details.
+//
+//  This software is distributed WITHOUT ANY WARRANTY; without even
+//  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+//  PURPOSE.  See the above copyright notice for more information.
+//============================================================================
+#include <vtkm/cont/ArrayCopy.h>
+#include <vtkm/cont/DataSet.h>
+#include <vtkm/cont/DataSetBuilderUniform.h>
+#include <vtkm/cont/Invoker.h>
+#include <vtkm/cont/testing/Testing.h>
+#include <vtkm/filter/entity_extraction/ThresholdPoints.h>
+#include <vtkm/filter/resampling/HistSampling.h>
+#include <vtkm/worklet/WorkletMapField.h>
+namespace
+{
+
+struct CreateFieldValueWorklet : public vtkm::worklet::WorkletMapField
+{
+  vtkm::Id SizePerDim;
+  VTKM_CONT CreateFieldValueWorklet(vtkm::Id sizePerDim)
+    : SizePerDim(sizePerDim)
+  {
+  }
+  using ControlSignature = void(FieldOut);
+  using ExecutionSignature = void(_1, InputIndex);
+  template <typename OutPutType>
+  VTKM_EXEC void operator()(OutPutType& val, vtkm::Id idx) const
+  {
+    vtkm::Id x = idx % this->SizePerDim;
+    vtkm::Id y = (idx / this->SizePerDim) % this->SizePerDim;
+    vtkm::Id z = idx / this->SizePerDim / this->SizePerDim;
+    vtkm::FloatDefault center = static_cast<vtkm::FloatDefault>((1.0 * this->SizePerDim) / 2.0);
+    vtkm::FloatDefault v = vtkm::Pow((static_cast<vtkm::FloatDefault>(x) - center), 2) +
+      vtkm::Pow((static_cast<vtkm::FloatDefault>(y) - center), 2) +
+      vtkm::Pow((static_cast<vtkm::FloatDefault>(z) - center), 2);
+    if (v < 0.5)
+    {
+      val = static_cast<vtkm::FloatDefault>(10.0);
+      return;
+    }
+    val = static_cast<vtkm::FloatDefault>(10.0 / vtkm::Sqrt(v));
+  };
+};
+
+void TestHistSamplingSingleBlock()
+{
+  //creating data set for testing
+  vtkm::cont::DataSetBuilderUniform dsb;
+  constexpr vtkm::Id sizePerDim = 20;
+  constexpr vtkm::Id3 dimensions(sizePerDim, sizePerDim, sizePerDim);
+  vtkm::cont::DataSet dataSet = dsb.Create(dimensions);
+
+  //creating data set for testing
+  vtkm::cont::ArrayHandle<vtkm::FloatDefault> scalarArray;
+  scalarArray.Allocate(sizePerDim * sizePerDim * sizePerDim);
+  vtkm::cont::Invoker invoker;
+  invoker(CreateFieldValueWorklet{ sizePerDim }, scalarArray);
+  dataSet.AddPointField("scalarField", scalarArray);
+
+  //calling filter
+  using AsscoType = vtkm::cont::Field::Association;
+  vtkm::filter::resampling::HistSampling histsample;
+  histsample.SetNumberOfBins(10);
+  histsample.SetActiveField("scalarField", AsscoType::Points);
+  auto outputDataSet = histsample.Execute(dataSet);
+
+  //checking data sets to make sure all rared region are sampled
+  //are there better way to test it?
+  vtkm::filter::entity_extraction::ThresholdPoints threshold;
+  threshold.SetActiveField("scalarField");
+  threshold.SetCompactPoints(true);
+  threshold.SetThresholdAbove(9.9);
+  vtkm::cont::DataSet thresholdDataSet = threshold.Execute(outputDataSet);
+  //There are 7 points that have the scalar value of 10
+  VTKM_TEST_ASSERT(thresholdDataSet.GetField("scalarField").GetNumberOfValues() == 7);
+}
+
+void TestHistSampling()
+{
+  TestHistSamplingSingleBlock();
+} // TestHistSampling
+
+}
+
+int UnitTestHistSampling(int argc, char* argv[])
+{
+  return vtkm::cont::testing::Testing::Run(TestHistSampling, argc, argv);
+}

vtkm/filter/resampling/vtkm.module

@@ -4,6 +4,8 @@ GROUPS
   Filters
 DEPENDS
   vtkm_filter_core
+  vtkm_filter_density_estimate
+  vtkm_filter_entity_extraction
 PRIVATE_DEPENDS
   vtkm_worklet
 TEST_DEPENDS

vtkm/filter/resampling/worklet/CMakeLists.txt

@@ -7,9 +7,9 @@
 ##  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
 ##  PURPOSE.  See the above copyright notice for more information.
 ##============================================================================
-
 set(headers
   Probe.h
+  HistSampling.h
   )
 
 vtkm_declare_headers(${headers})

vtkm/filter/resampling/worklet/HistSampling.h

@@ -0,0 +1,79 @@
+//============================================================================
+//  Copyright (c) Kitware, Inc.
+//  All rights reserved.
+//  See LICENSE.txt for details.
+//
+//  This software is distributed WITHOUT ANY WARRANTY; without even
+//  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+//  PURPOSE.  See the above copyright notice for more information.
+//============================================================================
+#ifndef vtk_m_worklet_HistSampling_h
+#define vtk_m_worklet_HistSampling_h
+
+#include <vtkm/worklet/WorkletMapField.h>
+
+namespace vtkm
+{
+namespace worklet
+{
+class AcceptanceProbsWorklet : public vtkm::worklet::WorkletMapField
+{
+public:
+  using ControlSignature = void(FieldIn, FieldIn, FieldIn, WholeArrayOut);
+  using ExecutionSignature = void(_1, _2, _3, _4);
+  template <typename TypeOutPortal>
+  VTKM_EXEC void operator()(const vtkm::FloatDefault& targetSampleNum,
+                            const vtkm::Id& binIndex,
+                            const vtkm::Id& binCount,
+                            TypeOutPortal arrayOutPortal) const
+  {
+    if (binCount < 1 || targetSampleNum < 0.000001)
+    {
+      arrayOutPortal.Set(binIndex, 0.0);
+    }
+    else
+    {
+      arrayOutPortal.Set(binIndex, targetSampleNum / static_cast<vtkm::FloatDefault>(binCount));
+    }
+  }
+};
+
+class LookupWorklet : public vtkm::worklet::WorkletMapField
+{
+private:
+  vtkm::Id m_num_bins;
+  vtkm::Float64 m_min;
+  vtkm::Float64 m_bin_delta;
+
+public:
+  LookupWorklet(vtkm::Id num_bins, vtkm::Float64 min_value, 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);
+  template <typename TablePortal, typename FieldType>
+  VTKM_EXEC vtkm::FloatDefault operator()(const FieldType& field_value,
+                                          TablePortal table,
+                                          const vtkm::FloatDefault& random) 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;
+    }
+    return random < table.Get(bin);
+  }
+};
+};
+} // vtkm::worklet
+
+
+#endif // vtk_m_worklet_HistSampling_h