mirror of
https://gitlab.kitware.com/vtk/vtk-m
synced 2024-09-19 18:45:43 +00:00
6f03f72b49
The two worklets for marching cubes use tables stored in arrays that have random access. Previously, they arrays were passed using the ExecObject tag in ControlSignature along with ExecutionWholeArrayConst. This changes to using a WholeArrayIn tag and just passing the ArrayHandle directly to the Invoke method. The end result is the same, but the code is a bit cleaner.
321 lines
12 KiB
C++
321 lines
12 KiB
C++
//============================================================================
|
|
// 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 2014 Sandia Corporation.
|
|
// Copyright 2014 UT-Battelle, LLC.
|
|
// Copyright 2014 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_MarchingCubes_h
|
|
#define vtk_m_worklet_MarchingCubes_h
|
|
|
|
#include <vtkm/VectorAnalysis.h>
|
|
|
|
#include <vtkm/exec/CellDerivative.h>
|
|
#include <vtkm/exec/ParametricCoordinates.h>
|
|
|
|
#include <vtkm/cont/ArrayHandle.h>
|
|
#include <vtkm/cont/ArrayHandleCompositeVector.h>
|
|
#include <vtkm/cont/ArrayHandleGroupVec.h>
|
|
#include <vtkm/cont/ArrayHandleIndex.h>
|
|
#include <vtkm/cont/ArrayHandlePermutation.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/ScatterCounting.h>
|
|
#include <vtkm/worklet/WorkletMapTopology.h>
|
|
|
|
#include <vtkm/worklet/MarchingCubesDataTables.h>
|
|
|
|
namespace vtkm {
|
|
namespace worklet {
|
|
|
|
/// \brief Compute the isosurface for a uniform grid data set
|
|
template <typename FieldType, typename DeviceAdapter>
|
|
class MarchingCubes
|
|
{
|
|
public:
|
|
typedef vtkm::cont::ArrayHandle<FieldType> WeightHandle;
|
|
typedef vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Id,2> > IdPairHandle;
|
|
|
|
class ClassifyCell : public vtkm::worklet::WorkletMapPointToCell
|
|
{
|
|
public:
|
|
typedef void ControlSignature(
|
|
FieldInPoint<Scalar> inNodes,
|
|
TopologyIn topology,
|
|
FieldOutCell<> outNumTriangles,
|
|
WholeArrayIn<IdComponentType> numTrianglesTable);
|
|
typedef void ExecutionSignature(_1, _3, _4);
|
|
typedef _2 InputDomain;
|
|
|
|
FieldType Isovalue;
|
|
|
|
VTKM_CONT_EXPORT
|
|
ClassifyCell(FieldType isovalue) :
|
|
Isovalue(isovalue)
|
|
{
|
|
}
|
|
|
|
template<typename InPointVecType,
|
|
typename NumTrianglesTablePortalType>
|
|
VTKM_EXEC_EXPORT
|
|
void operator()(const InPointVecType &fieldIn,
|
|
vtkm::IdComponent &numTriangles,
|
|
const NumTrianglesTablePortalType &numTrianglesTable) const
|
|
{
|
|
vtkm::IdComponent caseNumber =
|
|
( (fieldIn[0] > this->Isovalue)
|
|
| (fieldIn[1] > this->Isovalue)<<1
|
|
| (fieldIn[2] > this->Isovalue)<<2
|
|
| (fieldIn[3] > this->Isovalue)<<3
|
|
| (fieldIn[4] > this->Isovalue)<<4
|
|
| (fieldIn[5] > this->Isovalue)<<5
|
|
| (fieldIn[6] > this->Isovalue)<<6
|
|
| (fieldIn[7] > this->Isovalue)<<7 );
|
|
numTriangles = numTrianglesTable.Get(caseNumber);
|
|
}
|
|
};
|
|
|
|
/// \brief Compute isosurface vertices and scalars
|
|
class IsosurfaceGenerate : public vtkm::worklet::WorkletMapPointToCell
|
|
{
|
|
public:
|
|
typedef typename vtkm::cont::ArrayHandle<vtkm::IdComponent>::
|
|
ExecutionTypes<DeviceAdapter>::PortalConst IdPortalConstType;
|
|
IdPortalConstType EdgeTable;
|
|
|
|
typedef void ControlSignature(
|
|
TopologyIn topology, // Cell set
|
|
FieldInPoint<> fieldIn, // Input point field defining the contour
|
|
FieldInPoint<Vec3> pcoordIn, // Input point coordinates
|
|
FieldOutCell<> interpolationWeights,
|
|
FieldOutCell<> interpolationIds,
|
|
FieldOutCell<> vertexOut, // Vertices for output triangles
|
|
FieldOutCell<> normalsOut, // Estimated normals (one per tri vertex)
|
|
WholeArrayIn<IdComponentType> TriTable // An array portal with the triangle table
|
|
);
|
|
typedef void ExecutionSignature(
|
|
CellShape, _2, _3, _4, _5, _6, _7, _8, VisitIndex, FromIndices);
|
|
|
|
typedef vtkm::worklet::ScatterCounting ScatterType;
|
|
VTKM_CONT_EXPORT
|
|
ScatterType GetScatter() const
|
|
{
|
|
return this->Scatter;
|
|
}
|
|
|
|
template<typename CountArrayType, typename Device>
|
|
VTKM_CONT_EXPORT
|
|
IsosurfaceGenerate(FieldType isovalue,
|
|
const CountArrayType &countArray,
|
|
Device,
|
|
IdPortalConstType edgeTable) : EdgeTable(edgeTable),
|
|
Isovalue(isovalue),
|
|
Scatter(countArray,
|
|
Device()) { }
|
|
|
|
template<typename CellShapeTag,
|
|
typename FieldInType, // Vec-like, one per input point
|
|
typename CoordType, // Vec-like (one per input point) of Vec-3
|
|
typename WeightType,
|
|
typename IdType,
|
|
typename VertexOutType, // Vec-3 of Vec-3 coords (for triangle)
|
|
typename NormalOutType, // Vec-3 of Vec-3
|
|
typename TriTablePortalType, // Array portal
|
|
typename IndicesVecType>
|
|
VTKM_EXEC_EXPORT
|
|
void operator()(
|
|
CellShapeTag shape,
|
|
const FieldInType &fieldIn, // Input point field defining the contour
|
|
const CoordType &coords, // Input point coordinates
|
|
WeightType &interpolationWeights,
|
|
IdType &interpolationIds,
|
|
VertexOutType &vertexOut, // Vertices for output triangles
|
|
NormalOutType &normalsOut, // Estimated normals (one per tri vertex)
|
|
const TriTablePortalType &triTable, // An array portal with the triangle table
|
|
vtkm::IdComponent visitIndex,
|
|
const IndicesVecType &indices) const
|
|
{
|
|
// Compute the Marching Cubes case number for this cell
|
|
vtkm::IdComponent caseNumber =
|
|
( (fieldIn[0] > this->Isovalue)
|
|
| (fieldIn[1] > this->Isovalue)<<1
|
|
| (fieldIn[2] > this->Isovalue)<<2
|
|
| (fieldIn[3] > this->Isovalue)<<3
|
|
| (fieldIn[4] > this->Isovalue)<<4
|
|
| (fieldIn[5] > this->Isovalue)<<5
|
|
| (fieldIn[6] > this->Isovalue)<<6
|
|
| (fieldIn[7] > this->Isovalue)<<7 );
|
|
|
|
// Interpolate for vertex positions and associated scalar values
|
|
const vtkm::Id triTableOffset =
|
|
static_cast<vtkm::Id>(caseNumber*16 + visitIndex*3);
|
|
for (vtkm::IdComponent triVertex = 0; triVertex < 3; triVertex++)
|
|
{
|
|
const vtkm::IdComponent edgeIndex =
|
|
triTable.Get(triTableOffset + triVertex);
|
|
const vtkm::IdComponent edgeVertex0 =
|
|
this->EdgeTable.Get(2*edgeIndex + 0);
|
|
const vtkm::IdComponent edgeVertex1 =
|
|
this->EdgeTable.Get(2*edgeIndex + 1);
|
|
const FieldType fieldValue0 = fieldIn[edgeVertex0];
|
|
const FieldType fieldValue1 = fieldIn[edgeVertex1];
|
|
const FieldType interpolant =
|
|
(this->Isovalue - fieldValue0) / (fieldValue1 - fieldValue0);
|
|
vertexOut[triVertex] = vtkm::Lerp(coords[edgeVertex0],
|
|
coords[edgeVertex1],
|
|
interpolant);
|
|
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);
|
|
normalsOut[triVertex] =
|
|
vtkm::Normal(vtkm::exec::CellDerivative(
|
|
fieldIn, coords, interpPCoord, shape, *this));
|
|
interpolationIds[triVertex][0] = indices[edgeVertex0];
|
|
interpolationIds[triVertex][1] = indices[edgeVertex1];
|
|
interpolationWeights[triVertex] = interpolant;
|
|
}
|
|
}
|
|
|
|
private:
|
|
const FieldType Isovalue;
|
|
ScatterType Scatter;
|
|
};
|
|
|
|
class ApplyToField : public vtkm::worklet::WorkletMapField
|
|
{
|
|
public:
|
|
typedef void ControlSignature(FieldIn<Scalar> interpolationLow,
|
|
FieldIn<Scalar> interpolationHigh,
|
|
FieldIn<Scalar> interpolationWeight,
|
|
FieldOut<Scalar> interpolatedOutput);
|
|
typedef void ExecutionSignature(_1, _2, _3, _4);
|
|
typedef _1 InputDomain;
|
|
|
|
VTKM_CONT_EXPORT
|
|
ApplyToField() {}
|
|
|
|
template <typename Field>
|
|
VTKM_EXEC_EXPORT
|
|
void operator()(const Field& low,
|
|
const Field& high,
|
|
const FieldType &weight,
|
|
Field &result) const
|
|
{
|
|
result = vtkm::Lerp(low, high, weight);
|
|
}
|
|
};
|
|
|
|
|
|
MarchingCubes() {}
|
|
|
|
template<typename CellSetType,typename StorageTag, typename CoordinateType>
|
|
void Run(const float &isovalue,
|
|
const CellSetType& cellSet,
|
|
const vtkm::cont::CoordinateSystem& coordinateSystem,
|
|
const vtkm::cont::ArrayHandle<FieldType, StorageTag>& field,
|
|
vtkm::cont::ArrayHandle< vtkm::Vec<CoordinateType,3> > vertices,
|
|
vtkm::cont::ArrayHandle< vtkm::Vec<CoordinateType,3> > normals)
|
|
{
|
|
// Set up the Marching Cubes case tables
|
|
vtkm::cont::ArrayHandle<vtkm::IdComponent> edgeTable =
|
|
vtkm::cont::make_ArrayHandle(vtkm::worklet::internal::edgeTable,
|
|
24);
|
|
vtkm::cont::ArrayHandle<vtkm::IdComponent> numTrianglesTable =
|
|
vtkm::cont::make_ArrayHandle(vtkm::worklet::internal::numTrianglesTable,
|
|
256);
|
|
vtkm::cont::ArrayHandle<vtkm::IdComponent> triangleTableArray =
|
|
vtkm::cont::make_ArrayHandle(vtkm::worklet::internal::triTable,
|
|
256*16);
|
|
|
|
vtkm::cont::ArrayHandle<vtkm::IdComponent> numOutputTrisPerCell;
|
|
|
|
// Call the ClassifyCell functor to compute the Marching Cubes case numbers
|
|
// for each cell, and the number of vertices to be generated
|
|
ClassifyCell classifyCell(isovalue);
|
|
|
|
typedef typename vtkm::worklet::DispatcherMapTopology<
|
|
ClassifyCell,
|
|
DeviceAdapter> ClassifyCellDispatcher;
|
|
ClassifyCellDispatcher classifyCellDispatcher(classifyCell);
|
|
|
|
classifyCellDispatcher.Invoke(field,
|
|
cellSet,
|
|
numOutputTrisPerCell,
|
|
numTrianglesTable);
|
|
|
|
IsosurfaceGenerate isosurface(isovalue,
|
|
numOutputTrisPerCell,
|
|
DeviceAdapter(),
|
|
edgeTable.PrepareForInput(DeviceAdapter()));
|
|
|
|
vtkm::worklet::DispatcherMapTopology<IsosurfaceGenerate, DeviceAdapter>
|
|
isosurfaceDispatcher(isosurface);
|
|
isosurfaceDispatcher.Invoke(
|
|
cellSet,
|
|
field,
|
|
coordinateSystem.GetData(),
|
|
vtkm::cont::make_ArrayHandleGroupVec<3>(this->InterpolationWeights),
|
|
vtkm::cont::make_ArrayHandleGroupVec<3>(this->InterpolationIds),
|
|
vtkm::cont::make_ArrayHandleGroupVec<3>(vertices),
|
|
vtkm::cont::make_ArrayHandleGroupVec<3>(normals),
|
|
triangleTableArray);
|
|
}
|
|
|
|
template<typename ArrayHandleIn, typename ArrayHandleOut>
|
|
void MapFieldOntoIsosurface(const ArrayHandleIn& fieldIn,
|
|
ArrayHandleOut& fieldOut)
|
|
{
|
|
typedef typename vtkm::cont::ArrayHandleCompositeVectorType<IdPairHandle>
|
|
::type IdType;
|
|
typedef vtkm::cont::ArrayHandlePermutation<IdType,ArrayHandleIn>
|
|
FieldPermutationHandleType;
|
|
|
|
FieldPermutationHandleType
|
|
low(vtkm::cont::make_ArrayHandleCompositeVector(this->InterpolationIds,
|
|
0),fieldIn);
|
|
FieldPermutationHandleType
|
|
high(vtkm::cont::make_ArrayHandleCompositeVector(this->InterpolationIds,
|
|
1),fieldIn);
|
|
|
|
ApplyToField applyToField;
|
|
vtkm::worklet::DispatcherMapField<ApplyToField,DeviceAdapter>
|
|
applyToFieldDispatcher(applyToField);
|
|
|
|
applyToFieldDispatcher.Invoke(low,
|
|
high,
|
|
this->InterpolationWeights,
|
|
fieldOut);
|
|
}
|
|
|
|
private:
|
|
WeightHandle InterpolationWeights;
|
|
IdPairHandle InterpolationIds;
|
|
};
|
|
|
|
}
|
|
} // namespace vtkm::worklet
|
|
|
|
#endif // vtk_m_worklet_MarchingCubes_h
|