Shuenhoy/vtk-m
1
0
Fork 0
mirror of https://gitlab.kitware.com/vtk/vtk-m synced 2024-09-08 13:23:51 +00:00
Code Issues 2 Actions Packages Projects Releases Wiki Activity

Allow MarchingCubes to handle multiple iso-contour values.

Browse Source
This commit is contained in:
Robert Maynard 2017-03-14 08:56:44 -04:00
parent ed1f568d2d
commit 5566afdd8d
8 changed files with 262 additions and 126 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
examples/demo/Demo.cxx
View File

@ -124,7 +124,7 @@ int main(int argc, char* argv[])
vtkm::filter::MarchingCubes filter;
filter.SetGenerateNormals(false);
filter.SetMergeDuplicatePoints(false);
filter.SetIsoValue(isovalue);
filter.SetIsoValue(0, isovalue);
vtkm::filter::ResultDataSet result = filter.Execute( inputData,
inputData.GetField(fieldName) );
filter.MapFieldOntoOutput(result, inputData.GetField(fieldName));

4
examples/isosurface/IsosurfaceUniformGrid.cxx
View File

@ -234,8 +234,8 @@ int main(int argc, char* argv[])
vtkm::filter::MarchingCubes filter;
filter.SetGenerateNormals(true);
filter.SetMergeDuplicatePoints( false );
filter.SetIsoValue( 0.5 );
filter.SetMergeDuplicatePoints(false);
filter.SetIsoValue(0, 0.5);
vtkm::filter::ResultDataSet result =
filter.Execute( dataSet, dataSet.GetField("nodevar") );

6
vtkm/exec/ExecutionWholeArray.h
View File

@ -68,6 +68,9 @@ public:
VTKM_EXEC
T Get(vtkm::Id index) const { return this->Portal.Get(index); }
VTKM_EXEC
T operator[](vtkm::Id index) const { return this->Portal.Get(index); }
VTKM_EXEC
void Set(vtkm::Id index, const T& t) const { this->Portal.Set(index, t); }
@ -112,6 +115,9 @@ public:
VTKM_EXEC
T Get(vtkm::Id index) const { return this->Portal.Get(index); }
VTKM_EXEC
T operator[](vtkm::Id index) const { return this->Portal.Get(index); }
VTKM_EXEC
const PortalType& GetPortal() const { return this->Portal; }

18
vtkm/filter/MarchingCubes.h
View File

@ -45,10 +45,22 @@ public:
MarchingCubes();
VTKM_CONT
void SetIsoValue(vtkm::Float64 value){ this->IsoValue = value; }
void SetNumberOfIsoValues(vtkm::Id num);
VTKM_CONT
vtkm::Float64 GetIsoValue() const { return this->IsoValue; }
vtkm::Id GetNumberOfIsoValues() const;
VTKM_CONT
void SetIsoValue(vtkm::Float64 v) { this->SetIsoValue(0, v); }
VTKM_CONT
void SetIsoValue(vtkm::Id index, vtkm::Float64);
VTKM_CONT
void SetIsoValues(const std::vector<vtkm::Float64>& values);
VTKM_CONT
vtkm::Float64 GetIsoValue(vtkm::Id index) const;
VTKM_CONT
void SetMergeDuplicatePoints(bool on) { this->Worklet.SetMergeDuplicatePoints(on); }
@ -87,7 +99,7 @@ public:
const DeviceAdapter& tag);
private:
double IsoValue;
std::vector<vtkm::Float64> IsoValues;
bool GenerateNormals;
std::string NormalArrayName;
vtkm::worklet::MarchingCubes Worklet;

63
vtkm/filter/MarchingCubes.hxx
View File

@ -32,7 +32,7 @@ namespace filter {
inline VTKM_CONT
MarchingCubes::MarchingCubes():
vtkm::filter::FilterDataSetWithField<MarchingCubes>(),
IsoValue(0),
IsoValues(),
GenerateNormals(false),
NormalArrayName("normals"),
Worklet()
@ -40,6 +40,50 @@ MarchingCubes::MarchingCubes():
// todo: keep an instance of marching cubes worklet as a member variable
}
//-----------------------------------------------------------------------------
inline
void MarchingCubes::SetNumberOfIsoValues(vtkm::Id num)
{
if(num >= 0)
{
this->IsoValues.resize(static_cast<std::size_t>(num));
}
}
//-----------------------------------------------------------------------------
inline
vtkm::Id MarchingCubes::GetNumberOfIsoValues() const
{
return static_cast<vtkm::Id>(this->IsoValues.size());
}
//-----------------------------------------------------------------------------
inline
void MarchingCubes::SetIsoValue(vtkm::Id index, vtkm::Float64 v)
{
std::size_t i = static_cast<std::size_t>(index);
if(i >= this->IsoValues.size())
{
this->IsoValues.resize(i+1);
}
this->IsoValues[i] = v;
}
//-----------------------------------------------------------------------------
inline
void MarchingCubes::SetIsoValues(const std::vector<vtkm::Float64>& values)
{
this->IsoValues = values;
}
//-----------------------------------------------------------------------------
inline
vtkm::Float64 MarchingCubes::GetIsoValue(vtkm::Id index) const
{
return this->IsoValues[static_cast<std::size_t>(index)];
}
//-----------------------------------------------------------------------------
template<typename T,
typename StorageType,
@ -59,6 +103,11 @@ vtkm::filter::ResultDataSet MarchingCubes::DoExecute(const vtkm::cont::DataSet&
return vtkm::filter::ResultDataSet();
}
if(this->IsoValues.size() == 0)
{
return vtkm::filter::ResultDataSet();
}
//get the cells and coordinates of the dataset
const vtkm::cont::DynamicCellSet& cells =
input.GetCellSet(this->GetActiveCellSetIndex());
@ -73,6 +122,12 @@ vtkm::filter::ResultDataSet MarchingCubes::DoExecute(const vtkm::cont::DataSet&
vtkm::cont::DataSet output;
vtkm::cont::CellSetSingleType< > outputCells;
std::vector<T> ivalues(this->IsoValues.size());
for(std::size_t i = 0; i < ivalues.size(); ++i)
{
ivalues[i] = static_cast<T>(this->IsoValues[i]);
}
//not sold on this as we have to generate more signatures for the
//worklet with the design
//But I think we should get this to compile before we tinker with
@ -80,7 +135,8 @@ vtkm::filter::ResultDataSet MarchingCubes::DoExecute(const vtkm::cont::DataSet&
if(this->GenerateNormals)
{
outputCells =
this->Worklet.Run( static_cast<T>(this->IsoValue),
this->Worklet.Run( &ivalues[0],
static_cast<vtkm::Id>(ivalues.size()),
vtkm::filter::ApplyPolicy(cells, policy),
vtkm::filter::ApplyPolicy(coords, policy),
field,
@ -92,7 +148,8 @@ vtkm::filter::ResultDataSet MarchingCubes::DoExecute(const vtkm::cont::DataSet&
else
{
outputCells =
this->Worklet.Run( static_cast<T>(this->IsoValue),
this->Worklet.Run( &ivalues[0],
static_cast<vtkm::Id>(ivalues.size()),
vtkm::filter::ApplyPolicy(cells, policy),
vtkm::filter::ApplyPolicy(coords, policy),
field,

7
vtkm/filter/testing/UnitTestMarchingCubesFilter.cxx
View File

@ -267,7 +267,7 @@ void TestMarchingCubesUniformGrid()
vtkm::filter::MarchingCubes mc;
mc.SetGenerateNormals(true);
mc.SetIsoValue( 0.5 );
mc.SetIsoValue(0, 0.5);
result = mc.Execute( dataSet,
dataSet.GetField("nodevar") );
@ -339,7 +339,10 @@ void TestMarchingCubesCustomPolicy()
vtkm::filter::MarchingCubes mc;
mc.SetGenerateNormals( false );
mc.SetIsoValue( 0.45 );
mc.SetIsoValue(0, 0.45);
mc.SetIsoValue(1, 0.45);
mc.SetIsoValue(2, 0.45);
mc.SetIsoValue(3, 0.45);
//We specify a custom execution policy here, since the contourField is a
//custom field type

281
vtkm/worklet/MarchingCubes.h
View File

@ -48,6 +48,11 @@ namespace worklet {
namespace marchingcubes {
template<typename T> struct float_type { using type = vtkm::FloatDefault; };
template<> struct float_type< vtkm::Float32 > { using type = vtkm::Float32; };
template<> struct float_type< vtkm::Float64 > { using type = vtkm::Float64; };
// -----------------------------------------------------------------------------
template<typename S>
vtkm::cont::ArrayHandle<vtkm::Float32,S> make_ScalarField(const vtkm::cont::ArrayHandle<vtkm::Float32,S>& ah)
@ -67,21 +72,6 @@ vtkm::cont::ArrayHandleCast<vtkm::FloatDefault, vtkm::cont::ArrayHandle<vtkm::In
make_ScalarField(const vtkm::cont::ArrayHandle<vtkm::Int8,S>& ah)
{ return vtkm::cont::make_ArrayHandleCast(ah, vtkm::FloatDefault()); }
// -----------------------------------------------------------------------------
template<typename T, typename U>
VTKM_EXEC
int GetHexahedronClassification(const T& values, const U isoValue)
{
return ((values[0] > isoValue) |
(values[1] > isoValue) << 1 |
(values[2] > isoValue) << 2 |
(values[3] > isoValue) << 3 |
(values[4] > isoValue) << 4 |
(values[5] > isoValue) << 5 |
(values[6] > isoValue) << 6 |
(values[7] > isoValue) << 7);
}
// ---------------------------------------------------------------------------
template<typename T>
class ClassifyCell : public vtkm::worklet::WorkletMapPointToCell
@ -90,46 +80,45 @@ public:
struct ClassifyCellTagType : vtkm::ListTagBase<T> { };
typedef void ControlSignature(
FieldInPoint< ClassifyCellTagType > inNodes,
WholeArrayIn< ClassifyCellTagType > isoValues,
FieldInPoint< ClassifyCellTagType > fieldIn,
CellSetIn cellset,
FieldOutCell< IdComponentType > outNumTriangles,
WholeArrayIn< IdComponentType > numTrianglesTable);
typedef void ExecutionSignature(CellShape,_1, _3, _4);
typedef _2 InputDomain;
typedef void ExecutionSignature(CellShape,_1, _2, _4, _5);
typedef _3 InputDomain;
T Isovalue;
VTKM_CONT
ClassifyCell(T isovalue) :
Isovalue(isovalue)
{
}
template<typename FieldInType,
template<typename IsoValuesType,
typename FieldInType,
typename NumTrianglesTablePortalType>
VTKM_EXEC
void operator()(vtkm::CellShapeTagGeneric shape,
const IsoValuesType& isovalues,
const FieldInType &fieldIn,
vtkm::IdComponent &numTriangles,
const NumTrianglesTablePortalType &numTrianglesTable) const
{
if(shape.Id == CELL_SHAPE_HEXAHEDRON )
{
{
this->operator()(vtkm::CellShapeTagHexahedron(),
isovalues,
fieldIn,
numTriangles,
numTrianglesTable);
}
}
else
{
numTriangles = 0;
}
}
template<typename FieldInType,
template<typename IsoValuesType,
typename FieldInType,
typename NumTrianglesTablePortalType>
VTKM_EXEC
void operator()(vtkm::CellShapeTagQuad vtkmNotUsed(shape),
const IsoValuesType &vtkmNotUsed(isovalues),
const FieldInType &vtkmNotUsed(fieldIn),
vtkm::IdComponent &vtkmNotUsed(numTriangles),
const NumTrianglesTablePortalType
@ -137,26 +126,30 @@ public:
{
}
template<typename FieldInType,
template<typename IsoValuesType,
typename FieldInType,
typename NumTrianglesTablePortalType>
VTKM_EXEC
void operator()(vtkm::CellShapeTagHexahedron vtkmNotUsed(shape),
const IsoValuesType& isovalues,
const FieldInType &fieldIn,
vtkm::IdComponent &numTriangles,
const NumTrianglesTablePortalType &numTrianglesTable) const
{
typedef typename vtkm::VecTraits<FieldInType>::ComponentType FieldType;
const FieldType iso = static_cast<FieldType>(this->Isovalue);
const vtkm::IdComponent caseNumber = ((fieldIn[0] > iso) |
(fieldIn[1] > iso) << 1 |
(fieldIn[2] > iso) << 2 |
(fieldIn[3] > iso) << 3 |
(fieldIn[4] > iso) << 4 |
(fieldIn[5] > iso) << 5 |
(fieldIn[6] > iso) << 6 |
(fieldIn[7] > iso) << 7);
numTriangles = numTrianglesTable.Get(caseNumber);
vtkm::IdComponent sum = 0;
for(vtkm::Id i=0; i < isovalues.GetNumberOfValues(); ++i)
{
const vtkm::IdComponent 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);
sum += numTrianglesTable.Get(caseNumber);
}
numTriangles = sum;
}
};
@ -165,8 +158,8 @@ public:
/// This information is not passed as part of the arguments to the worklet as
/// that dramatically increase compile time by 200%
// -----------------------------------------------------------------------------
template< typename ScalarType,
typename NormalStorageType,
template< typename NormalType,
typename NormalStorage,
typename DeviceAdapter >
class EdgeWeightGenerateMetaData
{
@ -182,7 +175,7 @@ class EdgeWeightGenerateMetaData
struct NormalPortalTypes
{
typedef vtkm::cont::ArrayHandle<vtkm::Vec< ScalarType, 3>, NormalStorageType> HandleType;
typedef vtkm::cont::ArrayHandle<vtkm::Vec< NormalType, 3>, NormalStorage> HandleType;
typedef typename HandleType::template ExecutionTypes<DeviceAdapter> ExecutionTypes;
typedef typename ExecutionTypes::Portal Portal;
@ -192,9 +185,10 @@ public:
VTKM_CONT
EdgeWeightGenerateMetaData(
vtkm::Id size,
vtkm::cont::ArrayHandle< vtkm::Vec<ScalarType, 3>, NormalStorageType >& normals,
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::UInt8>& interpContourId,
const vtkm::cont::ArrayHandle< vtkm::IdComponent >& edgeTable,
const vtkm::cont::ArrayHandle< vtkm::IdComponent >& numTriTable,
const vtkm::cont::ArrayHandle< vtkm::IdComponent >& triTable,
@ -202,6 +196,7 @@ public:
NormalPortal( normals.PrepareForOutput( 3*size, DeviceAdapter() ) ),
InterpWeightsPortal( interpWeights.PrepareForOutput( 3*size, DeviceAdapter()) ),
InterpIdPortal( interpIds.PrepareForOutput( 3*size, DeviceAdapter() ) ),
InterpContourPortal( interpContourId.PrepareForOutput( 3*size, DeviceAdapter() ) ),
EdgeTable( edgeTable.PrepareForInput(DeviceAdapter()) ),
NumTriTable( numTriTable.PrepareForInput(DeviceAdapter()) ),
TriTable( triTable.PrepareForInput(DeviceAdapter()) ),
@ -217,6 +212,7 @@ public:
typename NormalPortalTypes::Portal NormalPortal;
typename PortalTypes<vtkm::FloatDefault>::Portal InterpWeightsPortal;
typename PortalTypes<vtkm::Id2>::Portal InterpIdPortal;
typename PortalTypes<vtkm::UInt8>::Portal InterpContourPortal;
typename PortalTypes<vtkm::IdComponent>::PortalConst EdgeTable;
typename PortalTypes<vtkm::IdComponent>::PortalConst NumTriTable;
typename PortalTypes<vtkm::IdComponent>::PortalConst TriTable;
@ -226,61 +222,70 @@ public:
/// \brief Compute the weights for each edge that is used to generate
/// a point in the resulting iso-surface
// -----------------------------------------------------------------------------
template< typename ScalarType,
typename NormalStorageType,
template< typename T,
typename NormalType,
typename NormalStorage,
typename DeviceAdapter >
class EdgeWeightGenerate : public vtkm::worklet::WorkletMapPointToCell
{
public:
struct ClassifyCellTagType : vtkm::ListTagBase< typename float_type<T>::type > { };
typedef vtkm::worklet::ScatterCounting ScatterType;
typedef void ControlSignature(
CellSetIn cellset, // Cell set
FieldInPoint<Scalar> fieldIn, // Input point field defining the contour
WholeArrayIn< ClassifyCellTagType > isoValues,
FieldInPoint< ClassifyCellTagType > fieldIn, // Input point field defining the contour
FieldInPoint<Vec3> pcoordIn // Input point coordinates
);
typedef void ExecutionSignature(CellShape, _2, _3, WorkIndex, VisitIndex, FromIndices);
typedef void ExecutionSignature(CellShape, _2, _3, _4, WorkIndex, VisitIndex, FromIndices);
typedef _1 InputDomain;
VTKM_CONT
EdgeWeightGenerate(vtkm::Float64 isovalue,
bool genNormals,
const EdgeWeightGenerateMetaData<ScalarType, NormalStorageType, DeviceAdapter>& meta) :
Isovalue(isovalue),
EdgeWeightGenerate(bool genNormals,
const EdgeWeightGenerateMetaData<NormalType, NormalStorage, DeviceAdapter>& meta) :
GenerateNormals(genNormals),
MetaData( meta )
{
}
template<typename FieldInType, // Vec-like, one per input point
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 outputCellId,
vtkm::IdComponent visitIndex,
const IndicesVecType & indices) const
{ //covers when we have hexs coming from unstructured data
VTKM_ASSUME( shape.Id == CELL_SHAPE_HEXAHEDRON );
this->operator()(vtkm::CellShapeTagHexahedron(),
fieldIn,
coords,
outputCellId,
visitIndex,
indices);
if(shape.Id == CELL_SHAPE_HEXAHEDRON )
{
this->operator()(vtkm::CellShapeTagHexahedron(),
isovalues,
fieldIn,
coords,
outputCellId,
visitIndex,
indices);
}
}
template<typename FieldInType, // Vec-like, one per input point
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(outputCellId),
@ -289,12 +294,14 @@ public:
{ //covers when we have quads coming from 2d structured data
}
template<typename FieldInType, // Vec-like, one per input point
template<typename IsoValuesType,
typename FieldInType, // Vec-like, one per input point
typename CoordType,
typename IndicesVecType>
VTKM_EXEC
void operator()(
vtkm::CellShapeTagHexahedron shape,
const IsoValuesType &isovalues,
const FieldInType &fieldIn, // Input point field defining the contour
const CoordType &coords, // Input point coordinates
vtkm::Id outputCellId,
@ -303,17 +310,31 @@ public:
{ //covers when we have hexs coming from 3d structured data
const vtkm::Id outputPointId = 3 * outputCellId;
typedef typename vtkm::VecTraits<FieldInType>::ComponentType FieldType;
const FieldType iso = static_cast<FieldType>(this->Isovalue);
// Compute the Marching Cubes case number for this cell
const vtkm::IdComponent caseNumber = ((fieldIn[0] > iso) |
(fieldIn[1] > iso) << 1 |
(fieldIn[2] > iso) << 2 |
(fieldIn[3] > iso) << 3 |
(fieldIn[4] > iso) << 4 |
(fieldIn[5] > iso) << 5 |
(fieldIn[6] > iso) << 6 |
(fieldIn[7] > iso) << 7);
vtkm::IdComponent sum = 0, caseNumber = 0;
vtkm::Id i=0;
for(i=0; i < isovalues.GetNumberOfValues(); ++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);
sum += MetaData.NumTriTable.Get(caseNumber);
if(sum > visitIndex)
{
break;
}
}
visitIndex = sum - visitIndex - 1;
// Interpolate for vertex positions and associated scalar values
const vtkm::Id triTableOffset =
@ -327,13 +348,14 @@ public:
const FieldType fieldValue0 = fieldIn[edgeVertex0];
const FieldType fieldValue1 = fieldIn[edgeVertex1];
//need to factor in outputCellId
MetaData.InterpContourPortal.Set(outputPointId+triVertex, static_cast<vtkm::UInt8>(i) );
MetaData.InterpIdPortal.Set(
outputPointId+triVertex,
vtkm::Id2(indices[edgeVertex0], indices[edgeVertex1]));
vtkm::Id2(indices[edgeVertex0],
indices[edgeVertex1]));
vtkm::FloatDefault interpolant =
static_cast<vtkm::FloatDefault>(iso - fieldValue0) /
static_cast<vtkm::FloatDefault>(isovalues[i] - fieldValue0) /
static_cast<vtkm::FloatDefault>(fieldValue1 - fieldValue0);
//need to factor in outputCellId
@ -368,11 +390,10 @@ public:
}
private:
const vtkm::Float64 Isovalue;
const bool GenerateNormals;
EdgeWeightGenerateMetaData<ScalarType, NormalStorageType, DeviceAdapter> MetaData;
EdgeWeightGenerateMetaData<NormalType, NormalStorage, DeviceAdapter> MetaData;
void operator=(const EdgeWeightGenerate<ScalarType,NormalStorageType,DeviceAdapter> &) = delete;
void operator=(const EdgeWeightGenerate<T,NormalType,NormalStorage,DeviceAdapter> &) = delete;
};
@ -465,7 +486,8 @@ template<typename ValueType,
typename StorageTagVertices,
typename DeviceAdapter>
vtkm::cont::CellSetSingleType< >
Run(const ValueType &isovalue,
Run(const ValueType* const isovalues,
const vtkm::Id numIsoValues,
const CellSetType& cells,
const CoordinateSystem& coordinateSystem,
const vtkm::cont::ArrayHandle<ValueType, StorageTagField>& input,
@ -473,7 +495,7 @@ vtkm::cont::CellSetSingleType< >
const DeviceAdapter& device)
{
vtkm::cont::ArrayHandle< vtkm::Vec<CoordinateType,3> > normals;
return this->DoRun(isovalue,cells,coordinateSystem,input,vertices,normals,false,device);
return this->DoRun(isovalues,numIsoValues,cells,coordinateSystem,input,vertices,normals,false,device);
}
//----------------------------------------------------------------------------
@ -486,7 +508,8 @@ template<typename ValueType,
typename StorageTagNormals,
typename DeviceAdapter>
vtkm::cont::CellSetSingleType< >
Run(const ValueType &isovalue,
Run(const ValueType* const isovalues,
const vtkm::Id numIsoValues,
const CellSetType& cells,
const CoordinateSystem& coordinateSystem,
const vtkm::cont::ArrayHandle<ValueType, StorageTagField>& input,
@ -494,7 +517,7 @@ vtkm::cont::CellSetSingleType< >
vtkm::cont::ArrayHandle< vtkm::Vec<CoordinateType,3>, StorageTagNormals > normals,
const DeviceAdapter& device)
{
return this->DoRun(isovalue,cells,coordinateSystem,input,vertices,normals,true,device);
return this->DoRun(isovalues,numIsoValues,cells,coordinateSystem,input,vertices,normals,true,device);
}
//----------------------------------------------------------------------------
@ -528,16 +551,18 @@ template<typename ValueType,
typename StorageTagVertices,
typename StorageTagNormals,
typename CoordinateType,
typename NormalType,
typename DeviceAdapter>
vtkm::cont::CellSetSingleType< >
DoRun(const ValueType &isovalue,
const CellSetType& cells,
const CoordinateSystem& coordinateSystem,
const vtkm::cont::ArrayHandle<ValueType, StorageTagField>& inputField,
vtkm::cont::ArrayHandle< vtkm::Vec<CoordinateType,3>, StorageTagVertices > vertices,
vtkm::cont::ArrayHandle< vtkm::Vec<CoordinateType,3>, StorageTagNormals > normals,
bool withNormals,
const DeviceAdapter& )
DoRun(const ValueType* isovalues,
const vtkm::Id numIsoValues,
const CellSetType& cells,
const CoordinateSystem& coordinateSystem,
const vtkm::cont::ArrayHandle<ValueType, StorageTagField>& inputField,
vtkm::cont::ArrayHandle< vtkm::Vec<CoordinateType,3>, StorageTagVertices > vertices,
vtkm::cont::ArrayHandle< vtkm::Vec<NormalType,3>, StorageTagNormals > normals,
bool withNormals,
const DeviceAdapter& )
{
using vtkm::worklet::marchingcubes::ApplyToField;
using vtkm::worklet::marchingcubes::EdgeWeightGenerate;
@ -551,21 +576,24 @@ vtkm::cont::CellSetSingleType< >
> ClassifyDispatcher;
typedef typename vtkm::worklet::DispatcherMapTopology<
EdgeWeightGenerate<CoordinateType,
EdgeWeightGenerate<ValueType,
NormalType,
StorageTagNormals,
DeviceAdapter
>,
DeviceAdapter
> GenerateDispatcher;
vtkm::cont::ArrayHandle<ValueType> isoValuesHandle =
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( isovalue );
ClassifyCell<ValueType> classifyCell;
ClassifyDispatcher classifyCellDispatcher(classifyCell);
vtkm::cont::ArrayHandle<vtkm::IdComponent> numOutputTrisPerCell;
classifyCellDispatcher.Invoke(inputField,
classifyCellDispatcher.Invoke(isoValuesHandle,
inputField,
cells,
numOutputTrisPerCell,
this->NumTrianglesTable);
@ -574,32 +602,40 @@ vtkm::cont::CellSetSingleType< >
//Pass 2 Generate the edges
vtkm::worklet::ScatterCounting scatter(numOutputTrisPerCell, DeviceAdapter());
EdgeWeightGenerateMetaData< CoordinateType,
vtkm::cont::ArrayHandle<vtkm::UInt8> contourIds;
EdgeWeightGenerateMetaData< NormalType,
StorageTagNormals,
DeviceAdapter
> metaData( scatter.GetOutputRange(numOutputTrisPerCell.GetNumberOfValues()),
normals,
this->InterpolationWeights,
this->InterpolationIds,
contourIds,
this->EdgeTable,
this->NumTrianglesTable,
this->TriangleTable,
scatter
);
EdgeWeightGenerate<CoordinateType,
EdgeWeightGenerate<ValueType,
CoordinateType,
StorageTagNormals,
DeviceAdapter
> weightGenerate( isovalue,
withNormals,
> weightGenerate( withNormals,
metaData);
GenerateDispatcher edgeDispatcher(weightGenerate);
edgeDispatcher.Invoke( cells,
//cast to a scalar field if not one, as cellderivative only works on those
marchingcubes::make_ScalarField(inputField),
coordinateSystem
);
//cast to a scalar field if not one, as cellderivative only works on those
marchingcubes::make_ScalarField(isoValuesHandle),
marchingcubes::make_ScalarField(inputField),
coordinateSystem
);
if(numIsoValues <= 1)
{ //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
//following:
@ -612,6 +648,7 @@ vtkm::cont::CellSetSingleType< >
vtkm::cont::ArrayHandle< vtkm::Id > connectivity;
typedef vtkm::cont::ArrayHandle< vtkm::Id2 > Id2HandleType;
typedef vtkm::cont::ArrayHandle< vtkm::UInt8 > ContourIdHandleType;
typedef vtkm::cont::ArrayHandle<vtkm::FloatDefault> WeightHandleType;
if(this->MergeDuplicatePoints)
{
@ -621,13 +658,23 @@ vtkm::cont::CellSetSingleType< >
Algorithm::Copy(this->InterpolationIds, uniqueIds);
if(withNormals)
{
{
typedef vtkm::cont::ArrayHandle< vtkm::Vec<CoordinateType,3>, StorageTagNormals > NormalHandlType;
typedef vtkm::cont::ArrayHandleZip<WeightHandleType, NormalHandlType> KeyType;
KeyType keys = vtkm::cont::make_ArrayHandleZip(this->InterpolationWeights, normals);
//2. now we need to do a sort by key, making duplicate ids be adjacent
Algorithm::SortByKey(uniqueIds, keys);
if(numIsoValues > 1)
{
vtkm::cont::ArrayHandleZip<
Id2HandleType, ContourIdHandleType> uniqueIdsWithContourId =
vtkm::cont::make_ArrayHandleZip(uniqueIds, contourIds);
Algorithm::SortByKey(uniqueIdsWithContourId, keys);
}
else
{
Algorithm::SortByKey(uniqueIds, keys);
}
//3. lastly we need to do a unique by key, but since vtkm doesn't
// offer that feature, we use a zip handle.
@ -636,11 +683,22 @@ vtkm::cont::CellSetSingleType< >
vtkm::cont::ArrayHandleZip<Id2HandleType, KeyType> zipped =
vtkm::cont::make_ArrayHandleZip(uniqueIds,keys);
Algorithm::Unique( zipped, marchingcubes::FirstValueSame());
}
}
else
{
{
//2. now we need to do a sort by key, making duplicate ids be adjacent
Algorithm::SortByKey(uniqueIds, this->InterpolationWeights);
if(numIsoValues > 1)
{
vtkm::cont::ArrayHandleZip<
Id2HandleType, ContourIdHandleType> uniqueIdsWithContourId =
vtkm::cont::make_ArrayHandleZip(uniqueIds, contourIds);
Algorithm::SortByKey(uniqueIdsWithContourId, this->InterpolationWeights);
}
else
{
Algorithm::SortByKey(uniqueIds, this->InterpolationWeights);
}
//3. lastly we need to do a unique by key, but since vtkm doesn't
// offer that feature, we use a zip handle.
@ -649,7 +707,7 @@ vtkm::cont::CellSetSingleType< >
vtkm::cont::ArrayHandleZip<Id2HandleType, WeightHandleType> zipped =
vtkm::cont::make_ArrayHandleZip(uniqueIds, this->InterpolationWeights);
Algorithm::Unique( zipped, marchingcubes::FirstValueSame());
}
}
//4.
//LowerBounds generates the output cell connections. It does this by
@ -670,11 +728,8 @@ 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 use a CellSetSingleType with explicit
//storage;
{
vtkm::cont::ArrayHandleIndex temp(this->InterpolationIds.GetNumberOfValues());
Algorithm::Copy(temp, connectivity);
}
}
//generate the vertices's

7
vtkm/worklet/testing/UnitTestMarchingCubes.cxx
View File

@ -251,10 +251,12 @@ void TestMarchingCubesUniformGrid()
vtkm::worklet::MarchingCubes isosurfaceFilter;
isosurfaceFilter.SetMergeDuplicatePoints(false);
vtkm::Float32 contourValue = 0.5f;
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float32,3> > verticesArray;
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::Float32,3> > normalsArray;
vtkm::cont::ArrayHandle<vtkm::Float32> scalarsArray;
isosurfaceFilter.Run(0.5f,
isosurfaceFilter.Run(&contourValue,
1,
cellSet,
dataSet.GetCoordinateSystem(),
fieldArray,
@ -309,7 +311,8 @@ void TestMarchingCubesExplicit()
vtkm::worklet::MarchingCubes marchingCubes;
marchingCubes.SetMergeDuplicatePoints(false);
marchingCubes.Run(contourValue,
marchingCubes.Run(&contourValue,
1,
cellSet,
dataSet.GetCoordinateSystem(),
contourArray,