Isosurface Uncertainty Visualization Filter
This commit is contained in:
parent
f4c6394264
commit
55a08d51b7
@ -33,7 +33,7 @@ void BenchParticleAdvection(::benchmark::State& state)
|
||||
|
||||
vtkm::Id numPoints = dims[0] * dims[1] * dims[2];
|
||||
|
||||
std::vector<vtkm::Vec3f> vectorField(static_cast<std::size_t>(numPoints));
|
||||
std::vector<vtkm::Vec3f> vectorField(static_cast<std::size_t>(numPoints)); // 3D
|
||||
for (std::size_t i = 0; i < static_cast<std::size_t>(numPoints); i++)
|
||||
vectorField[i] = vecX;
|
||||
|
||||
|
7
docs/changelog/isosurface-uncertainty.md
Normal file
7
docs/changelog/isosurface-uncertainty.md
Normal file
@ -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.
|
@ -7,6 +7,7 @@
|
||||
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
|
||||
// PURPOSE. See the above copyright notice for more information.
|
||||
//============================================================================
|
||||
|
||||
#include <vtkm/cont/testing/Testing.h>
|
||||
#include <vtkm/filter/field_transform/LogValues.h>
|
||||
|
||||
|
25
vtkm/filter/uncertainty/CMakeLists.txt
Normal file
25
vtkm/filter/uncertainty/CMakeLists.txt
Normal file
@ -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
|
||||
)
|
205
vtkm/filter/uncertainty/ContourUncertainUniform.cxx
Normal file
205
vtkm/filter/uncertainty/ContourUncertainUniform.cxx
Normal file
@ -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 <vtkm/cont/ArrayCopy.h>
|
||||
#include <vtkm/cont/DataSet.h>
|
||||
#include <vtkm/cont/Timer.h>
|
||||
#include <vtkm/filter/uncertainty/ContourUncertainUniform.h>
|
||||
#include <vtkm/worklet/WorkletMapTopology.h>
|
||||
|
||||
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 <typename InPointFieldMinType,
|
||||
typename InPointFieldMaxType,
|
||||
typename OutCellFieldType1,
|
||||
typename OutCellFieldType2,
|
||||
typename OutCellFieldType3>
|
||||
|
||||
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<vtkm::Vec2f, 8> ProbList;
|
||||
|
||||
constexpr vtkm::IdComponent totalNumCases = 256;
|
||||
vtkm::Vec<vtkm::FloatDefault, totalNumCases> probHistogram;
|
||||
|
||||
for (vtkm::IdComponent pointIndex = 0; pointIndex < numPoints; ++pointIndex)
|
||||
{
|
||||
vtkm::FloatDefault minV = static_cast<vtkm::FloatDefault>(inPointFieldVecMin[pointIndex]);
|
||||
vtkm::FloatDefault maxV = static_cast<vtkm::FloatDefault>(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<vtkm::FloatDefault>((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<vtkm::Vec2f, 8>& ProbList,
|
||||
vtkm::Vec<vtkm::FloatDefault, 256>& 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<vtkm::cont::CellSetStructured<3>>())
|
||||
{
|
||||
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<decltype(concreteMinField)>;
|
||||
using ValueType = typename ArrayType::ValueType;
|
||||
ArrayType concreteMaxField;
|
||||
vtkm::cont::ArrayCopyShallowIfPossible(maxField.GetData(), concreteMaxField);
|
||||
|
||||
vtkm::cont::ArrayHandle<ValueType> concreteCrossProb;
|
||||
vtkm::cont::ArrayHandle<vtkm::Id> concreteNumNonZeroProb;
|
||||
vtkm::cont::ArrayHandle<ValueType> 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;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
118
vtkm/filter/uncertainty/ContourUncertainUniform.h
Normal file
118
vtkm/filter/uncertainty/ContourUncertainUniform.h
Normal file
@ -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 <vtkm/filter/FilterField.h>
|
||||
#include <vtkm/filter/uncertainty/vtkm_filter_uncertainty_export.h>
|
||||
|
||||
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
|
222
vtkm/filter/uncertainty/ContourUncertainUniformMonteCarlo.cxx
Normal file
222
vtkm/filter/uncertainty/ContourUncertainUniformMonteCarlo.cxx
Normal file
@ -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 <vtkm/cont/ArrayCopy.h>
|
||||
#include <vtkm/cont/ArrayHandleGroupVecVariable.h>
|
||||
#include <vtkm/cont/ArrayHandleRandomUniformReal.h>
|
||||
#include <vtkm/cont/DataSet.h>
|
||||
#include <vtkm/cont/Timer.h>
|
||||
#include <vtkm/filter/uncertainty/ContourUncertainUniformMonteCarlo.h>
|
||||
#include <vtkm/worklet/WorkletMapTopology.h>
|
||||
|
||||
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 <typename InPointFieldMinType,
|
||||
typename InPointFieldMaxType,
|
||||
typename InRandomNumbersType,
|
||||
typename OutCellFieldType1,
|
||||
typename OutCellFieldType2,
|
||||
typename OutCellFieldType3>
|
||||
|
||||
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<vtkm::FloatDefault, totalNumCases> 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<vtkm::FloatDefault>(inPointFieldVecMin[pointIndex]);
|
||||
maxV = static_cast<vtkm::FloatDefault>(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<vtkm::IdComponent>(totalSum)] += 1;
|
||||
}
|
||||
}
|
||||
|
||||
for (vtkm::IdComponent i = 0; i < totalNumCases; i++)
|
||||
{
|
||||
templog = 0;
|
||||
value = static_cast<vtkm::FloatDefault>(probHistogram[i] /
|
||||
static_cast<vtkm::FloatDefault>(numSample));
|
||||
if (probHistogram[i] > 0.00001)
|
||||
{
|
||||
nonZeroCase++;
|
||||
templog = vtkm::Log2(value);
|
||||
}
|
||||
entropyValue = entropyValue + (-value) * templog;
|
||||
}
|
||||
|
||||
crossProb = numCrossing / static_cast<vtkm::FloatDefault>(numSample);
|
||||
outNonCrossProb = static_cast<vtkm::FloatDefault>(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<vtkm::cont::CellSetStructured<3>>())
|
||||
{
|
||||
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<decltype(concreteMinField)>;
|
||||
using ValueType = typename ArrayType::ValueType;
|
||||
ArrayType concreteMaxField;
|
||||
vtkm::cont::ArrayCopyShallowIfPossible(maxField.GetData(), concreteMaxField);
|
||||
|
||||
vtkm::cont::ArrayHandle<ValueType> concreteCrossProb;
|
||||
vtkm::cont::ArrayHandle<ValueType> concreteNonCrossProb;
|
||||
vtkm::cont::ArrayHandle<ValueType> concreteEntropyProb;
|
||||
|
||||
vtkm::cont::ArrayHandleRandomUniformReal<vtkm::FloatDefault> 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<vtkm::Id>(
|
||||
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;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
79
vtkm/filter/uncertainty/ContourUncertainUniformMonteCarlo.h
Normal file
79
vtkm/filter/uncertainty/ContourUncertainUniformMonteCarlo.h
Normal file
@ -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 <vtkm/filter/FilterField.h>
|
||||
#include <vtkm/filter/uncertainty/vtkm_filter_uncertainty_export.h>
|
||||
|
||||
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
|
25
vtkm/filter/uncertainty/testing/CMakeLists.txt
Normal file
25
vtkm/filter/uncertainty/testing/CMakeLists.txt
Normal file
@ -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
|
||||
)
|
@ -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 <iostream>
|
||||
#include <random>
|
||||
#include <vtkm/Math.h>
|
||||
#include <vtkm/cont/ArrayHandleRandomUniformReal.h>
|
||||
#include <vtkm/cont/DataSet.h>
|
||||
#include <vtkm/cont/DataSetBuilderUniform.h>
|
||||
#include <vtkm/cont/testing/Testing.h>
|
||||
#include <vtkm/filter/uncertainty/ContourUncertainUniform.h>
|
||||
#include <vtkm/filter/uncertainty/ContourUncertainUniformMonteCarlo.h>
|
||||
#include <vtkm/io/VTKDataSetReader.h>
|
||||
#include <vtkm/io/VTKDataSetWriter.h>
|
||||
|
||||
namespace
|
||||
{
|
||||
template <typename T>
|
||||
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<T> ensemble_max;
|
||||
std::vector<T> ensemble_min;
|
||||
|
||||
vtkm::IdComponent k = 0;
|
||||
vtkm::cont::ArrayHandleRandomUniformReal<vtkm::FloatDefault> 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<T>(vtkm::Max(value1, value2)));
|
||||
ensemble_min.push_back(static_cast<T>(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::FloatDefault>();
|
||||
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<vtkm::FloatDefault> crossProbArray;
|
||||
CrossProb.GetData().AsArrayHandle(crossProbArray);
|
||||
vtkm::cont::ArrayHandle<vtkm::FloatDefault>::ReadPortalType CrossPortal =
|
||||
crossProbArray.ReadPortal();
|
||||
|
||||
// Nonzero Value field (closed form)
|
||||
vtkm::cont::Field NonzeroProb = output.GetField("NonzeroProbablity");
|
||||
vtkm::cont::ArrayHandle<vtkm::Id> NonzeroProbArray;
|
||||
NonzeroProb.GetData().AsArrayHandle(NonzeroProbArray);
|
||||
vtkm::cont::ArrayHandle<vtkm::Id>::ReadPortalType NonzeroPortal = NonzeroProbArray.ReadPortal();
|
||||
|
||||
// Entropy field (closed form)
|
||||
vtkm::cont::Field entropy = output.GetField("Entropy");
|
||||
vtkm::cont::ArrayHandle<vtkm::FloatDefault> EntropyArray;
|
||||
entropy.GetData().AsArrayHandle(EntropyArray);
|
||||
vtkm::cont::ArrayHandle<vtkm::FloatDefault>::ReadPortalType EntropyPortal =
|
||||
EntropyArray.ReadPortal();
|
||||
|
||||
// Crossing Probablity field (Marching Cube)
|
||||
vtkm::cont::Field CrossProbMC = output_mc.GetField("CrossProbablityMC");
|
||||
vtkm::cont::ArrayHandle<vtkm::FloatDefault> crossProbMCArray;
|
||||
CrossProbMC.GetData().AsArrayHandle(crossProbMCArray);
|
||||
vtkm::cont::ArrayHandle<vtkm::FloatDefault>::ReadPortalType CrossMCPortal =
|
||||
crossProbMCArray.ReadPortal();
|
||||
|
||||
// Nonzero Value field (Marching Cube)
|
||||
vtkm::cont::Field NonzeroMCProb = output_mc.GetField("NonzeroProbablityMC");
|
||||
vtkm::cont::ArrayHandle<vtkm::FloatDefault> NonzeroProbMCArray;
|
||||
NonzeroMCProb.GetData().AsArrayHandle(NonzeroProbMCArray);
|
||||
vtkm::cont::ArrayHandle<vtkm::FloatDefault>::ReadPortalType NonzeroMCPortal =
|
||||
NonzeroProbMCArray.ReadPortal();
|
||||
|
||||
// Entropy field (Marching Cube)
|
||||
vtkm::cont::Field entropyMC = output_mc.GetField("EntropyMC");
|
||||
vtkm::cont::ArrayHandle<vtkm::FloatDefault> EntropyMCArray;
|
||||
entropyMC.GetData().AsArrayHandle(EntropyMCArray);
|
||||
vtkm::cont::ArrayHandle<vtkm::FloatDefault>::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);
|
||||
}
|
13
vtkm/filter/uncertainty/vtkm.module
Normal file
13
vtkm/filter/uncertainty/vtkm.module
Normal file
@ -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
|
Loading…
Reference in New Issue
Block a user