add particle density

vtkm/filter/CMakeLists.txt

@@ -50,6 +50,7 @@ set(headers
   MeshQuality.h
   NDEntropy.h
   NDHistogram.h
+  ParticleDensityNGP.h
   Pathline.h
   PointAverage.h
   PointElevation.h
@@ -116,6 +117,7 @@ set(header_template_sources
   MeshQuality.hxx
   NDEntropy.hxx
   NDHistogram.hxx
+  ParticleDensityNGP.hxx
   Pathline.hxx
   PointAverage.hxx
   PointElevation.hxx

vtkm/filter/ParticleDensityNGP.h

@@ -0,0 +1,48 @@
+//============================================================================
+//  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_particle_density_ngp_h
+#define vtk_m_filter_particle_density_ngp_h
+
+#include <vtkm/filter/FilterField.h>
+
+namespace vtkm
+{
+namespace filter
+{
+/// \brief Estimate the density of particles using the Nearest Grid Point method
+
+// We only need the CoordinateSystem of the input dataset thus a FilterField
+class ParticleDensityNGP : public vtkm::filter::FilterField<ParticleDensityNGP>
+{
+public:
+  // ParticleDensity only support turning 2D/3D particle positions into density
+  using SupportedTypes = vtkm::TypeListFloatVec;
+
+  //
+  ParticleDensityNGP(vtkm::Id3& dimension, vtkm::Vec3f& origin, vtkm::Vec3f& spacing);
+
+  template <typename T, typename StorageType, typename Policy>
+  VTKM_CONT vtkm::cont::DataSet DoExecute(const vtkm::cont::DataSet& input,
+                                          const vtkm::cont::ArrayHandle<T, StorageType>& field,
+                                          const vtkm::filter::FieldMetadata& fieldMeta,
+                                          vtkm::filter::PolicyBase<Policy> policy);
+
+private:
+  vtkm::Id3 Dimension;
+  vtkm::Vec3f Origin;
+  vtkm::Vec3f Spacing;
+};
+}
+}
+
+#include <vtkm/filter/ParticleDensityNGP.hxx>
+
+#endif //vtk_m_filter_particle_density_ngp_h

vtkm/filter/ParticleDensityNGP.hxx

@@ -0,0 +1,102 @@
+//============================================================================
+//  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_particle_density_ngp_hxx
+#define vtk_m_filter_particle_density_ngp_hxx
+
+#include "ParticleDensityNGP.h"
+#include <vtkm/cont/CellLocatorUniformGrid.h>
+#include <vtkm/cont/DataSetBuilderUniform.h>
+#include <vtkm/filter/PolicyBase.h>
+#include <vtkm/worklet/WorkletMapField.h>
+
+namespace vtkm
+{
+namespace worklet
+{
+class NGPWorklet : public vtkm::worklet::WorkletMapField
+{
+public:
+  using ControlSignature = void(FieldIn coords, ExecObject locator, AtomicArrayInOut density);
+  using ExecutionSignature = void(_1, _2, _3);
+
+  template <typename Point, typename CellLocatorExecObj, typename AtomicArray>
+  VTKM_EXEC void operator()(const Point& point,
+                            const CellLocatorExecObj& locator,
+                            AtomicArray& density)
+  {
+    vtkm::Id cellId;
+    vtkm::Vec3f parametric;
+
+    // Find the cell containing the point
+    locator.FindCell(point, cellId, parametric, *this);
+
+    // increment density
+    density.Add(cellId, 1);
+  }
+}; //NGPWorklet
+} //worklet
+} //vtkm
+
+
+namespace vtkm
+{
+namespace filter
+{
+inline VTKM_CONT ParticleDensityNGP::ParticleDensityNGP(vtkm::Id3& dimension,
+                                                        vtkm::Vec3f& origin,
+                                                        vtkm::Vec3f& spacing)
+  : Dimension(dimension)
+  , Origin(origin)
+  , Spacing(spacing)
+{
+}
+
+template <typename T, typename StorageType, typename Policy>
+inline VTKM_CONT vtkm::cont::DataSet ParticleDensityNGP::DoExecute(
+  const vtkm::cont::DataSet& input,
+  const vtkm::cont::ArrayHandle<T, StorageType>& field,
+  const vtkm::filter::FieldMetadata& fieldMeta,
+  vtkm::filter::PolicyBase<Policy> 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.
+
+  // 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);
+
+  // Create a CellSetLocator
+  vtkm::cont::CellLocatorUniformGrid locator;
+  locator.SetCellSet(surrogate.GetCellSet());
+  locator.SetCoordinates(surrogate.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
+  vtkm::cont::ArrayHandle<vtkm::Id> density;
+  vtkm::cont::ArrayCopy(vtkm::cont::ArrayHandleConstant<vtkm::Id>(0, surrogate.GetNumberOfCells()),
+                        density);
+
+  this->Invoke(input, locator, density);
+
+  surrogate.AddField(vtkm::cont::make_FieldCell("density", density));
+
+  return surrogate;
+}
+}
+}
+#endif //vtk_m_filter_particle_density_ngp_hxx

vtkm/filter/testing/UnitTestParticleDensity.cxx

@@ -0,0 +1,29 @@
+//============================================================================
+//  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/testing/Testing.h>
+#include <vtkm/filter/ParticleDensityNGP.h>
+
+#include <vector>
+
+void TestNGP()
+{
+  std::vector<vtkm::Vec3f> positions = { { 0.5, 0.5 } };
+}
+
+void TestParticleDensity()
+{
+  TestNGP();
+}
+
+int UnitTestParticleDensity(int argc, char* argv[])
+{
+  return vtkm::cont::testing::Testing::Run(TestParticleDensity, argc, argv);
+}