Shuenhoy/vtk-m2
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Neareast Grid Point particle density estimate

Browse Source
This commit is contained in:
Li-Ta Lo 2020-10-19 15:29:14 -06:00
parent 8c15745c04
commit 9c4129bfb1
4 changed files with 73 additions and 36 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

8
vtkm/filter/ParticleDensityNGP.h
View File

@ -24,10 +24,12 @@ class ParticleDensityNGP : public vtkm::filter::FilterField<ParticleDensityNGP>
{
public:
// ParticleDensity only support turning 2D/3D particle positions into density
using SupportedTypes = vtkm::TypeListFloatVec;
using SupportedTypes = vtkm::TypeListFieldVec3;
//
ParticleDensityNGP(vtkm::Id3& dimension, vtkm::Vec3f& origin, vtkm::Vec3f& spacing);
ParticleDensityNGP(const vtkm::Id3& dimension,
const vtkm::Vec3f& origin,
const vtkm::Vec3f& spacing);
template <typename T, typename StorageType, typename Policy>
VTKM_CONT vtkm::cont::DataSet DoExecute(const vtkm::cont::DataSet& input,
@ -36,7 +38,7 @@ public:
vtkm::filter::PolicyBase<Policy> policy);
private:
vtkm::Id3 Dimension;
vtkm::Id3 Dimension; // Cell dimension
vtkm::Vec3f Origin;
vtkm::Vec3f Spacing;
};

66
vtkm/filter/ParticleDensityNGP.hxx
View File

@ -30,16 +30,24 @@ public:
template <typename Point, typename CellLocatorExecObj, typename AtomicArray>
VTKM_EXEC void operator()(const Point& point,
const CellLocatorExecObj& locator,
AtomicArray& density)
AtomicArray& density) const
{
vtkm::Id cellId;
vtkm::Id cellId{};
vtkm::Vec3f parametric;
// Find the cell containing the point
locator.FindCell(point, cellId, parametric, *this);
if (locator->FindCell(point, cellId, parametric) == vtkm::ErrorCode::Success)
{
// increment density
density.Add(cellId, 1);
}
else
{
// increment density
density.Add(cellId, 1);
// FIXME: what does mean when it is not found?
// We simply ignore that particular particle.
std::cout << "WTF: " << point << std::endl;
}
}
}; //NGPWorklet
} //worklet
@ -50,9 +58,9 @@ namespace vtkm
{
namespace filter
{
inline VTKM_CONT ParticleDensityNGP::ParticleDensityNGP(vtkm::Id3& dimension,
vtkm::Vec3f& origin,
vtkm::Vec3f& spacing)
inline VTKM_CONT ParticleDensityNGP::ParticleDensityNGP(const vtkm::Id3& dimension,
const vtkm::Vec3f& origin,
const vtkm::Vec3f& spacing)
: Dimension(dimension)
, Origin(origin)
, Spacing(spacing)
@ -61,41 +69,39 @@ inline VTKM_CONT ParticleDensityNGP::ParticleDensityNGP(vtkm::Id3& dimension,
template <typename T, typename StorageType, typename Policy>
inline VTKM_CONT vtkm::cont::DataSet ParticleDensityNGP::DoExecute(
const vtkm::cont::DataSet& input,
const vtkm::cont::DataSet&,
const vtkm::cont::ArrayHandle<T, StorageType>& field,
const vtkm::filter::FieldMetadata& fieldMeta,
vtkm::filter::PolicyBase<Policy> policy)
const vtkm::filter::FieldMetadata&,
vtkm::filter::PolicyBase<Policy>)
{
// TODO: it really doesn't need to be a UniformGrid, any CellSet with CellLocator will work.
// Make it another input rather an output generated.
// Make it another input rather an output generated.
// We want to stores density as PointField which conforms to VTK/VTKm's idea of ImageDataset
// and works with the ImageConnectivity for segmentation purpose. We thus create a surrogate
// uniform dataset that has the cell dimension as the point dimension of the output. We use
// this dataset only for point in cell lookup purpose. This is a somewhat convolved way of
// doing some simple math.
auto surrogate = vtkm::cont::DataSetBuilderUniform::Create(
this->Dimension - vtkm::Id3{ 1, 1, 1 }, this->Origin, this->Spacing);
// We stores density as CellField which conforms to physicists' idea of particle density
// better. However, VTK/VTKm's idea of "Image" Dataset and the ImageConnectivity filter
// expect a PointField. For better separation of concerns, we create a uniform dataset
// that has the cell dimension as expected and later convert the dataset to its dual.
auto uniform = vtkm::cont::DataSetBuilderUniform::Create(
this->Dimension + vtkm::Id3{ 1, 1, 1 }, this->Origin, this->Spacing);
// Create a CellSetLocator
// Create a CellLocator
vtkm::cont::CellLocatorUniformGrid locator;
locator.SetCellSet(surrogate.GetCellSet());
locator.SetCoordinates(surrogate.GetCoordinateSystem());
locator.SetCellSet(uniform.GetCellSet());
locator.SetCoordinates(uniform.GetCoordinateSystem());
locator.Update();
// We still use an ArrayHandle<T> and pass it to the Worklet as AtomicArrayInOut
// it will be type converted automatically. However the ArrayHandle needs to be
// allocated and initialized. The easily way to do it is to copy from an
// ArrayHandleConstant
// We create an ArrayHandle and pass it to the Worklet as AtomicArrayInOut.
// However the ArrayHandle needs to be allocated and initialized first. The
// easily way to do it is to copy from an ArrayHandleConstant
vtkm::cont::ArrayHandle<vtkm::Id> density;
vtkm::cont::ArrayCopy(vtkm::cont::ArrayHandleConstant<vtkm::Id>(0, surrogate.GetNumberOfCells()),
vtkm::cont::ArrayCopy(vtkm::cont::ArrayHandleConstant<vtkm::Id>(0, uniform.GetNumberOfCells()),
density);
this->Invoke(input, locator, density);
this->Invoke(vtkm::worklet::NGPWorklet{}, field, locator, density);
surrogate.AddField(vtkm::cont::make_FieldCell("density", density));
uniform.AddField(vtkm::cont::make_FieldCell("density", density));
return surrogate;
return uniform;
}
}

1
vtkm/filter/testing/CMakeLists.txt
View File

@ -47,6 +47,7 @@ set(unit_tests
UnitTestNDEntropyFilter.cxx
UnitTestNDHistogramFilter.cxx
UnitTestParticleAdvectionFilter.cxx
UnitTestParticleDensity.cxx
UnitTestPartitionedDataSetFilters.cxx
UnitTestPartitionedDataSetHistogramFilter.cxx
UnitTestPointAverageFilter.cxx

34
vtkm/filter/testing/UnitTestParticleDensity.cxx
View File

@ -8,14 +8,42 @@
// PURPOSE. See the above copyright notice for more information.
//============================================================================
#include <vtkm/cont/ArrayHandleRandomUniformReal.h>
#include <vtkm/cont/DataSetBuilderExplicit.h>
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/filter/ParticleDensityNGP.h>
#include <vector>
#include <vtkm/worklet/DescriptiveStatistics.h>
void TestNGP()
{
std::vector<vtkm::Vec3f> positions = { { 0.5, 0.5 } };
const vtkm::Id N = 1000000;
auto composite = vtkm::cont::make_ArrayHandleCompositeVector(
vtkm::cont::ArrayHandleRandomUniformReal<vtkm::Float32>(N, 0xceed),
vtkm::cont::ArrayHandleRandomUniformReal<vtkm::Float32>(N, 0xdeed),
vtkm::cont::ArrayHandleRandomUniformReal<vtkm::Float32>(N, 0xabba));
vtkm::cont::ArrayHandle<vtkm::Vec3f> positions;
vtkm::cont::ArrayCopy(composite, positions);
vtkm::cont::ArrayHandle<vtkm::Id> connectivity;
vtkm::cont::ArrayCopy(vtkm::cont::make_ArrayHandleIndex(N), connectivity);
auto dataSet = vtkm::cont::DataSetBuilderExplicit::Create(
positions, vtkm::CellShapeTagVertex{}, 1, connectivity);
auto cellDims = vtkm::Id3{ 3, 3, 3 };
vtkm::filter::ParticleDensityNGP filter{ cellDims,
{ 0.f, 0.f, 0.f },
vtkm::Vec3f{ 1.f / 3.f, 1.f / 3.f, 1.f / 3.f } };
filter.SetUseCoordinateSystemAsField(true);
auto density = filter.Execute(dataSet);
vtkm::cont::ArrayHandle<vtkm::Id> field;
density.GetCellField("density").GetData().AsArrayHandle<vtkm::Id>(field);
auto field_f = vtkm::cont::make_ArrayHandleCast<vtkm::Float32>(field);
auto result = vtkm::worklet::DescriptiveStatistics::Run(field_f);
VTKM_TEST_ASSERT(test_equal(result.Sum(), N));
VTKM_TEST_ASSERT(test_equal(result.Mean(), N / density.GetNumberOfCells()));
}
void TestParticleDensity()