added Moments calcuation example

examples/moments/CMakeLists.txt

@@ -0,0 +1,35 @@
+##=============================================================================
+##
+##  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.
+##
+##  Copyright 2015 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
+##  Copyright 2015 UT-Battelle, LLC.
+##  Copyright 2015 Los Alamos National Security.
+##
+##  Under the terms of Contract DE-NA0003525 with NTESS,
+##  the U.S. Government retains certain rights in this software.
+##  Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
+##  Laboratory (LANL), the U.S. Government retains certain rights in
+##  this software.
+##
+##=============================================================================
+cmake_minimum_required(VERSION 3.3...3.12 FATAL_ERROR)
+project(HelloWorld CXX)
+
+#Find the VTK-m package
+find_package(VTKm REQUIRED QUIET)
+vtkm_find_gl(OPTIONAL GL GLUT GLEW)
+
+#if(TARGET OpenGL::GL AND
+#        TARGET GLUT::GLUT AND
+#        TARGET GLEW::GLEW)
+    add_executable(moments_SERIAL moments.cxx)
+    target_link_libraries(moments_SERIAL
+            PRIVATE vtkm_cont ${gl_libs})
+#endif()

examples/moments/moments.cxx

@@ -0,0 +1,151 @@
+//
+// Created by ollie on 10/10/18.
+//
+
+#include <vtkm/cont/ArrayHandle.h>
+#include <vtkm/io/reader/VTKDataSetReader.h>
+
+#include <vtkm/worklet/DispatcherPointNeighborhood.h>
+#include <vtkm/worklet/WorkletPointNeighborhood.h>
+
+#include <iostream>
+#include <vtkm/io/writer/VTKDataSetWriter.h>
+
+namespace detail
+{
+
+template <int... factors>
+vtkm::Float64 eval_functor(vtkm::Float64 x, vtkm::Float64 y);
+
+template <int factor, int... factors>
+vtkm::Float64 eval_functor_impl(vtkm::Float64 x, vtkm::Float64 y)
+{
+  return eval_functor<factor>(x, y) * eval_functor<factors...>(x, y);
+}
+
+template <int... factors>
+vtkm::Float64 eval_functor(vtkm::Float64 x, vtkm::Float64 y)
+{
+  return eval_functor_impl<factors...>(x, y);
+}
+
+template <>
+vtkm::Float64 eval_functor<0>(vtkm::Float64 x, vtkm::Float64 y)
+{
+  return x;
+}
+
+template <>
+vtkm::Float64 eval_functor<1>(vtkm::Float64 x, vtkm::Float64 y)
+{
+  return y;
+}
+
+template <>
+vtkm::Float64 eval_functor<>(vtkm::Float64 x, vtkm::Float64 y)
+{
+  return 1;
+}
+
+} // detail
+
+struct Monomial
+{
+public:
+  Monomial(int _p, int _q)
+    : p(_p)
+    , q(_q)
+  {
+  }
+
+  VTKM_EXEC
+  vtkm::Float64 eval(vtkm::Float64 x, vtkm::Float64 y) const
+  {
+    // return detail::eval_functor<factors...>(x, y);
+    return pow(x, p) * pow(y, q);
+  }
+
+private:
+  int p;
+  int q;
+};
+
+// FIXME: Why can't we pass the template parameter Radius to WorkletPointNeighborhood?
+// TODO: figure out how to make Radius runtime configurable. We need to update
+// VTKM's implementation of WorkletPointNeighborhood
+template <int Radius>
+struct Moments : public vtkm::worklet::WorkletPointNeighborhood<2>
+{
+public:
+  Moments(int _p, int _q)
+    : monomial(_p, _q)
+  {
+    assert(_p >= 0);
+    assert(_q >= 0);
+  }
+
+  using ControlSignature = void(CellSetIn, FieldInNeighborhood<>, FieldOut<>);
+
+  using ExecutionSignature = void(_2, _3);
+
+  template <typename NeighIn>
+  VTKM_EXEC void operator()(const NeighIn& image, vtkm::Float64& moment) const
+  {
+    vtkm::Float64 sum = 0;
+    vtkm::Float64 recp = 1.0 / Radius;
+
+    for (int j = -Radius; j <= Radius; j++)
+    {
+      for (int i = -Radius; i <= Radius; i++)
+      {
+        if (i * i + j * j <= Radius * Radius)
+          // TODO: there is several order of magnitude of difference between this and the
+          // the reference implementation. Scale by dx dy?
+          sum += monomial.eval(i * recp, j * recp) * image.Get(i, j, 0);
+      }
+    }
+
+    moment = sum;
+  }
+
+private:
+  Monomial monomial;
+};
+
+int main(int argc, char* argv[])
+{
+  vtkm::io::reader::VTKDataSetReader reader(argv[1]);
+  vtkm::cont::DataSet input = reader.ReadDataSet();
+
+  auto field = input.GetPointField("values");
+  vtkm::cont::ArrayHandle<vtkm::Float64> pixels;
+  field.GetData().CopyTo(pixels);
+
+  vtkm::cont::DataSet output;
+  output.AddCellSet(input.GetCellSet(0));
+  output.AddCoordinateSystem(input.GetCoordinateSystem(0));
+
+  for (int order = 0; order <= 2; ++order)
+  {
+    for (int p = order; p >= 0; --p)
+    {
+      std::cout << "(" << p << ", " << order - p << ")" << std::endl;
+
+      vtkm::cont::ArrayHandle<vtkm::Float64> moments;
+
+      using DispatcherType = vtkm::worklet::DispatcherPointNeighborhood<Moments<2>>;
+      DispatcherType dispatcher(Moments<2>{ p, order - p });
+      dispatcher.SetDevice(vtkm::cont::DeviceAdapterTagSerial());
+      dispatcher.Invoke(input.GetCellSet(0), pixels, moments);
+
+      std::string fieldName = "moments_";
+      fieldName += std::to_string(p) + std::to_string(order - p);
+
+      vtkm::cont::Field momentsField(fieldName, vtkm::cont::Field::Association::POINTS, moments);
+      output.AddField(momentsField);
+    }
+  }
+
+  vtkm::io::writer::VTKDataSetWriter writer("moments_all.vtk");
+  writer.WriteDataSet(output);
+}