Add ComputeMoments worklet and filter

This commit is contained in:
Li-Ta Lo 2019-04-23 08:57:35 -06:00 committed by Allison Vacanti
parent a1b1bcd7f6
commit 434f751de1
5 changed files with 232 additions and 0 deletions

@ -0,0 +1,48 @@
//
// Created by ollie on 4/18/19.
//
#ifndef vtk_m_filter_ComputeMoments_h
#define vtk_m_filter_ComputeMoments_h
#include <vtkm/filter/FilterCell.h>
namespace vtkm
{
namespace filter
{
class ComputeMoments : public vtkm::filter::FilterCell<ComputeMoments>
{
public:
VTKM_CONT ComputeMoments();
template <typename T, typename StorageType, typename DerivedPolicy>
VTKM_CONT vtkm::cont::DataSet DoExecute(const vtkm::cont::DataSet& input,
const vtkm::cont::ArrayHandle<T, StorageType>& field,
const vtkm::filter::FieldMetadata& fieldMetadata,
const vtkm::filter::PolicyBase<DerivedPolicy>&);
VTKM_CONT void SetRadius(int _radius) { this->radius = _radius; }
VTKM_CONT void SetOrder(int _order) { this->order = _order; }
private:
int radius;
int order;
};
template <>
class FilterTraits<vtkm::filter::ComputeMoments>
{
public:
struct InputFieldTypeList : vtkm::TypeListTagScalarAll
{
};
};
}
} // namespace vtkm::filter
#include <vtkm/filter/ComputeMoments.hxx>
#endif //vtk_m_filter_ComputeMoments_h

@ -0,0 +1,78 @@
//=============================================================================
//
// 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 2018 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
// Copyright 2018 UT-Battelle, LLC.
// Copyright 2018 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.
//
//=============================================================================
#ifndef vtk_m_filter_ComputeMoments_hxx
#define vtk_m_filter_ComputeMoments_hxx
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/filter/internal/CreateResult.h>
#include <vtkm/worklet/moments/ComputeMoments.h>
namespace vtkm
{
namespace filter
{
VTKM_CONT ComputeMoments::ComputeMoments()
{
this->SetOutputFieldName("moments_");
}
template <typename T, typename StorageType, typename DerivedPolicy>
inline VTKM_CONT vtkm::cont::DataSet ComputeMoments::DoExecute(
const vtkm::cont::DataSet& input,
const vtkm::cont::ArrayHandle<T, StorageType>& field,
const vtkm::filter::FieldMetadata& fieldMetadata,
const vtkm::filter::PolicyBase<DerivedPolicy>&)
{
if (fieldMetadata.GetAssociation() != vtkm::cont::Field::Association::POINTS)
{
throw vtkm::cont::ErrorBadValue("Active field for ComputeMoments must be a point field.");
}
vtkm::cont::DataSet output = internal::CreateResult(input);
for (int i = 0; i <= this->order; ++i)
{
for (int p = i; p >= 0; --p)
{
vtkm::cont::ArrayHandle<T> moments;
using DispatcherType =
vtkm::worklet::DispatcherPointNeighborhood<vtkm::worklet::moments::ComputeMoments>;
DispatcherType dispatcher(vtkm::worklet::moments::ComputeMoments{ this->radius, p, i - p });
dispatcher.SetDevice(vtkm::cont::DeviceAdapterTagSerial());
dispatcher.Invoke(input.GetCellSet(0), field, moments);
std::string fieldName = "moments_";
fieldName += std::to_string(p) + std::to_string(i - p);
vtkm::cont::Field momentsField(fieldName, vtkm::cont::Field::Association::POINTS, moments);
output.AddField(momentsField);
}
}
return output;
}
}
}
#endif //vtk_m_filter_ComputeMoments_hxx

@ -127,6 +127,7 @@ add_subdirectory(cosmotools)
add_subdirectory(gradient)
add_subdirectory(histogram)
add_subdirectory(lcs)
add_subdirectory(moments)
add_subdirectory(splatkernels)
add_subdirectory(spatialstructure)
add_subdirectory(tetrahedralize)

@ -0,0 +1,26 @@
##============================================================================
## 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 2018 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
## Copyright 2018 UT-Battelle, LLC.
## Copyright 2018 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.
##============================================================================
set(headers
ComputeMoments.h
)
#-----------------------------------------------------------------------------
vtkm_declare_headers(${headers})

@ -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.
//
// Copyright 2018 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
// Copyright 2018 UT-Battelle, LLC.
// Copyright 2018 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.
//
//=============================================================================
#ifndef vtk_m_worklet_moments_ComputeMoments_h
#define vtk_m_worklet_moments_ComputeMoments_h
#include <vtkm/worklet/DispatcherPointNeighborhood.h>
#include <vtkm/worklet/WorkletPointNeighborhood.h>
namespace vtkm
{
namespace worklet
{
namespace moments
{
struct ComputeMoments : public vtkm::worklet::WorkletPointNeighborhood
{
public:
ComputeMoments(int _r, int _p, int _q)
: Radius(_r)
, p(_p)
, q(_q)
{
assert(_r >= 0);
assert(_p >= 0);
assert(_q >= 0);
}
using ControlSignature = void(CellSetIn, FieldInNeighborhood, FieldOut);
using ExecutionSignature = void(_2, _3);
template <typename NeighIn, typename T>
VTKM_EXEC void operator()(const NeighIn& image, T& 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)
sum += pow(i * recp, p) * pow(j * recp, q) * image.Get(i, j, 0);
}
}
moment = T(sum * recp * recp);
}
private:
const int Radius;
const int p;
const int q;
};
}
}
}
#endif // vtk_m_worklet_moments_ComputeMoments_h