Adding ScatterCounting

This commit is contained in:
Kenneth Moreland 2015-09-24 16:38:42 -06:00
parent 7b6e6e4a66
commit e6a9c96c96
7 changed files with 570 additions and 29 deletions

@ -28,10 +28,6 @@
#include <vtkm/TypeTraits.h>
#include <vtkm/VecTraits.h>
VTKM_THIRDPARTY_PRE_INCLUDE
#include <boost/static_assert.hpp>
VTKM_THIRDPARTY_POST_INCLUDE
#include <exception>
#include <iostream>
#include <sstream>

@ -33,6 +33,7 @@ set(headers
IsosurfaceUniformGrid.h
MarchingCubesDataTables.h
PointElevation.h
ScatterCounting.h
ScatterIdentity.h
TetrahedralizeExplicitGrid.h
TetrahedralizeUniformGrid.h

@ -0,0 +1,279 @@
//=============================================================================
//
// 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 Sandia Corporation.
// Copyright 2015 UT-Battelle, LLC.
// Copyright 2015 Los Alamos National Security.
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// 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_ScatterCounting_h
#define vtk_m_worklet_ScatterCounting_h
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/ArrayHandleCast.h>
#include <vtkm/cont/ArrayHandleIndex.h>
#include <vtkm/cont/DeviceAdapterAlgorithm.h>
#include <vtkm/cont/ErrorControlBadValue.h>
#include <vtkm/exec/FunctorBase.h>
#include <sstream>
namespace vtkm {
namespace worklet {
namespace detail {
template<typename Device>
struct ReverseInputToOutputMapKernel : vtkm::exec::FunctorBase
{
typedef typename
vtkm::cont::ArrayHandle<vtkm::Id>::ExecutionTypes<Device>::PortalConst
InputMapType;
typedef typename
vtkm::cont::ArrayHandle<vtkm::Id>::ExecutionTypes<Device>::Portal
OutputMapType;
InputMapType InputToOutputMap;
OutputMapType OutputToInputMap;
vtkm::Id OutputSize;
VTKM_CONT_EXPORT
ReverseInputToOutputMapKernel(const InputMapType &inputToOutputMap,
const OutputMapType &outputToInputMap,
vtkm::Id outputSize)
: InputToOutputMap(inputToOutputMap),
OutputToInputMap(outputToInputMap),
OutputSize(outputSize)
{ }
VTKM_EXEC_EXPORT
void operator()(vtkm::Id inputIndex) const
{
vtkm::Id outputStartIndex;
if (inputIndex > 0)
{
outputStartIndex = this->InputToOutputMap.Get(inputIndex-1);
}
else
{
outputStartIndex = 0;
}
vtkm::Id outputEndIndex = this->InputToOutputMap.Get(inputIndex);
for (vtkm::Id outputIndex = outputStartIndex;
outputIndex < outputEndIndex;
outputIndex++)
{
this->OutputToInputMap.Set(outputIndex, inputIndex);
}
}
};
template<typename Device>
struct SubtractToVisitIndexKernel : vtkm::exec::FunctorBase
{
typedef typename
vtkm::cont::ArrayHandle<vtkm::Id>::ExecutionTypes<Device>::PortalConst
StartsOfGroupsType;
typedef typename
vtkm::cont::ArrayHandle<vtkm::IdComponent>::ExecutionTypes<Device>::Portal
VisitType;
StartsOfGroupsType StartsOfGroups;
VisitType Visit;
VTKM_CONT_EXPORT
SubtractToVisitIndexKernel(const StartsOfGroupsType &startsOfGroups,
const VisitType &visit)
: StartsOfGroups(startsOfGroups), Visit(visit)
{ }
VTKM_EXEC_EXPORT
void operator()(vtkm::Id inputIndex) const
{
vtkm::Id startOfGroup = this->StartsOfGroups.Get(inputIndex);
vtkm::IdComponent visitIndex =
static_cast<vtkm::IdComponent>(inputIndex - startOfGroup);
this->Visit.Set(inputIndex, visitIndex);
}
};
} // namespace detail
/// \brief A scatter that maps input to some numbers of output.
///
/// The \c Scatter classes are responsible for defining how much output is
/// generated based on some sized input. \c ScatterCounting establishes a 1 to
/// N mapping from input to output. That is, every input element generates 0 or
/// more output elements associated with it. The output elements are grouped by
/// the input associated.
///
/// A counting scatter takes an array of counts for each input. The data is
/// taken in the constructor and the index arrays are derived from that. So
/// changing the counts after the scatter is created will have no effect.
///
struct ScatterCounting
{
template<typename CountArrayType, typename Device>
VTKM_CONT_EXPORT
ScatterCounting(const CountArrayType &countArray, Device)
{
this->BuildArrays(countArray, Device());
}
typedef vtkm::cont::ArrayHandle<vtkm::Id> OutputToInputMapType;
template<typename RangeType>
VTKM_CONT_EXPORT
OutputToInputMapType GetOutputToInputMap(RangeType) const
{
return this->OutputToInputMap;
}
typedef vtkm::cont::ArrayHandle<vtkm::IdComponent> VisitArrayType;
template<typename RangeType>
VTKM_CONT_EXPORT
VisitArrayType GetVisitArray(RangeType) const
{
return this->VisitArray;
}
VTKM_CONT_EXPORT
vtkm::Id GetOutputRange(vtkm::Id inputRange) const
{
if (inputRange != this->InputRange)
{
std::stringstream msg;
msg << "ScatterCounting initialized with input domain of size "
<< this->InputRange
<< " but used with a worklet invoke of size "
<< inputRange << std::endl;
throw vtkm::cont::ErrorControlBadValue(msg.str());
}
return this->VisitArray.GetNumberOfValues();
}
VTKM_CONT_EXPORT
vtkm::Id GetOutputRange(vtkm::Id3 inputRange) const
{
return this->GetOutputRange(inputRange[0]*inputRange[1]*inputRange[2]);
}
private:
vtkm::Id InputRange;
OutputToInputMapType OutputToInputMap;
VisitArrayType VisitArray;
template<typename CountArrayType, typename Device>
VTKM_CONT_EXPORT
void BuildArrays(const CountArrayType &count, Device)
{
VTKM_IS_ARRAY_HANDLE(CountArrayType);
VTKM_IS_DEVICE_ADAPTER_TAG(Device);
this->InputRange = count.GetNumberOfValues();
// Currently we are treating the input to output map as a temporary
// variable. However, it is possible that this could, be useful elsewhere,
// so we may want to save this and make it available.
//
// The input to output map is actually built off by one. The first entry
// is actually for the second value. The last entry is the total number of
// output. This off-by-one is so that an upper bound find will work when
// building the output to input map.
vtkm::cont::ArrayHandle<vtkm::Id> inputToOutputMap;
vtkm::Id outputSize =
vtkm::cont::DeviceAdapterAlgorithm<Device>::ScanInclusive(
vtkm::cont::make_ArrayHandleCast(count, vtkm::Id()),
inputToOutputMap);
// We have implemented two different ways to compute the output to input
// map. The first way is to use a binary search on each output index into
// the input map. The second way is to schedule on each input and
// iteratively fill all the output indices for that input. The first way is
// faster for output sizes that are small relative to the input (typical in
// Marching Cubes, for example) and also tends to be well load balanced.
// The second way is faster for larger outputs (typical in triangulation,
// for example). We will use the first method for small output sizes and
// the second for large output sizes. Toying with this might be a good
// place for optimization.
if (outputSize < this->InputRange)
{
this->BuildOutputToInputMapWithFind(
outputSize, inputToOutputMap, Device());
}
else
{
this->BuildOutputToInputMapWithIterate(
outputSize, inputToOutputMap, Device());
}
// This builds the visit indices using a parallel find. A prefix sum by
// key could be more efficient, but that is not implemented in the device
// adapter at the time of this writing.
this->BuildVisitArrayWithFind(Device());
}
template<typename Device>
VTKM_CONT_EXPORT
void BuildOutputToInputMapWithFind(
vtkm::Id outputSize,
vtkm::cont::ArrayHandle<vtkm::Id> inputToOutputMap,
Device)
{
vtkm::cont::ArrayHandleIndex outputIndices(outputSize);
vtkm::cont::DeviceAdapterAlgorithm<Device>::UpperBounds(
inputToOutputMap, outputIndices, this->OutputToInputMap);
}
template<typename Device>
VTKM_CONT_EXPORT
void BuildOutputToInputMapWithIterate(
vtkm::Id outputSize,
vtkm::cont::ArrayHandle<vtkm::Id> inputToOutputMap,
Device)
{
detail::ReverseInputToOutputMapKernel<Device>
kernel(inputToOutputMap.PrepareForInput(Device()),
this->OutputToInputMap.PrepareForOutput(outputSize, Device()),
outputSize);
vtkm::cont::DeviceAdapterAlgorithm<Device>::Schedule(
kernel, inputToOutputMap.GetNumberOfValues());
}
template<typename Device>
VTKM_CONT_EXPORT
void BuildVisitArrayWithFind(Device)
{
vtkm::cont::ArrayHandle<vtkm::Id> startsOfGroups;
// This find gives the index of the start of a group.
vtkm::cont::DeviceAdapterAlgorithm<Device>::LowerBounds(
this->OutputToInputMap, this->OutputToInputMap, startsOfGroups);
vtkm::Id outputSize = this->OutputToInputMap.GetNumberOfValues();
detail::SubtractToVisitIndexKernel<Device>
kernel(startsOfGroups.PrepareForInput(Device()),
this->VisitArray.PrepareForOutput(outputSize, Device()));
vtkm::cont::DeviceAdapterAlgorithm<Device>::Schedule(kernel, outputSize);
}
};
}
} // namespace vtkm::worklet
#endif //vtk_m_worklet_ScatterCounting_h

@ -32,41 +32,40 @@ namespace worklet {
///
/// The \c Scatter classes are responsible for defining how much output is
/// generated based on some sized input. \c ScatterIdentity establishes a 1 to
/// 1 mapping from output to input (and vice versa). That is, every input
/// 1 mapping from input to output (and vice versa). That is, every input
/// element generates one output element associated with it. This is the
/// default for basic maps.
///
struct ScatterIdentity
{
typedef vtkm::cont::ArrayHandleIndex OutputToInputMapType;
OutputToInputMapType GetOutputToInputMap(vtkm::Id inputRange) const
VTKM_CONT_EXPORT
OutputToInputMapType GetOutputToInputMap(vtkm::Id outputRange) const
{
return OutputToInputMapType(inputRange);
return OutputToInputMapType(outputRange);
}
OutputToInputMapType GetOutputToInputMap(vtkm::Id2 inputRange) const
VTKM_CONT_EXPORT
OutputToInputMapType GetOutputToInputMap(vtkm::Id3 outputRange) const
{
return this->GetOutputToInputMap(inputRange[0]*inputRange[1]);
}
OutputToInputMapType GetOutputToInputMap(vtkm::Id3 inputRange) const
{
return this->GetOutputToInputMap(inputRange[0]*inputRange[1]*inputRange[2]);
return this->GetOutputToInputMap(
outputRange[0]*outputRange[1]*outputRange[2]);
}
typedef vtkm::cont::ArrayHandleConstant<vtkm::IdComponent> VisitArrayType;
VisitArrayType GetVisitArray(vtkm::Id inputRange) const
VTKM_CONT_EXPORT
VisitArrayType GetVisitArray(vtkm::Id outputRange) const
{
return VisitArrayType(1, inputRange);
return VisitArrayType(1, outputRange);
}
VisitArrayType GetVisitArray(vtkm::Id2 inputRange) const
VTKM_CONT_EXPORT
VisitArrayType GetVisitArray(vtkm::Id3 outputRange) const
{
return this->GetVisitArray(inputRange[0]*inputRange[1]);
}
VisitArrayType GetVisitArray(vtkm::Id3 inputRange) const
{
return this->GetVisitArray(inputRange[0]*inputRange[1]*inputRange[2]);
return this->GetVisitArray(outputRange[0]*outputRange[1]*outputRange[2]);
}
vtkm::Id GetOutputRange(vtkm::Id inputRange) const
template<typename RangeType>
VTKM_CONT_EXPORT
RangeType GetOutputRange(RangeType inputRange) const
{
return inputRange;
}

@ -441,19 +441,23 @@ protected:
VTKM_CONT_EXPORT
void BasicInvoke(const Invocation &invocation,
vtkm::Id numInstances,
DeviceAdapter tag) const
DeviceAdapter device) const
{
this->InvokeTransportParameters(invocation, numInstances, tag);
this->InvokeTransportParameters(
invocation,
this->Worklet.GetScatter().GetOutputRange(numInstances),
device);
}
template<typename Invocation, typename DeviceAdapter>
VTKM_CONT_EXPORT
void BasicInvoke(const Invocation &invocation,
vtkm::Id2 dimensions,
DeviceAdapter tag) const
DeviceAdapter device) const
{
vtkm::Id3 dim3d(dimensions[0], dimensions[1], 1);
this->InvokeTransportParameters(invocation, dim3d, tag);
this->BasicInvoke(invocation,
vtkm::Id3(dimensions[0], dimensions[1], 1),
device);
}
@ -461,9 +465,12 @@ protected:
VTKM_CONT_EXPORT
void BasicInvoke(const Invocation &invocation,
vtkm::Id3 dimensions,
DeviceAdapter tag) const
DeviceAdapter device) const
{
this->InvokeTransportParameters(invocation, dimensions, tag);
this->InvokeTransportParameters(
invocation,
this->Worklet.GetScatter().GetOutputRange(dimensions),
device);
}
WorkletType Worklet;

@ -26,6 +26,7 @@ set(unit_tests
UnitTestFieldStatistics.cxx
UnitTestIsosurfaceUniformGrid.cxx
UnitTestPointElevation.cxx
UnitTestScatterCounting.cxx
UnitTestSplatKernels.cxx
UnitTestTetrahedralizeExplicitGrid.cxx
UnitTestTetrahedralizeUniformGrid.cxx

@ -0,0 +1,258 @@
//=============================================================================
//
// 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 Sandia Corporation.
// Copyright 2015 UT-Battelle, LLC.
// Copyright 2015 Los Alamos National Security.
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// 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.
//
//=============================================================================
#include <vtkm/worklet/ScatterCounting.h>
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/DeviceAdapterAlgorithm.h>
#include <vtkm/worklet/DispatcherMapField.h>
#include <vtkm/worklet/WorkletMapField.h>
#include <vtkm/cont/testing/Testing.h>
#include <vector>
namespace {
struct TestScatterArrays
{
vtkm::cont::ArrayHandle<vtkm::IdComponent> CountArray;
vtkm::cont::ArrayHandle<vtkm::Id> OutputToInputMap;
vtkm::cont::ArrayHandle<vtkm::IdComponent> VisitArray;
};
TestScatterArrays MakeScatterArraysShort()
{
const vtkm::Id countArraySize = 18;
const vtkm::IdComponent countArray[countArraySize] = {
1, 2, 0, 0, 1, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0
};
const vtkm::Id outputSize = 6;
const vtkm::Id outputToInputMap[outputSize] = {
0, 1, 1, 4, 6, 14
};
const vtkm::IdComponent visitArray[outputSize] = {
0, 0, 1, 0, 0, 0
};
TestScatterArrays arrays;
typedef vtkm::cont::DeviceAdapterAlgorithm<VTKM_DEFAULT_DEVICE_ADAPTER_TAG>
Algorithm;
// Need to copy arrays so that the data does not go out of scope.
Algorithm::Copy(vtkm::cont::make_ArrayHandle(countArray, countArraySize),
arrays.CountArray);
Algorithm::Copy(vtkm::cont::make_ArrayHandle(outputToInputMap, outputSize),
arrays.OutputToInputMap);
Algorithm::Copy(vtkm::cont::make_ArrayHandle(visitArray, outputSize),
arrays.VisitArray);
return arrays;
}
TestScatterArrays MakeScatterArraysLong()
{
const vtkm::Id countArraySize = 6;
const vtkm::IdComponent countArray[countArraySize] = {
0, 1, 2, 3, 4, 5
};
const vtkm::Id outputSize = 15;
const vtkm::Id outputToInputMap[outputSize] = {
1, 2, 2, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5, 5
};
const vtkm::IdComponent visitArray[outputSize] = {
0, 0, 1, 0, 1, 2, 0, 1, 2, 3, 0, 1, 2, 3, 4
};
TestScatterArrays arrays;
typedef vtkm::cont::DeviceAdapterAlgorithm<VTKM_DEFAULT_DEVICE_ADAPTER_TAG>
Algorithm;
// Need to copy arrays so that the data does not go out of scope.
Algorithm::Copy(vtkm::cont::make_ArrayHandle(countArray, countArraySize),
arrays.CountArray);
Algorithm::Copy(vtkm::cont::make_ArrayHandle(outputToInputMap, outputSize),
arrays.OutputToInputMap);
Algorithm::Copy(vtkm::cont::make_ArrayHandle(visitArray, outputSize),
arrays.VisitArray);
return arrays;
}
TestScatterArrays MakeScatterArraysZero()
{
const vtkm::Id countArraySize = 6;
const vtkm::IdComponent countArray[countArraySize] = {
0, 0, 0, 0, 0, 0
};
TestScatterArrays arrays;
typedef vtkm::cont::DeviceAdapterAlgorithm<VTKM_DEFAULT_DEVICE_ADAPTER_TAG>
Algorithm;
// Need to copy arrays so that the data does not go out of scope.
Algorithm::Copy(vtkm::cont::make_ArrayHandle(countArray, countArraySize),
arrays.CountArray);
arrays.OutputToInputMap.Allocate(0);
arrays.VisitArray.Allocate(0);
return arrays;
}
struct TestScatterCountingWorklet : public vtkm::worklet::WorkletMapField
{
typedef void ControlSignature(FieldIn<> inputIndices,
FieldOut<> copyIndices,
FieldOut<> recordVisit,
FieldOut<> recordWorkId);
typedef void ExecutionSignature(_1, _2 ,_3, _4, VisitIndex, WorkIndex);
typedef vtkm::worklet::ScatterCounting ScatterType;
VTKM_CONT_EXPORT
ScatterType GetScatter() const { return this->Scatter; }
template<typename CountArrayType>
VTKM_CONT_EXPORT
TestScatterCountingWorklet(const CountArrayType &countArray)
: Scatter(countArray, VTKM_DEFAULT_DEVICE_ADAPTER_TAG()) { }
template<typename CountArrayType, typename Device>
VTKM_CONT_EXPORT
TestScatterCountingWorklet(const CountArrayType &countArray, Device)
: Scatter(countArray, Device()) { }
VTKM_CONT_EXPORT
TestScatterCountingWorklet(const vtkm::worklet::ScatterCounting &scatter)
: Scatter(scatter) { }
VTKM_EXEC_EXPORT
void operator()(vtkm::Id inputIndex,
vtkm::Id &indexCopy,
vtkm::IdComponent &writeVisit,
vtkm::Float32 &captureWorkId,
vtkm::IdComponent visitIndex,
vtkm::Id workId) const
{
indexCopy = inputIndex;
writeVisit = visitIndex;
captureWorkId = TestValue(workId, vtkm::Float32());
}
private:
ScatterType Scatter;
};
template<typename T>
void CompareArrays(vtkm::cont::ArrayHandle<T> array1,
vtkm::cont::ArrayHandle<T> array2)
{
typedef typename vtkm::cont::ArrayHandle<T>::PortalConstControl PortalType;
PortalType portal1 = array1.GetPortalConstControl();
PortalType portal2 = array2.GetPortalConstControl();
VTKM_TEST_ASSERT(portal1.GetNumberOfValues() == portal2.GetNumberOfValues(),
"Arrays are not the same length.");
for (vtkm::Id index = 0; index < portal1.GetNumberOfValues(); index++)
{
T value1 = portal1.Get(index);
T value2 = portal2.Get(index);
VTKM_TEST_ASSERT(value1 == value2, "Array values not equal.");
}
}
// This unit test makes sure the ScatterCounting generates the correct map
// and visit arrays.
void TestScatterArrayGeneration(const TestScatterArrays &arrays)
{
std::cout << " Testing array generation" << std::endl;
vtkm::worklet::ScatterCounting scatter(arrays.CountArray,
VTKM_DEFAULT_DEVICE_ADAPTER_TAG());
vtkm::Id inputSize = arrays.CountArray.GetNumberOfValues();
std::cout << " Checking output to input map." << std::endl;
CompareArrays(arrays.OutputToInputMap,
scatter.GetOutputToInputMap(inputSize));
std::cout << " Checking visit array." << std::endl;
CompareArrays(arrays.VisitArray,
scatter.GetVisitArray(inputSize));
}
// This is more of an integration test that makes sure the scatter works with a
// worklet invocation.
void TestScatterWorklet(const TestScatterArrays &arrays)
{
std::cout << " Testing scatter counting in a worklet." << std::endl;
vtkm::worklet::ScatterCounting scatter(arrays.CountArray,
VTKM_DEFAULT_DEVICE_ADAPTER_TAG());
TestScatterCountingWorklet worklet(scatter);
vtkm::worklet::DispatcherMapField<TestScatterCountingWorklet> dispatcher(
worklet);
vtkm::Id inputSize = arrays.CountArray.GetNumberOfValues();
vtkm::cont::ArrayHandleIndex inputIndices(inputSize);
vtkm::cont::ArrayHandle<vtkm::Id> outputToInputMapCopy;
vtkm::cont::ArrayHandle<vtkm::IdComponent> visitCopy;
vtkm::cont::ArrayHandle<vtkm::Float32> captureWorkId;
std::cout << " Invoke worklet" << std::endl;
dispatcher.Invoke(
inputIndices, outputToInputMapCopy, visitCopy, captureWorkId);
std::cout << " Check output to input map." << std::endl;
CompareArrays(outputToInputMapCopy, arrays.OutputToInputMap);
std::cout << " Check visit." << std::endl;
CompareArrays(visitCopy, arrays.VisitArray);
std::cout << " Check work id." << std::endl;
CheckPortal(captureWorkId.GetPortalConstControl());
}
void TestScatterCountingWithArrays(const TestScatterArrays &arrays)
{
TestScatterArrayGeneration(arrays);
TestScatterWorklet(arrays);
}
void TestScatterCounting()
{
std::cout << "Testing arrays with output smaller than input." << std::endl;
TestScatterCountingWithArrays(MakeScatterArraysShort());
std::cout << "Testing arrays with output larger than input." << std::endl;
TestScatterCountingWithArrays(MakeScatterArraysLong());
std::cout << "Testing arrays with zero output." << std::endl;
TestScatterCountingWithArrays(MakeScatterArraysZero());
}
} // anonymous namespace
int UnitTestScatterCounting(int, char *[])
{
return vtkm::cont::testing::Testing::Run(TestScatterCounting);
}