Shuenhoy/vtk-m
1
0
Fork 0
mirror of https://gitlab.kitware.com/vtk/vtk-m synced 2024-09-16 17:22:55 +00:00
Code Issues 2 Actions Packages Projects Releases Wiki Activity

Merge topic 'marching_cubes_normal_gen_in_post_pass'

Browse Source 
881a6b1c Correct warnings in ImplicitFunction when using std::copy
0e31418f Improve the quality of normals from MarchingCubes.
ea731728 Refactor MarchingCubes normal pass to use the ReduceWorklet.
b241299d MarchingCubes now generates normals in a post pass.
6ed4bc78 Permuted structured cellsets produce VecRectilinearPointCoordinates
6da48cf3 Remove unneeded methods from the Connectivity classes.
b56f1604 Add CastAndCall specializations for the concrete CellSet types.
9d75e7b7 Remove unneeded member variables from tbb ScheduleKernelId3

Acked-by: Kitware Robot <kwrobot@kitware.com>
Merge-request: !741
This commit is contained in:
Robert Maynard 2017-05-23 16:23:46 +00:00 committed by Kitware Robot
commit 40a0b5dcce
13 changed files with 578 additions and 319 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
vtkm/cont/ImplicitFunction.h
View File

@ -20,12 +20,13 @@
#ifndef vtk_m_cont_ImplicitFunction_h
#define vtk_m_cont_ImplicitFunction_h
#include <vtkm/internal/Configure.h>
#include <vtkm/cont/VirtualObjectCache.h>
#include <vtkm/exec/ImplicitFunction.h>
#include <memory>
namespace vtkm {
namespace cont {

2
vtkm/cont/VirtualObjectCache.h
View File

@ -20,10 +20,10 @@
#ifndef vtk_m_cont_VirtualObjectCache_h
#define vtk_m_cont_VirtualObjectCache_h
#include <vtkm/cont/DeviceAdapterListTag.h>
#include <vtkm/cont/ErrorBadType.h>
#include <vtkm/cont/ErrorBadValue.h>
#include <vtkm/cont/internal/DeviceAdapterTag.h>
#include <vtkm/cont/DeviceAdapterListTag.h>
#include <vtkm/cont/internal/VirtualObjectTransfer.h>
#include <array>

4
vtkm/cont/arg/testing/UnitTestTransportCellSetIn.cxx
View File

@ -45,8 +45,8 @@ struct TestKernel : public vtkm::exec::FunctorBase
this->RaiseError("Got bad number of shapes in exec cellset object.");
}
if (this->CellSet.GetNumberOfIndices(0) != 3 ||
this->CellSet.GetNumberOfIndices(1) != 4 )
if (this->CellSet.GetIndices(0).GetNumberOfComponents() != 3 ||
this->CellSet.GetIndices(1).GetNumberOfComponents() != 4 )
{
this->RaiseError("Got bad number of Indices in exec cellset object.");
}

41
vtkm/cont/internal/DynamicTransform.h
View File

@ -26,6 +26,10 @@ namespace vtkm {
namespace cont {
template<typename T, typename S> class ArrayHandle;
template<vtkm::IdComponent> class CellSetStructured;
template<typename T> class CellSetSingleType;
template<typename T, typename S, typename U, typename V> class CellSetExplicit;
template<typename T, typename S> class CellSetPermutation;
/// A Generic interface to CastAndCall. The default implementation simply calls
/// DynamicObject's CastAndCall, but specializations of this function exist for
@ -45,6 +49,43 @@ void CastAndCall(const vtkm::cont::ArrayHandle<T,U>& handle, const Functor &f)
f(handle);
}
/// A specialization of CastAndCall for basic CellSetStructured types,
/// Since the type is already known no deduction is needed.
/// This specialization is used to simplify numerous worklet algorithms
template<vtkm::IdComponent Dim, typename Functor>
void CastAndCall(const vtkm::cont::CellSetStructured<Dim>& cellset, const Functor &f)
{
f(cellset);
}
/// A specialization of CastAndCall for basic CellSetSingleType types,
/// Since the type is already known no deduction is needed.
/// This specialization is used to simplify numerous worklet algorithms
template<typename ConnectivityStorageTag, typename Functor>
void CastAndCall(const vtkm::cont::CellSetSingleType<ConnectivityStorageTag>& cellset, const Functor &f)
{
f(cellset);
}
/// A specialization of CastAndCall for basic CellSetExplicit types,
/// Since the type is already known no deduction is needed.
/// This specialization is used to simplify numerous worklet algorithms
template<typename T, typename S, typename U, typename V, typename Functor>
void CastAndCall(const vtkm::cont::CellSetExplicit<T,S,U,V>& cellset, const Functor &f)
{
f(cellset);
}
/// A specialization of CastAndCall for basic CellSetPermutation types,
/// Since the type is already known no deduction is needed.
/// This specialization is used to simplify numerous worklet algorithms
template<typename PermutationType, typename CellSetType, typename Functor>
void CastAndCall(const vtkm::cont::CellSetPermutation<PermutationType,CellSetType>& cellset, const Functor &f)
{
f(cellset);
}
namespace internal {
/// Tag used to identify an object that is a dynamic object that contains a

2
vtkm/cont/tbb/internal/DeviceAdapterAlgorithmTBB.h
View File

@ -152,7 +152,7 @@ public:
0, rangeMax[1], TBB_GRAIN_SIZE_3D[1],
0, rangeMax[0], TBB_GRAIN_SIZE_3D[2]);
tbb::ScheduleKernelId3<FunctorType> kernel(functor,rangeMax);
tbb::ScheduleKernelId3<FunctorType> kernel(functor);
kernel.SetErrorMessageBuffer(errorMessage);
::tbb::parallel_for(range, kernel);

30
vtkm/cont/tbb/internal/FunctorsTBB.h
View File

@ -508,10 +508,8 @@ template<class FunctorType>
class ScheduleKernelId3
{
public:
VTKM_CONT ScheduleKernelId3(const FunctorType &functor,
const vtkm::Id3& dims)
: Functor(functor),
Dims(dims)
VTKM_CONT ScheduleKernelId3(const FunctorType &functor)
: Functor(functor)
{ }
VTKM_CONT void SetErrorMessageBuffer(
@ -525,17 +523,24 @@ public:
void operator()(const ::tbb::blocked_range3d<vtkm::Id> &range) const {
try
{
for( vtkm::Id k=range.pages().begin(); k!=range.pages().end(); ++k)
const vtkm::Id kstart = range.pages().begin();
const vtkm::Id kend = range.pages().end();
const vtkm::Id jstart =range.rows().begin();
const vtkm::Id jend = range.rows().end();
const vtkm::Id istart =range.cols().begin();
const vtkm::Id iend = range.cols().end();
vtkm::Id3 index;
for( vtkm::Id k=kstart; k!=kend; ++k)
{
for( vtkm::Id j=range.rows().begin(); j!=range.rows().end(); ++j)
index[2]=k;
for( vtkm::Id j=jstart; j!=jend; ++j)
{
const vtkm::Id start =range.cols().begin();
const vtkm::Id end = range.cols().end();
VTKM_VECTORIZATION_PRE_LOOP
for( vtkm::Id i=start; i != end; ++i)
index[1]=j;
for( vtkm::Id i=istart; i != iend; ++i)
{
VTKM_VECTORIZATION_IN_LOOP
this->Functor(vtkm::Id3(i, j, k));
index[0]=i;
this->Functor(index);
}
}
}
@ -552,7 +557,6 @@ VTKM_VECTORIZATION_IN_LOOP
}
private:
FunctorType Functor;
vtkm::Id3 Dims;
vtkm::exec::internal::ErrorMessageBuffer ErrorMessage;
};

12
vtkm/exec/ConnectivityExplicit.h
View File

@ -36,6 +36,8 @@ template<typename ShapePortalType,
class ConnectivityExplicit
{
public:
typedef typename vtkm::Id SchedulingRangeType;
ConnectivityExplicit() {}
ConnectivityExplicit(const ShapePortalType& shapePortal,
@ -52,17 +54,11 @@ public:
}
VTKM_EXEC
vtkm::Id GetNumberOfElements() const
SchedulingRangeType GetNumberOfElements() const
{
return this->Shapes.GetNumberOfValues();
}
VTKM_EXEC
vtkm::IdComponent GetNumberOfIndices(vtkm::Id index) const
{
return static_cast<vtkm::IdComponent>(this->NumIndices.Get(index));
}
typedef vtkm::CellShapeTagGeneric CellShapeTag;
VTKM_EXEC
@ -82,7 +78,7 @@ public:
IndicesType GetIndices(vtkm::Id index) const
{
vtkm::Id offset = this->IndexOffset.Get(index);
vtkm::IdComponent length = this->GetNumberOfIndices(index);
vtkm::IdComponent length = this->NumIndices.Get(index);
return IndicesType(this->Connectivity, length, offset);
}

12
vtkm/exec/ConnectivityPermuted.h
View File

@ -34,7 +34,7 @@ template<typename PermutationPortal,
class ConnectivityPermuted
{
public:
typedef vtkm::Id SchedulingRangeType;
typedef typename OriginalConnectivity::SchedulingRangeType SchedulingRangeType;
VTKM_SUPPRESS_EXEC_WARNINGS
VTKM_EXEC_CONT
@ -60,12 +60,6 @@ public:
{
}
VTKM_EXEC
vtkm::IdComponent GetNumberOfIndices(vtkm::Id index) const {
return this->Connectivity.GetNumberOfIndices( this->Portal.Get(index) );
}
typedef typename OriginalConnectivity::CellShapeTag CellShapeTag;
VTKM_EXEC
@ -76,13 +70,13 @@ public:
typedef typename OriginalConnectivity::IndicesType IndicesType;
template<typename IndexType>
VTKM_EXEC
IndicesType GetIndices(vtkm::Id index) const
IndicesType GetIndices(const IndexType &index) const
{
return this->Connectivity.GetIndices( this->Portal.Get(index) );
}
private:
PermutationPortal Portal;
OriginalConnectivity Connectivity;
};

6
vtkm/exec/ConnectivityStructured.h
View File

@ -64,12 +64,6 @@ public:
{
}
template<typename IndexType>
VTKM_EXEC
vtkm::IdComponent GetNumberOfIndices(const IndexType &index) const {
return Helper::GetNumberOfIndices(this->Internals, index);
}
// This needs some thought. What does cell shape mean when the to topology
// is not a cell?
typedef typename InternalsType::CellShapeTag CellShapeTag;

36
vtkm/exec/arg/FetchTagArrayTopologyMapIn.h
View File

@ -167,6 +167,42 @@ struct FetchArrayTopologyMapInImplementation<
}
};
template<typename PermutationPortal, vtkm::IdComponent NumDimensions>
struct FetchArrayTopologyMapInImplementation<
vtkm::exec::ConnectivityPermuted<PermutationPortal,
vtkm::exec::ConnectivityStructured<vtkm::TopologyElementTagPoint,
vtkm::TopologyElementTagCell,
NumDimensions> >,
vtkm::internal::ArrayPortalUniformPointCoordinates>
{
typedef vtkm::exec::ConnectivityPermuted<PermutationPortal,
vtkm::exec::ConnectivityStructured<vtkm::TopologyElementTagPoint,
vtkm::TopologyElementTagCell,
NumDimensions> > ConnectivityType;
typedef vtkm::exec::arg::ThreadIndicesTopologyMap<ConnectivityType>
ThreadIndicesType;
typedef vtkm::VecRectilinearPointCoordinates<NumDimensions> ValueType;
VTKM_SUPPRESS_EXEC_WARNINGS
VTKM_EXEC
static ValueType Load(
const ThreadIndicesType &indices,
const vtkm::internal::ArrayPortalUniformPointCoordinates &field)
{
// This works because the logical cell index is the same as the logical
// point index of the first point on the cell.
// we have a flat index but we need 3d rectilinear coordiantes, so we
// need to take an flat index and convert to logical index
return vtkm::exec::arg::detail::make_VecRectilinearPointCoordinates(
field.GetOrigin(),
field.GetSpacing(),
indices.GetIndexLogical());
}
};
} // namespace detail
template<typename ConnectivityType, typename ExecObjectType>

195
vtkm/exec/arg/ThreadIndicesTopologyMap.h
View File

@ -23,6 +23,7 @@
#include <vtkm/exec/arg/ThreadIndicesBasic.h>
#include <vtkm/exec/ConnectivityStructured.h>
#include <vtkm/exec/ConnectivityPermuted.h>
namespace vtkm {
namespace exec {
@ -69,18 +70,18 @@ class ThreadIndicesTopologyMap : public vtkm::exec::arg::ThreadIndicesBasic
typedef vtkm::exec::arg::ThreadIndicesBasic Superclass;
public:
typedef typename ConnectivityType::IndicesType IndicesFromType;
typedef typename ConnectivityType::CellShapeTag CellShapeTag;
using IndicesFromType = typename ConnectivityType::IndicesType;
using CellShapeTag = typename ConnectivityType::CellShapeTag;
VTKM_SUPPRESS_EXEC_WARNINGS
template<typename OutToInArrayType, typename VisitArrayType>
VTKM_EXEC
ThreadIndicesTopologyMap(vtkm::Id threadIndex,
const OutToInArrayType& inToOut,
const OutToInArrayType& outToIn,
const VisitArrayType& visit,
const ConnectivityType& connectivity,
vtkm::Id globalThreadIndexOffset=0)
: Superclass(threadIndex, inToOut.Get(threadIndex), visit.Get(threadIndex),
: Superclass(threadIndex, outToIn.Get(threadIndex), visit.Get(threadIndex),
globalThreadIndexOffset),
CellShape(detail::CellShapeInitializer<CellShapeTag>::GetDefault())
{
@ -197,25 +198,23 @@ template<typename FromTopology,
class ThreadIndicesTopologyMap<
vtkm::exec::ConnectivityStructured<FromTopology,ToTopology,Dimension> >
{
typedef vtkm::exec::arg::ThreadIndicesBasic Superclass;
typedef vtkm::exec::ConnectivityStructured<FromTopology,ToTopology,Dimension>
ConnectivityType;
using ConnectivityType = vtkm::exec::ConnectivityStructured<FromTopology,ToTopology,Dimension>;
public:
typedef typename ConnectivityType::IndicesType IndicesFromType;
typedef typename ConnectivityType::CellShapeTag CellShapeTag;
typedef typename ConnectivityType::SchedulingRangeType LogicalIndexType;
using IndicesFromType = typename ConnectivityType::IndicesType;
using CellShapeTag = typename ConnectivityType::CellShapeTag;
using LogicalIndexType = typename ConnectivityType::SchedulingRangeType;
template<typename OutToInArrayType, typename VisitArrayType>
VTKM_EXEC
ThreadIndicesTopologyMap(vtkm::Id threadIndex,
const OutToInArrayType& inToOut,
const OutToInArrayType& outToIn,
const VisitArrayType& visit,
const ConnectivityType& connectivity,
vtkm::Id globalThreadIndexOffset=0)
{
this->InputIndex = inToOut.Get(threadIndex);
this->InputIndex = outToIn.Get(threadIndex);
this->OutputIndex = threadIndex;
this->VisitIndex = visit.Get(threadIndex);
this->LogicalIndex = connectivity.FlatToLogicalToIndex(this->InputIndex);
@ -373,6 +372,178 @@ private:
vtkm::Id GlobalThreadIndexOffset;
};
// Specialization for permuted structured connectivity types.
template<typename PermutationPortal,
typename FromTopology,
typename ToTopology,
vtkm::IdComponent Dimension>
class ThreadIndicesTopologyMap<
vtkm::exec::ConnectivityPermuted<PermutationPortal,
vtkm::exec::ConnectivityStructured<FromTopology,ToTopology,Dimension> >
>
{
using PermutedConnectivityType =
vtkm::exec::ConnectivityPermuted<PermutationPortal,
vtkm::exec::ConnectivityStructured<FromTopology,
ToTopology,
Dimension>
>;
using ConnectivityType = vtkm::exec::ConnectivityStructured<FromTopology,
ToTopology,
Dimension
>;
public:
using IndicesFromType = typename ConnectivityType::IndicesType;
using CellShapeTag = typename ConnectivityType::CellShapeTag;
using LogicalIndexType = typename ConnectivityType::SchedulingRangeType;
template<typename OutToInArrayType, typename VisitArrayType>
VTKM_EXEC
ThreadIndicesTopologyMap(vtkm::Id threadIndex,
const OutToInArrayType& outToIn,
const VisitArrayType& visit,
const PermutedConnectivityType& permutation,
vtkm::Id globalThreadIndexOffset=0)
{
this->InputIndex = outToIn.Get( threadIndex );
this->OutputIndex = threadIndex;
this->VisitIndex = visit.Get(threadIndex);
const vtkm::Id permutedIndex = permutation.Portal.Get( this->InputIndex );
this->LogicalIndex = permutation.Connectivity.FlatToLogicalToIndex(permutedIndex);
this->IndicesFrom = permutation.Connectivity.GetIndices(this->LogicalIndex);
this->CellShape = permutation.Connectivity.GetCellShape(permutedIndex);
this->GlobalThreadIndexOffset = globalThreadIndexOffset;
}
VTKM_SUPPRESS_EXEC_WARNINGS
VTKM_EXEC
ThreadIndicesTopologyMap(vtkm::Id threadIndex,
vtkm::Id vtkmNotUsed(inIndex),
vtkm::IdComponent visitIndex,
const PermutedConnectivityType& permutation,
vtkm::Id globalThreadIndexOffset=0)
{
this->InputIndex = threadIndex;
this->OutputIndex = threadIndex;
this->VisitIndex = visitIndex;
const vtkm::Id permutedIndex = permutation.Portal.Get( this->InputIndex );
this->LogicalIndex = permutation.Connectivity.FlatToLogicalToIndex(permutedIndex);
this->IndicesFrom = permutation.Connectivity.GetIndices(this->LogicalIndex);
this->CellShape = permutation.Connectivity.GetCellShape(permutedIndex);
this->GlobalThreadIndexOffset = globalThreadIndexOffset;
}
/// \brief The logical index into the input domain.
///
/// This is similar to \c GetIndex3D except the Vec size matches the actual
/// dimensions of the data.
///
VTKM_EXEC
LogicalIndexType GetIndexLogical() const
{
return this->LogicalIndex;
}
/// \brief The index into the input domain.
///
/// This index refers to the input element (array value, cell, etc.) that
/// this thread is being invoked for. This is the typical index used during
/// fetches.
///
VTKM_EXEC
vtkm::Id GetInputIndex() const
{
return this->InputIndex;
}
/// \brief The 3D index into the input domain.
///
/// Overloads the implementation in the base class to return the 3D index
/// for the input.
///
VTKM_EXEC
vtkm::Id3 GetInputIndex3D() const
{
return detail::InflateTo3D(this->GetIndexLogical());
}
/// \brief The index into the output domain.
///
/// This index refers to the output element (array value, cell, etc.) that
/// this thread is creating. This is the typical index used during
/// Fetch::Store.
///
VTKM_EXEC
vtkm::Id GetOutputIndex() const
{
return this->OutputIndex;
}
/// \brief The visit index.
///
/// When multiple output indices have the same input index, they are
/// distinguished using the visit index.
///
VTKM_EXEC
vtkm::IdComponent GetVisitIndex() const
{
return this->VisitIndex;
}
VTKM_EXEC
vtkm::Id GetGlobalIndex() const
{
return (this->GlobalThreadIndexOffset + this->OutputIndex);
}
/// \brief The input indices of the "from" elements.
///
/// A topology map has "from" and "to" elements (for example from points to
/// cells). For each worklet invocation, there is exactly one "to" element,
/// but can be several "from" element. This method returns a Vec-like object
/// containing the indices to the "from" elements.
///
VTKM_EXEC
const IndicesFromType &GetIndicesFrom() const { return this->IndicesFrom; }
/// \brief The input indices of the "from" elements in pointer form.
///
/// Returns the same object as GetIndicesFrom except that it returns a
/// pointer to the internally held object rather than a reference or copy.
/// Since the from indices can be a sizeable Vec (8 entries is common), it is
/// best not to have a bunch a copies. Thus, you can pass around a pointer
/// instead. However, care should be taken to make sure that this object does
/// not go out of scope, at which time the returned pointer becomes invalid.
///
VTKM_EXEC
const IndicesFromType *GetIndicesFromPointer() const
{
return &this->IndicesFrom;
}
/// \brief The shape of the input cell.
///
/// In topology maps that map from points to something, the indices make up
/// the structure of a cell. Although the shape tag is not technically and
/// index, it defines the meaning of the indices, so we put it here. (That
/// and this class is the only convenient place to store it.)
///
VTKM_EXEC
CellShapeTag GetCellShape() const { return this->CellShape; }
private:
vtkm::Id InputIndex;
vtkm::Id OutputIndex;
vtkm::IdComponent VisitIndex;
LogicalIndexType LogicalIndex;
IndicesFromType IndicesFrom;
CellShapeTag CellShape;
vtkm::Id GlobalThreadIndexOffset;
};
}
}
} // namespace vtkm::exec::arg

76
vtkm/internal/ConnectivityStructuredInternals.h
View File

@ -92,14 +92,6 @@ public:
return pointIds;
}
VTKM_EXEC_CONT
vtkm::IdComponent GetNumberOfCellsIncidentOnPoint(vtkm::Id pointIndex) const
{
return
(static_cast<vtkm::IdComponent>(pointIndex > 0)
+ static_cast<vtkm::IdComponent>(pointIndex < this->PointDimensions-1));
}
VTKM_EXEC_CONT
vtkm::VecVariable<vtkm::Id,MAX_CELL_TO_POINT>
GetCellsOfPoint(vtkm::Id index) const
@ -226,25 +218,6 @@ public:
return this->GetPointsOfCell(this->FlatToLogicalCellIndex(cellIndex));
}
VTKM_EXEC_CONT
vtkm::IdComponent
GetNumberOfCellsIncidentOnPoint(const SchedulingRangeType &ij) const
{
return
(static_cast<vtkm::IdComponent>((ij[0] > 0) && (ij[1] > 0))
+ static_cast<vtkm::IdComponent>((ij[0] < this->PointDimensions[0]-1) && (ij[1] > 0))
+ static_cast<vtkm::IdComponent>((ij[0] > 0) && (ij[1] < this->PointDimensions[1]-1))
+ static_cast<vtkm::IdComponent>(
(ij[0] < this->PointDimensions[0]-1) && (ij[1] < this->PointDimensions[1]-1)));
}
VTKM_EXEC_CONT
vtkm::IdComponent GetNumberOfCellsIncidentOnPoint(vtkm::Id pointIndex) const
{
return this->GetNumberOfCellsIncidentOnPoint(
this->FlatToLogicalPointIndex(pointIndex));
}
VTKM_EXEC_CONT
vtkm::VecVariable<vtkm::Id,MAX_CELL_TO_POINT>
GetCellsOfPoint(const SchedulingRangeType &ij) const
@ -406,37 +379,6 @@ public:
return this->GetPointsOfCell(this->FlatToLogicalCellIndex(cellIndex));
}
VTKM_EXEC_CONT
vtkm::IdComponent
GetNumberOfCellsIncidentOnPoint(const SchedulingRangeType &ijk) const
{
return (
static_cast<vtkm::IdComponent>((ijk[0] > 0) && (ijk[1] > 0) && (ijk[2] > 0))
+ static_cast<vtkm::IdComponent>((ijk[0] < this->PointDimensions[0]-1) && (ijk[1] > 0) && (ijk[2] > 0))
+ static_cast<vtkm::IdComponent>((ijk[0] > 0) && (ijk[1] < this->PointDimensions[1]-1) && (ijk[2] > 0))
+ static_cast<vtkm::IdComponent>((ijk[0] < this->PointDimensions[0]-1) &&
(ijk[1] < this->PointDimensions[1]-1) &&
(ijk[2] > 0))
+ static_cast<vtkm::IdComponent>((ijk[0] > 0) && (ijk[1] > 0) && (ijk[2] < this->PointDimensions[2]-1))
+ static_cast<vtkm::IdComponent>((ijk[0] < this->PointDimensions[0]-1) &&
(ijk[1] > 0) &&
(ijk[2] < this->PointDimensions[2]-1))
+ static_cast<vtkm::IdComponent>((ijk[0] > 0) &&
(ijk[1] < this->PointDimensions[1]-1) &&
(ijk[2] < this->PointDimensions[2]-1))
+ static_cast<vtkm::IdComponent>((ijk[0] < this->PointDimensions[0]-1) &&
(ijk[1] < this->PointDimensions[1]-1) &&
(ijk[2] < this->PointDimensions[2]-1))
);
}
VTKM_EXEC_CONT
vtkm::IdComponent GetNumberOfCellsIncidentOnPoint(vtkm::Id pointIndex) const
{
return this->GetNumberOfCellsIncidentOnPoint(
this->FlatToLogicalPointIndex(pointIndex));
}
VTKM_EXEC_CONT
vtkm::VecVariable<vtkm::Id,MAX_CELL_TO_POINT>
GetCellsOfPoint(const SchedulingRangeType &ijk) const
@ -578,15 +520,6 @@ struct ConnectivityStructuredIndexHelper<
return connectivity.GetPointsOfCell(cellIndex);
}
template<typename IndexType>
VTKM_EXEC_CONT
static vtkm::IdComponent GetNumberOfIndices(
const ConnectivityType &connectivity,
const IndexType &vtkmNotUsed(cellIndex))
{
return connectivity.GetNumberOfPointsInCell();
}
VTKM_EXEC_CONT
static LogicalIndexType
FlatToLogicalFromIndex(const ConnectivityType &connectivity,
@ -639,15 +572,6 @@ struct ConnectivityStructuredIndexHelper<
return connectivity.GetCellsOfPoint(pointIndex);
}
template<typename IndexType>
VTKM_EXEC_CONT
static vtkm::IdComponent GetNumberOfIndices(
const ConnectivityType &connectivity,
const IndexType &pointIndex)
{
return connectivity.GetNumberOfCellsIncidentOnPoint(pointIndex);
}
VTKM_EXEC_CONT
static LogicalIndexType
FlatToLogicalFromIndex(const ConnectivityType &connectivity,

478
vtkm/worklet/MarchingCubes.h
View File

@ -33,14 +33,18 @@
#include <vtkm/cont/ArrayHandleIndex.h>
#include <vtkm/cont/ArrayHandlePermutation.h>
#include <vtkm/cont/ArrayHandleZip.h>
#include <vtkm/cont/CellSetPermutation.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/cont/DeviceAdapter.h>
#include <vtkm/cont/DynamicArrayHandle.h>
#include <vtkm/cont/Field.h>
#include <vtkm/worklet/DispatcherMapTopology.h>
#include <vtkm/worklet/DispatcherReduceByKey.h>
#include <vtkm/worklet/Keys.h>
#include <vtkm/worklet/ScatterCounting.h>
#include <vtkm/worklet/WorkletMapTopology.h>
#include <vtkm/worklet/WorkletReduceByKey.h>
#include <vtkm/worklet/MarchingCubesDataTables.h>
@ -159,9 +163,7 @@ public:
/// This information is not passed as part of the arguments to the worklet as
/// that dramatically increase compile time by 200%
// -----------------------------------------------------------------------------
template< typename NormalType,
typename NormalStorage,
typename DeviceAdapter >
template< typename DeviceAdapter >
class EdgeWeightGenerateMetaData
{
template<typename FieldType>
@ -174,45 +176,33 @@ class EdgeWeightGenerateMetaData
typedef typename ExecutionTypes::PortalConst PortalConst;
};
struct NormalPortalTypes
{
typedef vtkm::cont::ArrayHandle<vtkm::Vec< NormalType, 3>, NormalStorage> HandleType;
typedef typename HandleType::template ExecutionTypes<DeviceAdapter> ExecutionTypes;
typedef typename ExecutionTypes::Portal Portal;
};
public:
VTKM_CONT
EdgeWeightGenerateMetaData(
vtkm::Id size,
vtkm::cont::ArrayHandle< vtkm::Vec<NormalType, 3>, NormalStorage >& normals,
vtkm::cont::ArrayHandle< vtkm::FloatDefault >& interpWeights,
vtkm::cont::ArrayHandle<vtkm::Id2>& interpIds,
vtkm::cont::ArrayHandle<vtkm::Id>& interpCellIds,
vtkm::cont::ArrayHandle<vtkm::UInt8>& interpContourId,
const vtkm::cont::ArrayHandle< vtkm::IdComponent >& edgeTable,
const vtkm::cont::ArrayHandle< vtkm::IdComponent >& numTriTable,
const vtkm::cont::ArrayHandle< vtkm::IdComponent >& triTable,
const vtkm::worklet::ScatterCounting& scatter):
NormalPortal( normals.PrepareForOutput( 3*size, DeviceAdapter() ) ),
InterpWeightsPortal( interpWeights.PrepareForOutput( 3*size, DeviceAdapter()) ),
InterpIdPortal( interpIds.PrepareForOutput( 3*size, DeviceAdapter() ) ),
InterpCellIdPortal( interpCellIds .PrepareForOutput( 3*size, DeviceAdapter() ) ),
InterpContourPortal( interpContourId.PrepareForOutput( 3*size, DeviceAdapter() ) ),
EdgeTable( edgeTable.PrepareForInput(DeviceAdapter()) ),
NumTriTable( numTriTable.PrepareForInput(DeviceAdapter()) ),
TriTable( triTable.PrepareForInput(DeviceAdapter()) ),
Scatter(scatter)
{
//any way we can easily build an interface so that we don't need to hold
//onto a billion portals?
//Normal and Interp need to be 3 times longer than size as they
//are per point of the output triangle
// Interp needs to be 3 times longer than size as they are per point of the
// output triangle
}
typename NormalPortalTypes::Portal NormalPortal;
typename PortalTypes<vtkm::FloatDefault>::Portal InterpWeightsPortal;
typename PortalTypes<vtkm::Id2>::Portal InterpIdPortal;
typename PortalTypes<vtkm::Id>::Portal InterpCellIdPortal;
typename PortalTypes<vtkm::UInt8>::Portal InterpContourPortal;
typename PortalTypes<vtkm::IdComponent>::PortalConst EdgeTable;
typename PortalTypes<vtkm::IdComponent>::PortalConst NumTriTable;
@ -223,46 +213,39 @@ public:
/// \brief Compute the weights for each edge that is used to generate
/// a point in the resulting iso-surface
// -----------------------------------------------------------------------------
template< typename T,
typename NormalType,
typename NormalStorage,
typename DeviceAdapter >
template< typename T, typename DeviceAdapter >
class EdgeWeightGenerate : public vtkm::worklet::WorkletMapPointToCell
{
public:
struct ClassifyCellTagType : vtkm::ListTagBase< typename float_type<T>::type > { };
struct ClassifyCellTagType : vtkm::ListTagBase<T> { };
typedef vtkm::worklet::ScatterCounting ScatterType;
typedef void ControlSignature(
CellSetIn cellset, // Cell set
WholeArrayIn< ClassifyCellTagType > isoValues,
FieldInPoint< ClassifyCellTagType > fieldIn, // Input point field defining the contour
FieldInPoint<Vec3> pcoordIn // Input point coordinates
FieldInPoint< ClassifyCellTagType > fieldIn // Input point field defining the contour
);
typedef void ExecutionSignature(CellShape, _2, _3, _4, WorkIndex, VisitIndex, FromIndices);
typedef void ExecutionSignature(CellShape, _2, _3, InputIndex, WorkIndex, VisitIndex, FromIndices);
typedef _1 InputDomain;
VTKM_CONT
EdgeWeightGenerate(bool genNormals,
const EdgeWeightGenerateMetaData<NormalType, NormalStorage, DeviceAdapter>& meta) :
GenerateNormals(genNormals),
EdgeWeightGenerate(const EdgeWeightGenerateMetaData<DeviceAdapter>& meta) :
MetaData( meta )
{
}
template<typename IsoValuesType,
typename FieldInType, // Vec-like, one per input point
typename CoordType,
typename IndicesVecType>
VTKM_EXEC
void operator()(
vtkm::CellShapeTagGeneric shape,
const IsoValuesType &isovalues,
const FieldInType & fieldIn, // Input point field defining the contour
const CoordType & coords, // Input point coordinates
vtkm::Id inputCellId,
vtkm::Id outputCellId,
vtkm::IdComponent visitIndex,
const IndicesVecType & indices) const
@ -272,7 +255,7 @@ public:
this->operator()(vtkm::CellShapeTagHexahedron(),
isovalues,
fieldIn,
coords,
inputCellId,
outputCellId,
visitIndex,
indices);
@ -281,14 +264,13 @@ public:
template<typename IsoValuesType,
typename FieldInType, // Vec-like, one per input point
typename CoordType,
typename IndicesVecType>
VTKM_EXEC
void operator()(
CellShapeTagQuad vtkmNotUsed(shape),
const IsoValuesType &vtkmNotUsed(isovalues),
const FieldInType & vtkmNotUsed(fieldIn), // Input point field defining the contour
const CoordType & vtkmNotUsed(coords), // Input point coordinates
vtkm::Id vtkmNotUsed(inputCellId),
vtkm::Id vtkmNotUsed(outputCellId),
vtkm::IdComponent vtkmNotUsed(visitIndex),
const IndicesVecType & vtkmNotUsed(indices) ) const
@ -297,14 +279,13 @@ public:
template<typename IsoValuesType,
typename FieldInType, // Vec-like, one per input point
typename CoordType,
typename IndicesVecType>
VTKM_EXEC
void operator()(
vtkm::CellShapeTagHexahedron shape,
vtkm::CellShapeTagHexahedron,
const IsoValuesType &isovalues,
const FieldInType &fieldIn, // Input point field defining the contour
const CoordType &coords, // Input point coordinates
vtkm::Id inputCellId,
vtkm::Id outputCellId,
vtkm::IdComponent visitIndex,
const IndicesVecType &indices) const
@ -313,20 +294,21 @@ public:
typedef typename vtkm::VecTraits<FieldInType>::ComponentType FieldType;
vtkm::IdComponent sum = 0, caseNumber = 0;
vtkm::Id i=0;
for(i=0; i < isovalues.GetNumberOfValues(); ++i)
vtkm::IdComponent i=0, size = static_cast<vtkm::IdComponent>(isovalues.GetNumberOfValues());
for(i=0; i < size; ++i)
{
const FieldType ivalue = isovalues[i];
// Compute the Marching Cubes case number for this cell. We need to iterate
// the isovalues until the sum >= our visit index. But we need to make
// sure the caseNumber is correct before stoping
caseNumber = ((fieldIn[0] > isovalues[i]) |
(fieldIn[1] > isovalues[i]) << 1 |
(fieldIn[2] > isovalues[i]) << 2 |
(fieldIn[3] > isovalues[i]) << 3 |
(fieldIn[4] > isovalues[i]) << 4 |
(fieldIn[5] > isovalues[i]) << 5 |
(fieldIn[6] > isovalues[i]) << 6 |
(fieldIn[7] > isovalues[i]) << 7);
caseNumber = ((fieldIn[0] > ivalue) |
(fieldIn[1] > ivalue) << 1 |
(fieldIn[2] > ivalue) << 2 |
(fieldIn[3] > ivalue) << 3 |
(fieldIn[4] > ivalue) << 4 |
(fieldIn[5] > ivalue) << 5 |
(fieldIn[6] > ivalue) << 6 |
(fieldIn[7] > ivalue) << 7);
sum += MetaData.NumTriTable.Get(caseNumber);
if(sum > visitIndex)
{
@ -349,7 +331,12 @@ public:
const FieldType fieldValue0 = fieldIn[edgeVertex0];
const FieldType fieldValue1 = fieldIn[edgeVertex1];
// Store the input cell id so that we can properly generate the normals
// in a subsequent call, after we have merged duplicate points
MetaData.InterpCellIdPortal.Set(outputPointId+triVertex, inputCellId);
MetaData.InterpContourPortal.Set(outputPointId+triVertex, static_cast<vtkm::UInt8>(i) );
MetaData.InterpIdPortal.Set(
outputPointId+triVertex,
vtkm::Id2(indices[edgeVertex0],
@ -359,28 +346,7 @@ public:
static_cast<vtkm::FloatDefault>(isovalues[i] - fieldValue0) /
static_cast<vtkm::FloatDefault>(fieldValue1 - fieldValue0);
//need to factor in outputCellId
MetaData.InterpWeightsPortal.Set(outputPointId+triVertex, interpolant);
if(this->GenerateNormals)
{
const vtkm::Vec<vtkm::FloatDefault,3> edgePCoord0 =
vtkm::exec::ParametricCoordinatesPoint(
fieldIn.GetNumberOfComponents(), edgeVertex0, shape, *this);
const vtkm::Vec<vtkm::FloatDefault,3> edgePCoord1 =
vtkm::exec::ParametricCoordinatesPoint(
fieldIn.GetNumberOfComponents(), edgeVertex1, shape, *this);
const vtkm::Vec<vtkm::FloatDefault,3> interpPCoord =
vtkm::Lerp(edgePCoord0, edgePCoord1, interpolant);
//need to factor in outputCellId
MetaData.NormalPortal.Set(outputPointId+triVertex,
vtkm::Normal(vtkm::exec::CellDerivative(
fieldIn, coords, interpPCoord, shape, *this))
);
}
}
}
@ -391,13 +357,11 @@ public:
}
private:
const bool GenerateNormals;
EdgeWeightGenerateMetaData<NormalType, NormalStorage, DeviceAdapter> MetaData;
EdgeWeightGenerateMetaData<DeviceAdapter> MetaData;
void operator=(const EdgeWeightGenerate<T,NormalType,NormalStorage,DeviceAdapter> &) = delete;
void operator=(const EdgeWeightGenerate<T,DeviceAdapter> &) = delete;
};
// ---------------------------------------------------------------------------
class ApplyToField : public vtkm::worklet::WorkletMapField
{
@ -427,17 +391,6 @@ public:
}
};
// ---------------------------------------------------------------------------
struct FirstValueSame
{
template<typename T, typename U>
VTKM_EXEC_CONT bool operator()(const vtkm::Pair<T,U>& a,
const vtkm::Pair<T,U>& b) const
{
return (a.first == b.first);
}
};
// ---------------------------------------------------------------------------
struct MultiContourLess
{
@ -464,45 +417,213 @@ struct MultiContourLess
};
// ---------------------------------------------------------------------------
template<typename KeyType, typename KeyStorage,
typename ValueType, typename ValueStorage,
typename DeviceAdapterTag>
vtkm::cont::ArrayHandle<KeyType, KeyStorage>
MergeDuplicates(const vtkm::cont::ArrayHandle<KeyType, KeyStorage>& input_keys,
vtkm::cont::ArrayHandle<ValueType, ValueStorage> values,
vtkm::cont::ArrayHandle<vtkm::Id>& connectivity,
DeviceAdapterTag)
struct MergeDuplicateValues : vtkm::worklet::WorkletReduceByKey
{
typedef void ControlSignature(KeysIn keys,
ValuesIn<> valuesIn1,
ValuesIn<> valuesIn2,
ReducedValuesOut<> valueOut1,
ReducedValuesOut<> valueOut2);
typedef void ExecutionSignature(_1, _2, _3, _4, _5);
typedef _1 InputDomain;
template<typename T,
typename ValuesInType,
typename Values2InType,
typename ValuesOutType,
typename Values2OutType>
VTKM_EXEC
void operator()(const T &,
const ValuesInType &values1,
const Values2InType &values2,
ValuesOutType &valueOut1,
Values2OutType &valueOut2) const
{
valueOut1 = values1[0];
valueOut2 = values2[0];
}
};
// ---------------------------------------------------------------------------
struct CopyEdgeIds : vtkm::worklet::WorkletMapField
{
typedef void ControlSignature(FieldIn<>, FieldOut<>);
typedef void ExecutionSignature(_1, _2);
typedef _1 InputDomain;
VTKM_EXEC
void operator()(const vtkm::Id2& input, vtkm::Id2& output) const
{
output = input;
}
template<typename T>
VTKM_EXEC
void operator()(const vtkm::Pair<T, vtkm::Id2>& input, vtkm::Id2& output) const
{
output = input.second;
}
};
// ---------------------------------------------------------------------------
template<typename KeyType, typename KeyStorage, typename DeviceAdapterTag>
void
MergeDuplicates(const vtkm::cont::ArrayHandle<KeyType, KeyStorage>& original_keys,
vtkm::cont::ArrayHandle<vtkm::FloatDefault>& weights,
vtkm::cont::ArrayHandle<vtkm::Id2>& edgeIds,
vtkm::cont::ArrayHandle<vtkm::Id>& cellids,
vtkm::cont::ArrayHandle<vtkm::Id>& connectivity,
DeviceAdapterTag)
{
using Algorithm = vtkm::cont::DeviceAdapterAlgorithm<DeviceAdapterTag>;
//1. Copy the input keys as we need both a sorted & unique version and
// the original version to
vtkm::cont::ArrayHandle<KeyType, KeyStorage> keys;
Algorithm::Copy(input_keys, keys);
vtkm::cont::ArrayHandle<KeyType> input_keys;
Algorithm::Copy(original_keys, input_keys);
vtkm::worklet::Keys<KeyType> keys(input_keys, DeviceAdapterTag());
input_keys.ReleaseResources();
//2. Sort by key, making duplicate ids be adjacent so they are eligable
// to be removed by unique
Algorithm::SortByKey(keys, values, marchingcubes::MultiContourLess());
{
vtkm::worklet::DispatcherReduceByKey<MergeDuplicateValues> dispatcher;
vtkm::cont::ArrayHandle<vtkm::Id> writeCells;
vtkm::cont::ArrayHandle<vtkm::FloatDefault> writeWeights;
dispatcher.Invoke(keys,
weights,
cellids,
writeWeights,
writeCells);
weights = writeWeights;
cellids = writeCells;
}
//3. lastly we need to do a unique by key, but since vtkm doesn't
// offer that feature, we use a zip handle.
// We use a custom comparison operator as we only want to compare
// the keys
auto zipped_kv = vtkm::cont::make_ArrayHandleZip(keys, values);
Algorithm::Unique( zipped_kv, marchingcubes::FirstValueSame());
//need to build the new connectivity
auto uniqueKeys = keys.GetUniqueKeys();
Algorithm::LowerBounds(uniqueKeys, original_keys, connectivity, marchingcubes::MultiContourLess());
//4. LowerBounds generates the output cell connections. It does this by
// finding for each interpolationId where it would be inserted in the
// sorted & unique subset, which generates an index value aka the lookup
// value.
//
Algorithm::LowerBounds(keys, input_keys, connectivity, marchingcubes::MultiContourLess());
//5. We need to return the sorted-unique keys as the caller will need
// to hold onto it for interpolation of other fields
return keys;
//update the edge ids
vtkm::worklet::DispatcherMapField<CopyEdgeIds> edgeDispatcher;
edgeDispatcher.Invoke(uniqueKeys, edgeIds);
}
/// \brief Compute the weights for each edge that is used to generate
/// a point in the resulting iso-surface
// -----------------------------------------------------------------------------
template<typename T>
class NormalWorklet : public vtkm::worklet::WorkletMapPointToCell
{
public:
struct ClassifyCellTagType : vtkm::ListTagBase< typename float_type<T>::type > { };
typedef void ControlSignature(
CellSetIn cellset, // Cell set
FieldInPoint< ClassifyCellTagType > field,
FieldInPoint< Vec3 > originalCellPoints,
FieldInCell< Id2Type > edges,
FieldInCell< Scalar > weights,
FieldOutCell< Vec3 > normal
);
typedef void ExecutionSignature(CellShape, FromIndices, _2, _3, _4, _5, _6);
typedef _1 InputDomain;
template <typename CellType, typename IndicesVecType, typename FieldType,
typename CoordType, typename EdgeType, typename WeightType,
typename NormalType>
VTKM_EXEC void operator()(CellType shape, const IndicesVecType& indices,
const FieldType& fieldIn, const CoordType& coords, const EdgeType& edges,
const WeightType& weight, NormalType& normal) const
{
// steps to get the normal calculation to work.
// we need to back compute the correct parametric point
vtkm::IdComponent edgeVertex0 = 0;
vtkm::IdComponent edgeVertex1 = 1;
for (vtkm::IdComponent i = 0; i < indices.GetNumberOfComponents(); ++i)
{
if (indices[i] == edges[0])
{
edgeVertex0 = i;
}
if (indices[i] == edges[1])
{
edgeVertex1 = i;
}
}
auto grad0 = vtkm::exec::CellDerivative(
fieldIn, coords,
vtkm::exec::ParametricCoordinatesPoint(
fieldIn.GetNumberOfComponents(), edgeVertex0, shape, *this),
shape, *this);
auto grad1 = vtkm::exec::CellDerivative(
fieldIn, coords,
vtkm::exec::ParametricCoordinatesPoint(
fieldIn.GetNumberOfComponents(), edgeVertex1, shape, *this),
shape, *this);
normal = vtkm::Normal(vtkm::Lerp(grad0,grad1,weight));
}
};
// ---------------------------------------------------------------------------
//missing we need the original cells field
//missing we need the original cells coordinates
template< typename InputFieldType,
typename InputStorageType,
typename CoordinateSystem,
typename NormalCType,
typename DeviceAdapter >
struct GenerateNormals
{
GenerateNormals(vtkm::cont::ArrayHandle< vtkm::Vec<NormalCType,3> >& normals,
const vtkm::cont::ArrayHandle<InputFieldType, InputStorageType>& field,
const CoordinateSystem& coordinates,
vtkm::cont::ArrayHandle<vtkm::Id>& cellIds,
vtkm::cont::ArrayHandle<vtkm::Id2>& edges,
vtkm::cont::ArrayHandle<vtkm::FloatDefault>& weights):
Field(field),
Coordinates(coordinates),
CellIds(cellIds),
Edges(edges),
Weights(weights),
OutputNormals(normals)
{
}
template <typename CellSetType>
void operator()(const CellSetType &cellset) const
{
vtkm::cont::CellSetPermutation<CellSetType> permutation;
permutation.Fill(this->CellIds, cellset);
using NormalDispatcher = typename vtkm::worklet::DispatcherMapTopology<
NormalWorklet<InputFieldType>,
DeviceAdapter
>;
NormalDispatcher dispatcher;
dispatcher.Invoke( permutation,
marchingcubes::make_ScalarField(this->Field),
this->Coordinates,
this->Edges,
this->Weights,
this->OutputNormals );
}
const vtkm::cont::ArrayHandle<InputFieldType, InputStorageType>& Field;
const CoordinateSystem& Coordinates;
vtkm::cont::ArrayHandle<vtkm::Id>& CellIds;
vtkm::cont::ArrayHandle<vtkm::Id2>& Edges;
vtkm::cont::ArrayHandle<vtkm::FloatDefault>& Weights;
vtkm::cont::ArrayHandle< vtkm::Vec<NormalCType,3> >& OutputNormals;
};
}
/// \brief Compute the isosurface for a uniform grid data set
@ -517,7 +638,7 @@ MarchingCubes(bool mergeDuplicates=true):
NumTrianglesTable(),
TriangleTable(),
InterpolationWeights(),
InterpolationIds()
InterpolationEdgeIds()
{
// Set up the Marching Cubes case tables as part of the filter so that
// we cache these tables in the execution environment between execution runs
@ -601,7 +722,7 @@ void MapFieldOntoIsosurface(const ArrayHandleIn& input,
//todo: need to use the policy to get the correct storage tag for output
applyFieldDispatcher.Invoke(this->InterpolationIds,
applyFieldDispatcher.Invoke(this->InterpolationEdgeIds,
this->InterpolationWeights,
input,
output);
@ -643,8 +764,6 @@ vtkm::cont::CellSetSingleType< >
using GenerateDispatcher = typename vtkm::worklet::DispatcherMapTopology<
EdgeWeightGenerate<ValueType,
NormalType,
StorageTagNormals,
DeviceAdapter
>,
DeviceAdapter
@ -654,28 +773,30 @@ vtkm::cont::CellSetSingleType< >
vtkm::cont::make_ArrayHandle(isovalues, numIsoValues);
// Call the ClassifyCell functor to compute the Marching Cubes case numbers
// for each cell, and the number of vertices to be generated
ClassifyCell<ValueType> classifyCell;
ClassifyDispatcher classifyCellDispatcher(classifyCell);
vtkm::cont::ArrayHandle<vtkm::IdComponent> numOutputTrisPerCell;
{
ClassifyCell<ValueType> classifyCell;
ClassifyDispatcher classifyCellDispatcher(classifyCell);
classifyCellDispatcher.Invoke(isoValuesHandle,
inputField,
cells,
numOutputTrisPerCell,
this->NumTrianglesTable);
}
//Pass 2 Generate the edges
vtkm::cont::ArrayHandle<vtkm::UInt8> contourIds;
vtkm::cont::ArrayHandle<vtkm::Id> originalCellIds;
{
vtkm::worklet::ScatterCounting scatter(numOutputTrisPerCell, DeviceAdapter());
vtkm::cont::ArrayHandle<vtkm::UInt8> contourIds;
EdgeWeightGenerateMetaData< NormalType,
StorageTagNormals,
DeviceAdapter
EdgeWeightGenerateMetaData< DeviceAdapter
> metaData( scatter.GetOutputRange(numOutputTrisPerCell.GetNumberOfValues()),
normals,
this->InterpolationWeights,
this->InterpolationIds,
this->InterpolationEdgeIds,
originalCellIds,
contourIds,
this->EdgeTable,
this->NumTrianglesTable,
@ -683,86 +804,47 @@ vtkm::cont::CellSetSingleType< >
scatter
);
EdgeWeightGenerate<ValueType,
CoordinateType,
StorageTagNormals,
DeviceAdapter
> weightGenerate( withNormals,
metaData);
EdgeWeightGenerate<ValueType, DeviceAdapter> weightGenerate( metaData );
GenerateDispatcher edgeDispatcher(weightGenerate);
edgeDispatcher.Invoke( cells,
//cast to a scalar field if not one, as cellderivative only works on those
marchingcubes::make_ScalarField(isoValuesHandle),
marchingcubes::make_ScalarField(inputField),
coordinateSystem
isoValuesHandle,
inputField
);
}
if(numIsoValues <= 1 || !this->MergeDuplicatePoints)
{ //release memory early that we are not going to need again
contourIds.ReleaseResources();
}
// Now that we have the edge interpolation finished we can generate the
// coordinates, connectivity and resolve duplicate points.
// Given that normals, and point merging are optional it generates the
// following permutations that we need to support
//
//[0] 1 iso-contour
//[1] 1 iso-contour + point merging
//[2] 1 iso-contour + point merging + normals
//[3] 1 iso-contour + normals
//[4] 2+ iso-contour
//[5] 2+ iso-contour + point merging
//[6] 2+ iso-contour + point merging + normals
//[7] 2+ iso-contour + normals
//
// [0], [3], [4], and [7] are easy to implement as they require zero logic
// other than simple connectivity generation. The challenge is the other
// 4 options
vtkm::cont::DataSet output;
vtkm::cont::ArrayHandle< vtkm::Id > connectivity;
if(this->MergeDuplicatePoints)
{
// In all the below cases you will notice that only interpolation ids
// are updated. That is because MergeDuplicates will internally update
// the InterpolationWeights and normals arrays to be the correct for the
// output. But for InterpolationIds we need to do it manually once done
if(withNormals && numIsoValues == 1)
// the InterpolationWeights and InterpolationOriginCellIds arrays to be the correct for the
// output. But for InterpolationEdgeIds we need to do it manually once done
if(numIsoValues == 1)
{
auto&& result = marchingcubes::MergeDuplicates(
this->InterpolationIds, //keys
vtkm::cont::make_ArrayHandleZip(this->InterpolationWeights, normals), //values
connectivity,
DeviceAdapter() );
this->InterpolationIds = result;
}
else if(withNormals && numIsoValues > 1)
{
auto&& result = marchingcubes::MergeDuplicates(
vtkm::cont::make_ArrayHandleZip(contourIds, this->InterpolationIds), //keys
vtkm::cont::make_ArrayHandleZip(this->InterpolationWeights, normals), //values
connectivity,
DeviceAdapter() );
this->InterpolationIds = result.GetStorage().GetSecondArray();
}
else if(!withNormals && numIsoValues == 1)
{
auto&& result = marchingcubes::MergeDuplicates(
this->InterpolationIds, //keys
marchingcubes::MergeDuplicates(
this->InterpolationEdgeIds, //keys
this->InterpolationWeights, //values
connectivity,
this->InterpolationEdgeIds, //values
originalCellIds, //values
connectivity, // computed using lower bounds
DeviceAdapter() );
this->InterpolationIds = result;
}
else if(!withNormals && numIsoValues >= 1)
else if(numIsoValues > 1)
{
auto&& result = marchingcubes::MergeDuplicates(
vtkm::cont::make_ArrayHandleZip(contourIds, this->InterpolationIds), //keys
marchingcubes::MergeDuplicates(
vtkm::cont::make_ArrayHandleZip(contourIds, this->InterpolationEdgeIds), //keys
this->InterpolationWeights, //values
connectivity,
this->InterpolationEdgeIds, //values
originalCellIds, //values
connectivity, // computed using lower bounds
DeviceAdapter() );
this->InterpolationIds = result.GetStorage().GetSecondArray();
}
}
else
@ -771,17 +853,16 @@ vtkm::cont::CellSetSingleType< >
//by a counting array. The danger of doing it this way is that the output
//type is unknown. That is why we copy it into an explicit array
using Algorithm = vtkm::cont::DeviceAdapterAlgorithm<DeviceAdapter>;
vtkm::cont::ArrayHandleIndex temp(this->InterpolationIds.GetNumberOfValues());
vtkm::cont::ArrayHandleIndex temp(this->InterpolationEdgeIds.GetNumberOfValues());
Algorithm::Copy(temp, connectivity);
}
//generate the vertices's
ApplyToField applyToField;
vtkm::worklet::DispatcherMapField<ApplyToField,
DeviceAdapter> applyFieldDispatcher(applyToField);
applyFieldDispatcher.Invoke(this->InterpolationIds,
applyFieldDispatcher.Invoke(this->InterpolationEdgeIds,
this->InterpolationWeights,
coordinateSystem,
vertices);
@ -793,6 +874,23 @@ vtkm::cont::CellSetSingleType< >
3,
connectivity );
//now that the vertices have been generated we can generate the normals
if(withNormals)
{
using GenNormals = marchingcubes::GenerateNormals<ValueType, StorageTagField,
CoordinateSystem, NormalType,
DeviceAdapter>;
GenNormals gen( normals,
inputField,
coordinateSystem,
originalCellIds,
this->InterpolationEdgeIds,
this->InterpolationWeights
);
CastAndCall(cells,gen);
}
return outputCells;
}
@ -803,7 +901,7 @@ vtkm::cont::CellSetSingleType< >
vtkm::cont::ArrayHandle<vtkm::IdComponent> TriangleTable;
vtkm::cont::ArrayHandle<vtkm::FloatDefault> InterpolationWeights;
vtkm::cont::ArrayHandle<vtkm::Id2> InterpolationIds;
vtkm::cont::ArrayHandle<vtkm::Id2> InterpolationEdgeIds;
};
}