diff --git a/docs/changelog/hist-sampling.md b/docs/changelog/hist-sampling.md new file mode 100644 index 000000000..bb7dabb74 --- /dev/null +++ b/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. diff --git a/vtkm/filter/resampling/CMakeLists.txt b/vtkm/filter/resampling/CMakeLists.txt index 726d079a0..9376ae3e5 100644 --- a/vtkm/filter/resampling/CMakeLists.txt +++ b/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 diff --git a/vtkm/filter/resampling/HistSampling.cxx b/vtkm/filter/resampling/HistSampling.cxx new file mode 100644 index 000000000..561520225 --- /dev/null +++ b/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 +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace vtkm +{ +namespace filter +{ +namespace resampling +{ +namespace +{ +vtkm::cont::ArrayHandle CalculatPdf(vtkm::Id TotalPoints, + vtkm::FloatDefault SamplePercent, + vtkm::cont::ArrayHandle BinCount) +{ + vtkm::Id NumBins = BinCount.GetNumberOfValues(); + vtkm::cont::ArrayHandleIndex indexArray(NumBins); + vtkm::cont::ArrayHandle BinIndices; + vtkm::cont::Algorithm::Copy(indexArray, BinIndices); + vtkm::cont::Algorithm::SortByKey(BinCount, BinIndices); + + vtkm::FloatDefault remainingSamples = SamplePercent * TotalPoints; + vtkm::FloatDefault remainingBins = static_cast(NumBins); + vtkm::cont::ArrayHandle 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(binCountPortal.Get(i)); + vtkm::FloatDefault samplesTaken = vtkm::Min(curCount, targetNeededSamples); + targetWritePortal.Set(i, samplesTaken); + remainingBins = remainingBins - 1; + remainingSamples = remainingSamples - samplesTaken; + } + + vtkm::cont::ArrayHandle 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 binCountArray; + vtkm::cont::ArrayCopyShallowIfPossible( + histogramOutput.GetField(histogram.GetOutputFieldName()).GetData(), binCountArray); + vtkm::Id TotalPoints = input.GetNumberOfPoints(); + //computing pdf + vtkm::cont::ArrayHandle 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 outputArray; + auto resolveType = [&](const auto& concrete) { + vtkm::Id NumFieldValues = concrete.GetNumberOfValues(); + auto randArray = vtkm::cont::ArrayHandleRandomUniformReal( + NumFieldValues, { this->GetSeed() }); + vtkm::worklet::DispatcherMapField( + 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 diff --git a/vtkm/filter/resampling/HistSampling.h b/vtkm/filter/resampling/HistSampling.h new file mode 100644 index 000000000..1d037ab65 --- /dev/null +++ b/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 +#include + +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(0.1); + vtkm::UInt32 Seed = 0; +}; +} // namespace resampling +} // namespace filter +} // namespace vtkm + +#endif // vtk_m_filter_resampling_HistSampling_h diff --git a/vtkm/filter/resampling/testing/CMakeLists.txt b/vtkm/filter/resampling/testing/CMakeLists.txt index 299260f36..e354afe46 100644 --- a/vtkm/filter/resampling/testing/CMakeLists.txt +++ b/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 ) diff --git a/vtkm/filter/resampling/testing/UnitTestHistSampling.cxx b/vtkm/filter/resampling/testing/UnitTestHistSampling.cxx new file mode 100644 index 000000000..1f7add472 --- /dev/null +++ b/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 +#include +#include +#include +#include +#include +#include +#include +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 + 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((1.0 * this->SizePerDim) / 2.0); + vtkm::FloatDefault v = vtkm::Pow((static_cast(x) - center), 2) + + vtkm::Pow((static_cast(y) - center), 2) + + vtkm::Pow((static_cast(z) - center), 2); + if (v < 0.5) + { + val = static_cast(10.0); + return; + } + val = static_cast(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 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); +} diff --git a/vtkm/filter/resampling/vtkm.module b/vtkm/filter/resampling/vtkm.module index 6b17adf26..74a04254f 100644 --- a/vtkm/filter/resampling/vtkm.module +++ b/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 diff --git a/vtkm/filter/resampling/worklet/CMakeLists.txt b/vtkm/filter/resampling/worklet/CMakeLists.txt index 3aa9b2d67..007e6d74b 100644 --- a/vtkm/filter/resampling/worklet/CMakeLists.txt +++ b/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}) diff --git a/vtkm/filter/resampling/worklet/HistSampling.h b/vtkm/filter/resampling/worklet/HistSampling.h new file mode 100644 index 000000000..08c25648b --- /dev/null +++ b/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 + +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 + 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(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 + VTKM_EXEC vtkm::FloatDefault operator()(const FieldType& field_value, + TablePortal table, + const vtkm::FloatDefault& random) const + { + vtkm::Id bin = static_cast((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