From 55a08d51b72470799183c5a78ce6019bb9c434bc Mon Sep 17 00:00:00 2001 From: Nrushad Joshi Date: Wed, 7 Jun 2023 11:29:46 -0400 Subject: [PATCH] Isosurface Uncertainty Visualization Filter --- benchmarking/BenchmarkODEIntegrators.cxx | 2 +- docs/changelog/isosurface-uncertainty.md | 7 + .../testing/UnitTestLogValues.cxx | 1 + vtkm/filter/uncertainty/CMakeLists.txt | 25 ++ .../uncertainty/ContourUncertainUniform.cxx | 205 ++++++++++++++++ .../uncertainty/ContourUncertainUniform.h | 118 ++++++++++ .../ContourUncertainUniformMonteCarlo.cxx | 222 ++++++++++++++++++ .../ContourUncertainUniformMonteCarlo.h | 79 +++++++ .../filter/uncertainty/testing/CMakeLists.txt | 25 ++ .../UnitTestContourUncertainUniform.cxx | 158 +++++++++++++ vtkm/filter/uncertainty/vtkm.module | 13 + 11 files changed, 854 insertions(+), 1 deletion(-) create mode 100644 docs/changelog/isosurface-uncertainty.md create mode 100644 vtkm/filter/uncertainty/CMakeLists.txt create mode 100644 vtkm/filter/uncertainty/ContourUncertainUniform.cxx create mode 100644 vtkm/filter/uncertainty/ContourUncertainUniform.h create mode 100644 vtkm/filter/uncertainty/ContourUncertainUniformMonteCarlo.cxx create mode 100644 vtkm/filter/uncertainty/ContourUncertainUniformMonteCarlo.h create mode 100644 vtkm/filter/uncertainty/testing/CMakeLists.txt create mode 100644 vtkm/filter/uncertainty/testing/UnitTestContourUncertainUniform.cxx create mode 100644 vtkm/filter/uncertainty/vtkm.module diff --git a/benchmarking/BenchmarkODEIntegrators.cxx b/benchmarking/BenchmarkODEIntegrators.cxx index c3848a289..ac64eb6db 100644 --- a/benchmarking/BenchmarkODEIntegrators.cxx +++ b/benchmarking/BenchmarkODEIntegrators.cxx @@ -33,7 +33,7 @@ void BenchParticleAdvection(::benchmark::State& state) vtkm::Id numPoints = dims[0] * dims[1] * dims[2]; - std::vector vectorField(static_cast(numPoints)); + std::vector vectorField(static_cast(numPoints)); // 3D for (std::size_t i = 0; i < static_cast(numPoints); i++) vectorField[i] = vecX; diff --git a/docs/changelog/isosurface-uncertainty.md b/docs/changelog/isosurface-uncertainty.md new file mode 100644 index 000000000..403bb8d4c --- /dev/null +++ b/docs/changelog/isosurface-uncertainty.md @@ -0,0 +1,7 @@ +# New Isosurface Uncertainty Visualization Analysis filter + +This new filter, designed for uncertainty analysis of the marching cubes algorithm for isosurface visualization, computes three uncertainty metrics, namely, level-crossing probability, topology case count, and entropy-based uncertainty. This filter requires a 3D cell-structured dataset with ensemble minimum and ensemble maximum values denoting uncertain data range at each grid vertex. + +This uncertainty analysis filter analyzes the uncertainty in the marching cubes topology cases. The output dataset consists of 3D point fields corresponding to the level-crossing probability, topology case count, and entropy-based uncertainty. + +This VTK-m implementation is based on the algorithm presented in the paper "FunMC2: A Filter for Uncertainty Visualization of Marching Cubes on Multi-Core Devices" by Wang, Z., Athawale, T. M., Moreland, K., Chen, J., Johnson, C. R., & Pugmire, D. diff --git a/vtkm/filter/field_transform/testing/UnitTestLogValues.cxx b/vtkm/filter/field_transform/testing/UnitTestLogValues.cxx index bdd748547..cd13c6c84 100644 --- a/vtkm/filter/field_transform/testing/UnitTestLogValues.cxx +++ b/vtkm/filter/field_transform/testing/UnitTestLogValues.cxx @@ -7,6 +7,7 @@ // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR // PURPOSE. See the above copyright notice for more information. //============================================================================ + #include #include diff --git a/vtkm/filter/uncertainty/CMakeLists.txt b/vtkm/filter/uncertainty/CMakeLists.txt new file mode 100644 index 000000000..76d25fa35 --- /dev/null +++ b/vtkm/filter/uncertainty/CMakeLists.txt @@ -0,0 +1,25 @@ +##============================================================================ +## 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. +##============================================================================ + +set(uncertainty_headers + ContourUncertainUniform.h + ContourUncertainUniformMonteCarlo.h + ) +set(uncertainty_sources_device + ContourUncertainUniform.cxx + ContourUncertainUniformMonteCarlo.cxx + ) + +vtkm_library( + NAME vtkm_filter_uncertainty + HEADERS ${uncertainty_headers} + DEVICE_SOURCES ${uncertainty_sources_device} + USE_VTKM_JOB_POOL +) diff --git a/vtkm/filter/uncertainty/ContourUncertainUniform.cxx b/vtkm/filter/uncertainty/ContourUncertainUniform.cxx new file mode 100644 index 000000000..904100f48 --- /dev/null +++ b/vtkm/filter/uncertainty/ContourUncertainUniform.cxx @@ -0,0 +1,205 @@ +//============================================================================ +// 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. +//============================================================================ + +// This code is based on the algorithm presented in the following papers: +// Wang, J., Athawale, T., Moreland, K., Chen, J., Johnson, C., & Pugmire, +// D. (2023). FunMC^ 2: A Filter for Uncertainty Visualization of Marching +// Cubes on Multi-Core Devices. Oak Ridge National Laboratory (ORNL), +// Oak Ridge, TN (United States). +// +// Athawale, T. M., Sane, S., & Johnson, C. R. (2021, October). Uncertainty +// Visualization of the Marching Squares and Marching Cubes Topology Cases. +// In 2021 IEEE Visualization Conference (VIS) (pp. 106-110). IEEE. + +#include +#include +#include +#include +#include + +namespace +{ +class ClosedFormUniform : public vtkm::worklet::WorkletVisitCellsWithPoints +{ +public: + ClosedFormUniform(double isovalue) + : m_isovalue(isovalue){}; + using ControlSignature = + void(CellSetIn, FieldInPoint, FieldInPoint, FieldOutCell, FieldOutCell, FieldOutCell); + + using ExecutionSignature = void(_2, _3, _4, _5, _6); + + using InputDomain = _1; + template + + VTKM_EXEC void operator()(const InPointFieldMinType& inPointFieldVecMin, + const InPointFieldMaxType& inPointFieldVecMax, + OutCellFieldType1& outCellFieldCProb, + OutCellFieldType2& outCellFieldNumNonzeroProb, + OutCellFieldType3& outCellFieldEntropy) const + { + vtkm::IdComponent numPoints = inPointFieldVecMin.GetNumberOfComponents(); + + if (numPoints != 8) + { + this->RaiseError("This is the 3D version for 8 vertices\n"); + return; + } + + vtkm::FloatDefault allPositiveProb = 1.0; + vtkm::FloatDefault allNegativeProb = 1.0; + vtkm::FloatDefault allCrossProb = 0.0; + vtkm::FloatDefault positiveProb; + vtkm::FloatDefault negativeProb; + vtkm::Vec ProbList; + + constexpr vtkm::IdComponent totalNumCases = 256; + vtkm::Vec probHistogram; + + for (vtkm::IdComponent pointIndex = 0; pointIndex < numPoints; ++pointIndex) + { + vtkm::FloatDefault minV = static_cast(inPointFieldVecMin[pointIndex]); + vtkm::FloatDefault maxV = static_cast(inPointFieldVecMax[pointIndex]); + + if (this->m_isovalue <= minV) + { + positiveProb = 1.0; + negativeProb = 0.0; + } + else if (this->m_isovalue >= maxV) + { + positiveProb = 0.0; + negativeProb = 1.0; + } + else + { + positiveProb = static_cast((maxV - (this->m_isovalue)) / (maxV - minV)); + negativeProb = 1 - positiveProb; + } + + allNegativeProb *= negativeProb; + allPositiveProb *= positiveProb; + + ProbList[pointIndex][0] = negativeProb; + ProbList[pointIndex][1] = positiveProb; + } + + allCrossProb = 1 - allPositiveProb - allNegativeProb; + outCellFieldCProb = allCrossProb; + + TraverseBit(ProbList, probHistogram); + + vtkm::FloatDefault entropyValue = 0; + vtkm::Id nonzeroCases = 0; + vtkm::FloatDefault templog = 0; + + for (vtkm::IdComponent i = 0; i < totalNumCases; i++) + { + templog = 0; + if (probHistogram[i] > 0.00001) + { + nonzeroCases++; + templog = vtkm::Log2(probHistogram[i]); + } + entropyValue = entropyValue + (-probHistogram[i]) * templog; + } + + outCellFieldNumNonzeroProb = nonzeroCases; + outCellFieldEntropy = entropyValue; + } + + VTKM_EXEC inline void TraverseBit(vtkm::Vec& ProbList, + vtkm::Vec& probHistogram) const + { + + for (vtkm::IdComponent i = 0; i < 256; i++) + { + vtkm::FloatDefault currProb = 1.0; + for (vtkm::IdComponent j = 0; j < 8; j++) + { + if (i & (1 << j)) + { + currProb *= ProbList[j][1]; + } + else + { + currProb *= ProbList[j][0]; + } + } + probHistogram[i] = currProb; + } + } + +private: + double m_isovalue; +}; +} + +namespace vtkm +{ +namespace filter +{ +namespace uncertainty +{ +ContourUncertainUniform::ContourUncertainUniform() +{ + this->SetCrossProbabilityName("cross_probability"); +} +VTKM_CONT vtkm::cont::DataSet ContourUncertainUniform::DoExecute(const vtkm::cont::DataSet& input) +{ + vtkm::cont::Field minField = this->GetFieldFromDataSet(0, input); + vtkm::cont::Field maxField = this->GetFieldFromDataSet(1, input); + + vtkm::cont::UnknownArrayHandle crossProbability; + vtkm::cont::UnknownArrayHandle numNonZeroProbability; + vtkm::cont::UnknownArrayHandle entropy; + + if (!input.GetCellSet().IsType>()) + { + throw vtkm::cont::ErrorBadType("Uncertain contour only works for CellSetStructured<3>."); + } + vtkm::cont::CellSetStructured<3> cellSet; + input.GetCellSet().AsCellSet(cellSet); + + auto resolveType = [&](auto concreteMinField) { + using ArrayType = std::decay_t; + using ValueType = typename ArrayType::ValueType; + ArrayType concreteMaxField; + vtkm::cont::ArrayCopyShallowIfPossible(maxField.GetData(), concreteMaxField); + + vtkm::cont::ArrayHandle concreteCrossProb; + vtkm::cont::ArrayHandle concreteNumNonZeroProb; + vtkm::cont::ArrayHandle concreteEntropy; + this->Invoke(ClosedFormUniform{ this->IsoValue }, + cellSet, + concreteMinField, + concreteMaxField, + concreteCrossProb, + concreteNumNonZeroProb, + concreteEntropy); + crossProbability = concreteCrossProb; + numNonZeroProbability = concreteNumNonZeroProb; + entropy = concreteEntropy; + }; + this->CastAndCallScalarField(minField, resolveType); + + vtkm::cont::DataSet result = this->CreateResult(input); + result.AddCellField(this->GetCrossProbabilityName(), crossProbability); + result.AddCellField(this->GetNumberNonzeroProbabilityName(), numNonZeroProbability); + result.AddCellField(this->GetEntropyName(), entropy); + return result; +} +} +} +} diff --git a/vtkm/filter/uncertainty/ContourUncertainUniform.h b/vtkm/filter/uncertainty/ContourUncertainUniform.h new file mode 100644 index 000000000..ed08ba5ed --- /dev/null +++ b/vtkm/filter/uncertainty/ContourUncertainUniform.h @@ -0,0 +1,118 @@ +//============================================================================ +// 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. +//============================================================================ + +// This code is based on the algorithm presented in the following papers: +// Wang, J., Athawale, T., Moreland, K., Chen, J., Johnson, C., & Pugmire, +// D. (2023). FunMC^ 2: A Filter for Uncertainty Visualization of Marching +// Cubes on Multi-Core Devices. Oak Ridge National Laboratory (ORNL), +// Oak Ridge, TN (United States). +// +// Athawale, T. M., Sane, S., & Johnson, C. R. (2021, October). Uncertainty +// Visualization of the Marching Squares and Marching Cubes Topology Cases. +// In 2021 IEEE Visualization Conference (VIS) (pp. 106-110). IEEE. + +#ifndef vtk_m_filter_uncertainty_ContourUncertainUniform_h +#define vtk_m_filter_uncertainty_ContourUncertainUniform_h + +#include +#include + +namespace vtkm +{ +namespace filter +{ +namespace uncertainty +{ +/// @brief Visualize isosurface uncertainty for uniform distributed data. +/// +/// This filter computes the positional uncertainty of isosurfaces as a +/// function of uncertainty in input data, where the data are assumed to +/// be uniformly distributed and sampled on a regular grid. The uniform +/// distribution range is given through the input datasets via the minimum +/// and maximum fields. Given the uniform distribution range, the computed +/// isosurface uncertainty corresponds to uncertainty in topology cases in +/// the marching cubes algorithm. +/// +class VTKM_FILTER_UNCERTAINTY_EXPORT ContourUncertainUniform : public vtkm::filter::FilterField +{ + + std::string NumberNonzeroProbabilityName = "num_nonzero_probability"; + std::string EntropyName = "entropy"; + vtkm::Float64 IsoValue = 0.0; + +public: + /// @brief Constructor + VTKM_CONT ContourUncertainUniform(); + + /// @brief Sets minimum field. + /// Sets minimum value of uniform distribution at each grid point. + VTKM_CONT void SetMinField(const std::string& fieldName) + { + this->SetActiveField(0, fieldName, vtkm::cont::Field::Association::Points); + } + + /// @brief Sets maximum field. + /// Sets maximum value of uniform distribution at each grid point. + VTKM_CONT void SetMaxField(const std::string& fieldName) + { + this->SetActiveField(1, fieldName, vtkm::cont::Field::Association::Points); + } + + /// @brief Sets isovalue. + /// Sets isovalue for extracting isosurfaces. + VTKM_CONT void SetIsoValue(vtkm::Float64 value) { this->IsoValue = value; } + + /// @brief Gets isovalue + /// Gets isovalue used for visualizing isosurfaces + VTKM_CONT vtkm::Float64 GetIsoValue() const { return this->IsoValue; } + + /// @brief Sets crossing probability field (uncertainty field type 1). + /// Sets the output field name that stores isosurface crossing probabiliy for each grid cell. + VTKM_CONT void SetCrossProbabilityName(const std::string& name) + { + this->SetOutputFieldName(name); + } + + /// @brief Gets crossing probability field (uncertainty field type 1). + /// Gets the output field name that stores isosurface crossing probability for each grid cell. + VTKM_CONT const std::string& GetCrossProbabilityName() const + { + return this->GetOutputFieldName(); + } + + /// @brief Sets toplogy case count field (uncertainty field type 2). + /// Sets the output field name that stores the number of marching cubes toplogy cases for each grid cell. + VTKM_CONT void SetNumberNonzeroProbabilityName(const std::string& name) + { + this->NumberNonzeroProbabilityName = name; + } + + /// @brief Gets toplogy case count field (uncertainty field type 2. + /// Gets the output field name that stores the number of marching cubes toplogy cases for each grid cell. + VTKM_CONT const std::string& GetNumberNonzeroProbabilityName() const + { + return this->NumberNonzeroProbabilityName; + } + + /// @brief Sets entropy field. (uncertainty field type 3) + /// Sets the output field name that stores the entropy of a histogram of marching cubes toplogy cases. + VTKM_CONT void SetEntropyName(const std::string& name) { this->EntropyName = name; } + + /// @brief Gets entropy field. (uncertainty field type 3) + /// Gets the output field name that stores the entropy of a histogram of marching cubes toplogy cases. + VTKM_CONT const std::string& GetEntropyName() const { return this->EntropyName; } + +protected: + VTKM_CONT vtkm::cont::DataSet DoExecute(const vtkm::cont::DataSet& input) override; +}; +} +} +} +#endif diff --git a/vtkm/filter/uncertainty/ContourUncertainUniformMonteCarlo.cxx b/vtkm/filter/uncertainty/ContourUncertainUniformMonteCarlo.cxx new file mode 100644 index 000000000..ca7ba70c4 --- /dev/null +++ b/vtkm/filter/uncertainty/ContourUncertainUniformMonteCarlo.cxx @@ -0,0 +1,222 @@ +//============================================================================ +// 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 + +namespace vtkm +{ +namespace worklet +{ +namespace uniform +{ +class ContourUncertainUniformMonteCarlo : public vtkm::worklet::WorkletVisitCellsWithPoints +{ + +public: + ContourUncertainUniformMonteCarlo(double isovalue, int itervalue) + : m_isovalue(isovalue) + , m_numsample(itervalue) + { + } + using ControlSignature = void(CellSetIn, + FieldInPoint, + FieldInPoint, + FieldInCell, + FieldOutCell, + FieldOutCell, + FieldOutCell); + using ExecutionSignature = void(_2, _3, _4, _5, _6, _7); + using InputDomain = _1; + + template + + VTKM_EXEC void operator()(const InPointFieldMinType& inPointFieldVecMin, + const InPointFieldMaxType& inPointFieldVecMax, + const InRandomNumbersType& randomNumbers, + OutCellFieldType1& outNonCrossProb, + OutCellFieldType2& outCrossProb, + OutCellFieldType3& outEntropyProb) const + { + vtkm::IdComponent numPoints = inPointFieldVecMin.GetNumberOfComponents(); + + if (numPoints != 8) + { + this->RaiseError("This is the 3D version for 8 vertices\n"); + return; + } + + vtkm::FloatDefault minV = 0.0; + vtkm::FloatDefault maxV = 0.0; + vtkm::FloatDefault uniformDistValue = 0.0; + vtkm::IdComponent numSample = this->m_numsample; + vtkm::FloatDefault numCrossing = 0; + vtkm::FloatDefault crossProb = 0; + + vtkm::IdComponent zeroFlag; + vtkm::IdComponent oneFlag; + vtkm::Float64 base = 2.0; + vtkm::Float64 totalSum = 0.0; + vtkm::IdComponent nonZeroCase = 0; + vtkm::FloatDefault entropyValue = 0; + vtkm::FloatDefault templog = 0; + vtkm::FloatDefault value = 0.0; + vtkm::IdComponent k = 0; + + constexpr vtkm::IdComponent totalNumCases = 256; + vtkm::Vec probHistogram; + + for (vtkm::IdComponent j = 0; j < totalNumCases; j++) + { + probHistogram[j] = 0; + } + + for (vtkm::IdComponent i = 0; i < numSample; ++i) + { + zeroFlag = 0; + oneFlag = 0; + totalSum = 0.0; + for (vtkm::IdComponent pointIndex = 0; pointIndex < numPoints; ++pointIndex) + { + minV = static_cast(inPointFieldVecMin[pointIndex]); + maxV = static_cast(inPointFieldVecMax[pointIndex]); + + auto tmp = randomNumbers[k]; + + uniformDistValue = minV + tmp * (maxV - minV); + + if (uniformDistValue <= this->m_isovalue) + { + zeroFlag = 1; + } + else + { + oneFlag = 1; + totalSum += vtkm::Pow(base, pointIndex); + } + + k += 1; + } + + if ((oneFlag == 1) && (zeroFlag == 1)) + { + numCrossing += 1; + } + + if ((totalSum >= 0) && (totalSum <= 255)) + { + probHistogram[static_cast(totalSum)] += 1; + } + } + + for (vtkm::IdComponent i = 0; i < totalNumCases; i++) + { + templog = 0; + value = static_cast(probHistogram[i] / + static_cast(numSample)); + if (probHistogram[i] > 0.00001) + { + nonZeroCase++; + templog = vtkm::Log2(value); + } + entropyValue = entropyValue + (-value) * templog; + } + + crossProb = numCrossing / static_cast(numSample); + outNonCrossProb = static_cast(nonZeroCase); + outCrossProb = crossProb; + outEntropyProb = entropyValue; + } + +private: + double m_isovalue; + int m_numsample; +}; +} +} +} + +namespace vtkm +{ +namespace filter +{ +namespace uncertainty +{ +ContourUncertainUniformMonteCarlo::ContourUncertainUniformMonteCarlo() +{ + this->SetCrossProbabilityName("cross_probability"); +} +VTKM_CONT vtkm::cont::DataSet ContourUncertainUniformMonteCarlo::DoExecute( + const vtkm::cont::DataSet& input) +{ + vtkm::cont::Field minField = this->GetFieldFromDataSet(0, input); + vtkm::cont::Field maxField = this->GetFieldFromDataSet(1, input); + + vtkm::cont::UnknownArrayHandle crossProbability; + vtkm::cont::UnknownArrayHandle nonCrossProbability; + vtkm::cont::UnknownArrayHandle entropyProbability; + + if (!input.GetCellSet().IsType>()) + { + throw vtkm::cont::ErrorBadType("Uncertain contour only works for CellSetStructured<3>."); + } + vtkm::cont::CellSetStructured<3> cellSet; + input.GetCellSet().AsCellSet(cellSet); + + auto resolveType = [&](auto concreteMinField) { + using ArrayType = std::decay_t; + using ValueType = typename ArrayType::ValueType; + ArrayType concreteMaxField; + vtkm::cont::ArrayCopyShallowIfPossible(maxField.GetData(), concreteMaxField); + + vtkm::cont::ArrayHandle concreteCrossProb; + vtkm::cont::ArrayHandle concreteNonCrossProb; + vtkm::cont::ArrayHandle concreteEntropyProb; + + vtkm::cont::ArrayHandleRandomUniformReal randomArray( + cellSet.GetNumberOfCells() * this->IterValue * 8, { 0xceed }); + + this->Invoke( + vtkm::worklet::uniform::ContourUncertainUniformMonteCarlo{ this->IsoValue, this->IterValue }, + cellSet, + concreteMinField, + concreteMaxField, + vtkm::cont::make_ArrayHandleGroupVecVariable( + randomArray, + vtkm::cont::ArrayHandleCounting( + 0, this->IterValue * 8, cellSet.GetNumberOfCells() + 1)), + concreteNonCrossProb, + concreteCrossProb, + concreteEntropyProb); + + crossProbability = concreteCrossProb; + nonCrossProbability = concreteNonCrossProb; + entropyProbability = concreteEntropyProb; + }; + this->CastAndCallScalarField(minField, resolveType); + vtkm::cont::DataSet result = this->CreateResult(input); + result.AddCellField(this->GetCrossProbabilityName(), crossProbability); + result.AddCellField(this->GetNumberNonzeroProbabilityName(), nonCrossProbability); + result.AddCellField(this->GetEntropyName(), entropyProbability); + return result; +} +} +} +} diff --git a/vtkm/filter/uncertainty/ContourUncertainUniformMonteCarlo.h b/vtkm/filter/uncertainty/ContourUncertainUniformMonteCarlo.h new file mode 100644 index 000000000..673b1b94d --- /dev/null +++ b/vtkm/filter/uncertainty/ContourUncertainUniformMonteCarlo.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_filter_uncertainty_ContourUncertainUniformMonteCarlo_h +#define vtk_m_filter_uncertainty_ContourUncertainUniformMonteCarlo_h +#include +#include + +namespace vtkm +{ +namespace filter +{ +namespace uncertainty +{ +/// \brief Visualize isosurface uncertainty using Monte Carlo approach for uniformly distributed data. +/// +/// This filter is implemented to validate the correctness of the ContourUncertainUniform filter. +/// We encourage usage of the ContourUncertainUniform filter because the Monte Carlo approach implemented +/// in this filter is computationally inefficient. +/// +class VTKM_FILTER_UNCERTAINTY_EXPORT ContourUncertainUniformMonteCarlo + : public vtkm::filter::FilterField +{ + std::string NumberNonzeroProbabilityName = "num_nonzero_probability"; + std::string EntropyName = "entropy"; + vtkm::Float64 IsoValue = 0.0; + vtkm::IdComponent IterValue = 1; + +public: + VTKM_CONT ContourUncertainUniformMonteCarlo(); + + VTKM_CONT void SetMinField(const std::string& fieldName) + { + this->SetActiveField(0, fieldName, vtkm::cont::Field::Association::Points); + } + VTKM_CONT void SetMaxField(const std::string& fieldName) + { + this->SetActiveField(1, fieldName, vtkm::cont::Field::Association::Points); + } + VTKM_CONT void SetIsoValue(vtkm::Float64 value) { this->IsoValue = value; } + VTKM_CONT vtkm::Float64 GetIsoValue() const { return this->IsoValue; } + + VTKM_CONT void SetNumSample(vtkm::IdComponent value) { this->IterValue = value; } + VTKM_CONT vtkm::IdComponent GetNumSample() const { return this->IterValue; } + + VTKM_CONT void SetCrossProbabilityName(const std::string& name) + { + this->SetOutputFieldName(name); + } + VTKM_CONT const std::string& GetCrossProbabilityName() const + { + return this->GetOutputFieldName(); + } + + VTKM_CONT void SetNumberNonzeroProbabilityName(const std::string& name) + { + this->NumberNonzeroProbabilityName = name; + } + VTKM_CONT const std::string& GetNumberNonzeroProbabilityName() const + { + return this->NumberNonzeroProbabilityName; + } + VTKM_CONT void SetEntropyName(const std::string& name) { this->EntropyName = name; } + VTKM_CONT const std::string& GetEntropyName() const { return this->EntropyName; } + +protected: + VTKM_CONT vtkm::cont::DataSet DoExecute(const vtkm::cont::DataSet& input) override; +}; +} +} +} +#endif diff --git a/vtkm/filter/uncertainty/testing/CMakeLists.txt b/vtkm/filter/uncertainty/testing/CMakeLists.txt new file mode 100644 index 000000000..16b9583e7 --- /dev/null +++ b/vtkm/filter/uncertainty/testing/CMakeLists.txt @@ -0,0 +1,25 @@ +##============================================================================ +## 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. +##============================================================================ + +set(libraries + vtkm_filter_uncertainty + vtkm_filter_core + vtkm_filter_vector_analysis + ) + +set(unit_tests + UnitTestContourUncertainUniform.cxx + ) + +vtkm_unit_tests( + SOURCES ${unit_tests} + LIBRARIES ${libraries} + USE_VTKM_JOB_POOL + ) diff --git a/vtkm/filter/uncertainty/testing/UnitTestContourUncertainUniform.cxx b/vtkm/filter/uncertainty/testing/UnitTestContourUncertainUniform.cxx new file mode 100644 index 000000000..241624ff2 --- /dev/null +++ b/vtkm/filter/uncertainty/testing/UnitTestContourUncertainUniform.cxx @@ -0,0 +1,158 @@ +// +//============================================================================ +// 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 +#include + +namespace +{ +template +vtkm::cont::DataSet MakeContourUncertainUniformTestDataSet() +{ + const vtkm::Id3 dims(25, 25, 25); + + vtkm::Id numPoints = dims[0] * dims[1] * dims[2]; + vtkm::cont::DataSetBuilderUniform dataSetBuilder; + vtkm::cont::DataSet dataSet = dataSetBuilder.Create(dims); + + std::vector ensemble_max; + std::vector ensemble_min; + + vtkm::IdComponent k = 0; + vtkm::cont::ArrayHandleRandomUniformReal randomArray(25 * 25 * 25 * 2, + { 0xceed }); + auto portal = randomArray.ReadPortal(); + + for (vtkm::Id i = 0; i < numPoints; ++i) + { + vtkm::FloatDefault value1 = -20 + (40 * portal.Get(k)); + vtkm::FloatDefault value2 = -20 + (40 * portal.Get(k + 1)); + ensemble_max.push_back(static_cast(vtkm::Max(value1, value2))); + ensemble_min.push_back(static_cast(vtkm::Min(value1, value2))); + k += 2; + } + + dataSet.AddPointField("ensemble_max", ensemble_max); + dataSet.AddPointField("ensemble_min", ensemble_min); + return dataSet; +} + +void TestUncertaintyGeneral(vtkm::FloatDefault isoValue) +{ + // Isosurface uncertainty computation using closed form solution + vtkm::cont::DataSet input = MakeContourUncertainUniformTestDataSet(); + vtkm::filter::uncertainty::ContourUncertainUniform filter; + filter.SetIsoValue(isoValue); + filter.SetCrossProbabilityName("CrossProbablity"); + filter.SetNumberNonzeroProbabilityName("NonzeroProbablity"); + filter.SetEntropyName("Entropy"); + filter.SetMinField("ensemble_min"); + filter.SetMaxField("ensemble_max"); + vtkm::cont::DataSet output = filter.Execute(input); + + // Isosurface uncertainty computation using Monte Carlo solution + vtkm::filter::uncertainty::ContourUncertainUniformMonteCarlo filter_mc; + filter_mc.SetIsoValue(isoValue); + filter_mc.SetNumSample(1000); + filter_mc.SetCrossProbabilityName("CrossProbablityMC"); + filter_mc.SetNumberNonzeroProbabilityName("NonzeroProbablityMC"); + filter_mc.SetEntropyName("EntropyMC"); + filter_mc.SetMinField("ensemble_min"); + filter_mc.SetMaxField("ensemble_max"); + vtkm::cont::DataSet output_mc = filter_mc.Execute(input); + + // Crossing Probablity field (closed form) + vtkm::cont::Field CrossProb = output.GetField("CrossProbablity"); + vtkm::cont::ArrayHandle crossProbArray; + CrossProb.GetData().AsArrayHandle(crossProbArray); + vtkm::cont::ArrayHandle::ReadPortalType CrossPortal = + crossProbArray.ReadPortal(); + + // Nonzero Value field (closed form) + vtkm::cont::Field NonzeroProb = output.GetField("NonzeroProbablity"); + vtkm::cont::ArrayHandle NonzeroProbArray; + NonzeroProb.GetData().AsArrayHandle(NonzeroProbArray); + vtkm::cont::ArrayHandle::ReadPortalType NonzeroPortal = NonzeroProbArray.ReadPortal(); + + // Entropy field (closed form) + vtkm::cont::Field entropy = output.GetField("Entropy"); + vtkm::cont::ArrayHandle EntropyArray; + entropy.GetData().AsArrayHandle(EntropyArray); + vtkm::cont::ArrayHandle::ReadPortalType EntropyPortal = + EntropyArray.ReadPortal(); + + // Crossing Probablity field (Marching Cube) + vtkm::cont::Field CrossProbMC = output_mc.GetField("CrossProbablityMC"); + vtkm::cont::ArrayHandle crossProbMCArray; + CrossProbMC.GetData().AsArrayHandle(crossProbMCArray); + vtkm::cont::ArrayHandle::ReadPortalType CrossMCPortal = + crossProbMCArray.ReadPortal(); + + // Nonzero Value field (Marching Cube) + vtkm::cont::Field NonzeroMCProb = output_mc.GetField("NonzeroProbablityMC"); + vtkm::cont::ArrayHandle NonzeroProbMCArray; + NonzeroMCProb.GetData().AsArrayHandle(NonzeroProbMCArray); + vtkm::cont::ArrayHandle::ReadPortalType NonzeroMCPortal = + NonzeroProbMCArray.ReadPortal(); + + // Entropy field (Marching Cube) + vtkm::cont::Field entropyMC = output_mc.GetField("EntropyMC"); + vtkm::cont::ArrayHandle EntropyMCArray; + entropyMC.GetData().AsArrayHandle(EntropyMCArray); + vtkm::cont::ArrayHandle::ReadPortalType EntropyMCPortal = + EntropyMCArray.ReadPortal(); + + + // Comparision of Closed-form and Marching Cube + for (vtkm::Id i = 0; i < crossProbArray.GetNumberOfValues(); ++i) + { + vtkm::FloatDefault CrossProbValue = CrossPortal.Get(i); + vtkm::Id NonzeroProbValue = NonzeroPortal.Get(i); // long long not float + vtkm::FloatDefault EntropyValue = EntropyPortal.Get(i); + + vtkm::FloatDefault CrossProbMCValue = CrossMCPortal.Get(i); + vtkm::FloatDefault NonzeroProbMCValue = NonzeroMCPortal.Get(i); + vtkm::FloatDefault EntropyMCValue = EntropyMCPortal.Get(i); + + // Maximum Entropy Difference: 8 + // Maximum Cross Probability Difference: 1 + // Maximum Nonzero Difference: 256 + VTKM_TEST_ASSERT((vtkm::Abs(CrossProbMCValue - CrossProbValue) < 0.1) || + (vtkm::Abs(NonzeroProbMCValue - NonzeroProbValue) < 50) || + (vtkm::Abs(EntropyMCValue - EntropyValue) < 0.5), + vtkm::Abs(CrossProbMCValue - CrossProbValue), + ' ', + vtkm::Abs(NonzeroProbMCValue - NonzeroProbValue) < 50, + ' ', + vtkm::Abs(EntropyMCValue - EntropyValue), + " No Match With Monte Carlo Sampling"); + } +} + +void TestContourUncertainUniform() +{ + vtkm::FloatDefault isoValue = 0; + TestUncertaintyGeneral(isoValue); +} +} +int UnitTestContourUncertainUniform(int argc, char* argv[]) +{ + return vtkm::cont::testing::Testing::Run(TestContourUncertainUniform, argc, argv); +} diff --git a/vtkm/filter/uncertainty/vtkm.module b/vtkm/filter/uncertainty/vtkm.module new file mode 100644 index 000000000..ddeb736f5 --- /dev/null +++ b/vtkm/filter/uncertainty/vtkm.module @@ -0,0 +1,13 @@ +NAME + vtkm_filter_uncertainty +GROUPS + FiltersCommon + Filters +DEPENDS + vtkm_filter_core +PRIVATE_DEPENDS + vtkm_worklet +TEST_DEPENDS + vtkm_filter_uncertainty + vtkm_filter_core + vtkm_filter_vector_analysis