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

Significantly improve FlyingEdges performance across all devices

Browse Source 
We now use SumYAxis when executing with CUDA for better memory patterns.
Instead of using the heavy Pass4/Pass4WithNormals, CUDA now uses a
2 pass approach with the second pass outputting the normals and
coordinates using with significantly less warp divergence
This commit is contained in:
Robert Maynard 2020-05-01 13:03:36 -04:00
parent fa93738010
commit 251bd82b80
13 changed files with 1290 additions and 612 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
vtkm/cont/ArrayHandleView.h
View File

@ -274,7 +274,7 @@ public:
{
throw vtkm::cont::ErrorBadValue(
"The value array must be pre-allocated before it is used for the "
"output of ArrayHandlePermutation.");
"output of ArrayHandleView.");
}
return PortalExecution(

70
vtkm/filter/testing/UnitTestContourFilterNormals.cxx
View File

@ -17,6 +17,8 @@
#include <vtkm/filter/Contour.h>
#include <vtkm/io/VTKDataSetWriter.h>
namespace vtkm_ut_mc_normals
{
@ -70,6 +72,10 @@ void TestNormals(const vtkm::cont::DataSet& dataset, bool structured)
{ 0.770536f, -0.421248f, -0.478356f }, { -0.736036f, -0.445244f, -0.509910f },
{ 0.123446f, -0.887088f, -0.444788f }, { 0.133328f, -0.397444f, -0.907889f }
};
//Calculated using FlyingEdges and Y axis iteration which causes
//the points to be in a different order
const vtkm::Id fe_y_alg_ordering[numVerts] = { 0, 1, 3, 5, 4, 6, 2, 7,
9, 12, 10, 13, 8, 14, 11, 15 };
//Calculated using normals of the output triangles
const vtkm::Vec3f fast[numVerts] = {
@ -83,6 +89,19 @@ void TestNormals(const vtkm::cont::DataSet& dataset, bool structured)
{ 0.2164f, -0.9401f, -0.2635f }, { -0.1589f, -0.1642f, -0.9735f }
};
//When using the Y axis algorithm the cells are generated in a different
//order.
const vtkm::Vec3f fast_fe_y[numVerts] = {
{ -0.243433f, -0.429741f, -0.869519f }, { 0.158904f, 0.164214f, -0.973542f },
{ -0.895292f, -0.390217f, -0.214903f }, { -0.895057f, 0.134692f, -0.425125f },
{ 0.829547f, -0.418793f, -0.36941f }, { 0.846705f, 0.425787f, -0.319054f },
{ 0.253811f, -0.853394f, -0.4553f }, { -0.216381f, 0.940084f, -0.263478f },
{ -0.848579f, -0.35602f, 0.391362f }, { -0.93948f, 0.252957f, 0.231065f },
{ 0.831549f, -0.472663f, 0.291744f }, { 0.910494f, 0.0298277f, 0.412446f },
{ -0.362862f, -0.815464f, 0.450944f }, { 0.107848f, 0.958544f, 0.263748f },
{ 0.135131f, -0.437674f, 0.888921f }, { -0.286251f, 0.172078f, 0.942576f }
};
vtkm::cont::ArrayHandle<vtkm::Vec3f> normals;
vtkm::filter::Contour mc;
@ -97,16 +116,23 @@ void TestNormals(const vtkm::cont::DataSet& dataset, bool structured)
result.GetField("normals").GetData().CopyTo(normals);
VTKM_TEST_ASSERT(normals.GetNumberOfValues() == numVerts,
"Wrong number of values in normals field");
auto normalsPortal = normals.ReadPortal();
for (vtkm::Id i = 0; i < numVerts; ++i)
//determine if we are using flying edge Y axis algorithm by checking the first normal value that differs
const bool using_fe_y_alg_ordering =
test_equal(normals.ReadPortal().Get(2), expected[fe_y_alg_ordering[2]], 0.001);
{
VTKM_TEST_ASSERT(test_equal(normalsPortal.Get(i), expected[i], 0.001),
"Result (",
normalsPortal.Get(i),
") does not match expected value (",
expected[i],
") vert ",
i);
auto normalPotals = normals.ReadPortal();
for (vtkm::Id i = 0; i < numVerts; ++i)
{
auto expected_v = !using_fe_y_alg_ordering ? expected[i] : expected[fe_y_alg_ordering[i]];
VTKM_TEST_ASSERT(test_equal(normalPotals.Get(i), expected_v, 0.001),
"Result (",
normalPotals.Get(i),
") does not match expected value (",
expected_v,
") vert ",
i);
}
}
// Test the other normals generation method
@ -114,6 +140,10 @@ void TestNormals(const vtkm::cont::DataSet& dataset, bool structured)
{
mc.SetComputeFastNormalsForStructured(true);
expected = fast;
if (using_fe_y_alg_ordering)
{
expected = fast_fe_y;
}
}
else
{
@ -125,16 +155,20 @@ void TestNormals(const vtkm::cont::DataSet& dataset, bool structured)
result.GetField("normals").GetData().CopyTo(normals);
VTKM_TEST_ASSERT(normals.GetNumberOfValues() == numVerts,
"Wrong number of values in normals field");
normalsPortal = normals.ReadPortal();
for (vtkm::Id i = 0; i < numVerts; ++i)
{
VTKM_TEST_ASSERT(test_equal(normalsPortal.Get(i), expected[i], 0.001),
"Result (",
normalsPortal.Get(i),
") does not match expected value (",
expected[i],
") vert ",
i);
auto normalPotals = normals.ReadPortal();
for (vtkm::Id i = 0; i < numVerts; ++i)
{
bool equal = test_equal(normalPotals.Get(i), expected[i], 0.001);
VTKM_TEST_ASSERT(equal,
"Result (",
normalPotals.Get(i),
") does not match expected value (",
expected[i],
") vert ",
i);
}
}
}

5
vtkm/worklet/contour/CMakeLists.txt
View File

@ -16,7 +16,10 @@ set(headers
FlyingEdgesPass1.h
FlyingEdgesPass2.h
FlyingEdgesPass4.h
FlyingEdgesPass4WithNormals.h
FlyingEdgesPass4Common.h
FlyingEdgesPass4X.h
FlyingEdgesPass4XWithNormals.h
FlyingEdgesPass4Y.h
FlyingEdgesTables.h
MarchingCellTables.h
MarchingCells.h

97
vtkm/worklet/contour/FlyingEdges.h
View File

@ -16,7 +16,6 @@
#include <vtkm/worklet/contour/FlyingEdgesPass1.h>
#include <vtkm/worklet/contour/FlyingEdgesPass2.h>
#include <vtkm/worklet/contour/FlyingEdgesPass4.h>
#include <vtkm/worklet/contour/FlyingEdgesPass4WithNormals.h>
#include <vtkm/cont/ArrayHandleGroupVec.h>
#include <vtkm/cont/Invoker.h>
@ -30,21 +29,6 @@ namespace flying_edges
namespace detail
{
inline vtkm::cont::CellSetStructured<2> make_metaDataMesh2D(SumXAxis, const vtkm::Id3& pdims)
{
vtkm::cont::CellSetStructured<2> metaDataMesh;
metaDataMesh.SetPointDimensions(vtkm::Id2{ pdims[1], pdims[2] });
return metaDataMesh;
}
inline vtkm::cont::CellSetStructured<2> make_metaDataMesh2D(SumYAxis, const vtkm::Id3& pdims)
{
vtkm::cont::CellSetStructured<2> metaDataMesh;
metaDataMesh.SetPointDimensions(vtkm::Id2{ pdims[0], pdims[2] });
return metaDataMesh;
}
template <typename T, typename S>
vtkm::Id extend_by(vtkm::cont::ArrayHandle<T, S>& handle, vtkm::Id size)
{
@ -64,7 +48,6 @@ vtkm::Id extend_by(vtkm::cont::ArrayHandle<T, S>& handle, vtkm::Id size)
}
}
//----------------------------------------------------------------------------
template <typename ValueType,
typename StorageTagField,
@ -81,12 +64,6 @@ vtkm::cont::CellSetSingleType<> execute(
vtkm::cont::ArrayHandle<vtkm::Vec<NormalType, 3>, StorageTagNormals>& normals,
vtkm::worklet::contour::CommonState& sharedState)
{
//Tasks:
//3. Support switching AxisToSum by running this whole thing in a TryExecute
// Passes 5 can ignore this
using AxisToSum = SumXAxis;
vtkm::cont::Invoker invoke;
vtkm::Vec3f origin, spacing;
@ -98,12 +75,10 @@ vtkm::cont::CellSetSingleType<> execute(
}
auto pdims = cells.GetPointDimensions();
vtkm::cont::ArrayHandle<vtkm::UInt8> edgeCases;
edgeCases.Allocate(coordinateSystem.GetNumberOfValues());
vtkm::cont::CellSetStructured<2> metaDataMesh2D = detail::make_metaDataMesh2D(AxisToSum{}, pdims);
vtkm::cont::CellSetStructured<2> metaDataMesh2D;
vtkm::cont::ArrayHandle<vtkm::Id> metaDataLinearSums; //per point of metaDataMesh
vtkm::cont::ArrayHandle<vtkm::Id> metaDataMin; //per point of metaDataMesh
vtkm::cont::ArrayHandle<vtkm::Id> metaDataMax; //per point of metaDataMesh
@ -141,7 +116,7 @@ vtkm::cont::CellSetSingleType<> execute(
// Additionally CUDA does significantly better when you do an initial fill
// and write only non-below values
//
ComputePass1<ValueType, AxisToSum> worklet1(isoval, pdims);
ComputePass1<ValueType> worklet1(isoval, pdims);
vtkm::cont::TryExecuteOnDevice(invoke.GetDevice(),
launchComputePass1{},
worklet1,
@ -160,7 +135,7 @@ vtkm::cont::CellSetSingleType<> execute(
// row. Use computational trimming to reduce work.
{
VTKM_LOG_SCOPE(vtkm::cont::LogLevel::Perf, "FlyingEdges Pass2");
ComputePass2<AxisToSum> worklet2(pdims);
ComputePass2 worklet2(pdims);
invoke(worklet2,
metaDataMesh2D,
metaDataSums,
@ -182,56 +157,40 @@ vtkm::cont::CellSetSingleType<> execute(
detail::extend_by(sharedState.CellIdMap, sumTris);
auto newPointSize =
vtkm::Id newPointSize =
vtkm::cont::Algorithm::ScanExclusive(metaDataLinearSums, metaDataLinearSums);
detail::extend_by(sharedState.InterpolationEdgeIds, newPointSize);
detail::extend_by(sharedState.InterpolationWeights, newPointSize);
//----------------------------------------------------------------------------
// PASS 4: Process voxel rows and generate topology, and interpolation state
if (sharedState.GenerateNormals)
{
VTKM_LOG_SCOPE(vtkm::cont::LogLevel::Perf, "FlyingEdges ComputePass4WithNormals");
detail::extend_by(points, newPointSize);
detail::extend_by(normals, newPointSize);
VTKM_LOG_SCOPE(vtkm::cont::LogLevel::Perf, "FlyingEdges Pass4");
ComputePass4WithNormals<ValueType, AxisToSum> worklet4(
isoval, pdims, origin, spacing, multiContourCellOffset, multiContourPointOffset);
invoke(worklet4,
metaDataMesh2D,
metaDataSums,
metaDataMin,
metaDataMax,
metaDataNumTris,
edgeCases,
inputField,
triangle_topology,
sharedState.InterpolationEdgeIds,
sharedState.InterpolationWeights,
sharedState.CellIdMap,
points,
normals);
}
else
{
VTKM_LOG_SCOPE(vtkm::cont::LogLevel::Perf, "FlyingEdges ComputePass4");
detail::extend_by(points, newPointSize);
launchComputePass4 pass4(
pdims, origin, spacing, multiContourCellOffset, multiContourPointOffset);
ComputePass4<ValueType, AxisToSum> worklet4(
isoval, pdims, origin, spacing, multiContourCellOffset, multiContourPointOffset);
invoke(worklet4,
metaDataMesh2D,
metaDataSums,
metaDataMin,
metaDataMax,
metaDataNumTris,
edgeCases,
inputField,
triangle_topology,
sharedState.InterpolationEdgeIds,
sharedState.InterpolationWeights,
sharedState.CellIdMap,
points);
detail::extend_by(points, newPointSize);
if (sharedState.GenerateNormals)
{
detail::extend_by(normals, newPointSize);
}
vtkm::cont::TryExecuteOnDevice(invoke.GetDevice(),
pass4,
newPointSize,
isoval,
inputField,
edgeCases,
metaDataMesh2D,
metaDataSums,
metaDataMin,
metaDataMax,
metaDataNumTris,
sharedState,
triangle_topology,
points,
normals);
}
}
}

26
vtkm/worklet/contour/FlyingEdgesHelpers.h
View File

@ -55,6 +55,32 @@ struct SumYAxis
static constexpr vtkm::Id zindex = 2;
};
template <typename Device>
struct select_AxisToSum
{
using type = SumXAxis;
};
template <>
struct select_AxisToSum<vtkm::cont::DeviceAdapterTagCuda>
{
using type = SumYAxis;
};
inline vtkm::cont::CellSetStructured<2> make_metaDataMesh2D(SumXAxis, const vtkm::Id3& pdims)
{
vtkm::cont::CellSetStructured<2> metaDataMesh;
metaDataMesh.SetPointDimensions(vtkm::Id2{ pdims[1], pdims[2] });
return metaDataMesh;
}
inline vtkm::cont::CellSetStructured<2> make_metaDataMesh2D(SumYAxis, const vtkm::Id3& pdims)
{
vtkm::cont::CellSetStructured<2> metaDataMesh;
metaDataMesh.SetPointDimensions(vtkm::Id2{ pdims[0], pdims[2] });
return metaDataMesh;
}
VTKM_EXEC inline vtkm::Id3 compute_ijk(SumXAxis, const vtkm::Id3& executionSpaceIJK)
{

26
vtkm/worklet/contour/FlyingEdgesPass1.h
View File

@ -65,7 +65,7 @@ inline VTKM_EXEC void write_edge(vtkm::cont::DeviceAdapterTagCuda,
}
}
template <typename T, typename AxisToSum>
template <typename T>
struct ComputePass1 : public vtkm::worklet::WorkletVisitPointsWithCells
{
vtkm::Id3 PointDims;
@ -99,6 +99,8 @@ struct ComputePass1 : public vtkm::worklet::WorkletVisitPointsWithCells
const WholeDataField& field,
Device device) const
{
using AxisToSum = typename select_AxisToSum<Device>::type;
const vtkm::Id3 ijk = compute_ijk(AxisToSum{}, threadIndices.GetInputIndex3D());
const vtkm::Id3 dims = this->PointDims;
const vtkm::Id startPos = compute_start(AxisToSum{}, ijk, dims);
@ -143,32 +145,34 @@ struct ComputePass1 : public vtkm::worklet::WorkletVisitPointsWithCells
struct launchComputePass1
{
template <typename DeviceAdapterTag,
typename T,
typename AxisType,
typename StorageTagField,
typename... Args>
template <typename DeviceAdapterTag, typename T, typename StorageTagField, typename... Args>
VTKM_CONT bool operator()(DeviceAdapterTag device,
const ComputePass1<T, AxisType>& worklet,
const ComputePass1<T>& worklet,
const vtkm::cont::ArrayHandle<T, StorageTagField>& inputField,
vtkm::cont::ArrayHandle<vtkm::UInt8> edgeCases,
vtkm::cont::CellSetStructured<2>& metaDataMesh2D,
Args&&... args) const
{
vtkm::cont::Invoker invoke(device);
invoke(worklet, std::forward<Args>(args)..., edgeCases, inputField);
metaDataMesh2D = make_metaDataMesh2D(SumXAxis{}, worklet.PointDims);
invoke(worklet, metaDataMesh2D, std::forward<Args>(args)..., edgeCases, inputField);
return true;
}
template <typename T, typename AxisType, typename StorageTagField, typename... Args>
template <typename T, typename StorageTagField, typename... Args>
VTKM_CONT bool operator()(vtkm::cont::DeviceAdapterTagCuda device,
const ComputePass1<T, AxisType>& worklet,
const ComputePass1<T>& worklet,
const vtkm::cont::ArrayHandle<T, StorageTagField>& inputField,
vtkm::cont::ArrayHandle<vtkm::UInt8> edgeCases,
vtkm::cont::CellSetStructured<2>& metaDataMesh2D,
Args&&... args) const
{
vtkm::cont::Invoker invoke(device);
metaDataMesh2D = make_metaDataMesh2D(SumYAxis{}, worklet.PointDims);
vtkm::cont::Algorithm::Fill(edgeCases, static_cast<vtkm::UInt8>(FlyingEdges3D::Below));
invoke(worklet, std::forward<Args>(args)..., edgeCases, inputField);
invoke(worklet, metaDataMesh2D, std::forward<Args>(args)..., edgeCases, inputField);
return true;
}
};

18
vtkm/worklet/contour/FlyingEdgesPass2.h
View File

@ -23,7 +23,6 @@ namespace worklet
namespace flying_edges
{
template <typename AxisToSum>
struct ComputePass2 : public vtkm::worklet::WorkletVisitCellsWithPoints
{
vtkm::Id3 PointDims;
@ -40,20 +39,24 @@ struct ComputePass2 : public vtkm::worklet::WorkletVisitCellsWithPoints
FieldInPoint axis_maxs,
FieldOutCell cell_tri_count,
WholeArrayIn edgeData);
using ExecutionSignature = void(ThreadIndices, _2, _3, _4, _5, _6);
using ExecutionSignature = void(ThreadIndices, _2, _3, _4, _5, _6, Device);
using InputDomain = _1;
template <typename ThreadIndices,
typename WholeSumField,
typename FieldInPointId,
typename WholeEdgeField>
typename WholeEdgeField,
typename Device>
VTKM_EXEC void operator()(const ThreadIndices& threadIndices,
const WholeSumField& axis_sums,
const FieldInPointId& axis_mins,
const FieldInPointId& axis_maxs,
vtkm::Int32& cell_tri_count,
const WholeEdgeField& edges) const
const WholeEdgeField& edges,
Device) const
{
using AxisToSum = typename select_AxisToSum<Device>::type;
// Pass 2. Traverse all cells in the meta data plane. This allows us to
// easily grab the four edge cases bounding this voxel-row
@ -134,7 +137,8 @@ struct ComputePass2 : public vtkm::worklet::WorkletVisitCellsWithPoints
sums[AxisToSum::zindex] += edgeUses[8];
// handle boundary
this->CountBoundaryEdgeUses(onBoundary, edgeUses, sums, adj_row_sum, adj_col_sum);
this->CountBoundaryEdgeUses(
AxisToSum{}, onBoundary, edgeUses, sums, adj_row_sum, adj_col_sum);
}
}
@ -157,7 +161,9 @@ struct ComputePass2 : public vtkm::worklet::WorkletVisitCellsWithPoints
//
// Only on these boundaries do we write to the metaData of our neighbor
// as it is safe as those
VTKM_EXEC inline void CountBoundaryEdgeUses(vtkm::Vec<bool, 3> onBoundary,
template <typename AxisToSum>
VTKM_EXEC inline void CountBoundaryEdgeUses(AxisToSum,
vtkm::Vec<bool, 3> onBoundary,
vtkm::UInt8 const* const edgeUses,
vtkm::Id3& sums,
vtkm::Id3& adj_row_sum,

612
vtkm/worklet/contour/FlyingEdgesPass4.h
View File

@ -13,9 +13,10 @@
#ifndef vtk_m_worklet_contour_flyingedges_pass4_h
#define vtk_m_worklet_contour_flyingedges_pass4_h
#include <vtkm/worklet/contour/FlyingEdgesHelpers.h>
#include <vtkm/worklet/contour/FlyingEdgesTables.h>
#include <vtkm/worklet/contour/FlyingEdgesPass4Common.h>
#include <vtkm/worklet/contour/FlyingEdgesPass4X.h>
#include <vtkm/worklet/contour/FlyingEdgesPass4XWithNormals.h>
#include <vtkm/worklet/contour/FlyingEdgesPass4Y.h>
namespace vtkm
{
@ -24,519 +25,152 @@ namespace worklet
namespace flying_edges
{
VTKM_EXEC inline vtkm::Id3 compute_incs3d(const vtkm::Id3& dims)
struct launchComputePass4
{
return vtkm::Id3{ 1, dims[0], (dims[0] * dims[1]) };
}
VTKM_EXEC inline bool case_includes_axes(vtkm::UInt8 const* const edgeUses)
{
return (edgeUses[0] != 0 || edgeUses[4] != 0 || edgeUses[8] != 0);
}
template <typename WholeConnField, typename WholeCellIdField>
VTKM_EXEC inline void generate_tris(vtkm::Id inputCellId,
vtkm::UInt8 edgeCase,
vtkm::UInt8 numTris,
vtkm::Id* edgeIds,
vtkm::Id& triId,
const WholeConnField& conn,
const WholeCellIdField& cellIds)
{
auto* edges = data::GetTriEdgeCases(edgeCase);
vtkm::Id edgeIndex = 1;
vtkm::Id index = static_cast<vtkm::Id>(triId) * 3;
for (vtkm::UInt8 i = 0; i < numTris; ++i)
{
cellIds.Set(triId + i, inputCellId);
//We use edgeIndex, edgeIndex+2, edgeIndex+1 to keep
//the same winding for the triangles that marching celss
//produced. By keeping the winding the same we make sure
//that 'fast' normals are consistent with the marching
//cells version
conn.Set(index, edgeIds[edges[edgeIndex]]);
conn.Set(index + 1, edgeIds[edges[edgeIndex + 2]]);
conn.Set(index + 2, edgeIds[edges[edgeIndex + 1]]);
index += 3;
edgeIndex += 3;
}
triId += numTris;
}
// Helper function to set up the point ids on voxel edges.
//----------------------------------------------------------------------------
template <typename AxisToSum, typename FieldInPointId3>
VTKM_EXEC inline void init_voxelIds(AxisToSum,
vtkm::Id writeOffset,
vtkm::UInt8 edgeCase,
const FieldInPointId3& axis_sums,
vtkm::Id* edgeIds)
{
auto* edgeUses = data::GetEdgeUses(edgeCase);
edgeIds[0] = writeOffset + axis_sums[0][AxisToSum::xindex]; // x-edges
edgeIds[1] = writeOffset + axis_sums[1][AxisToSum::xindex];
edgeIds[2] = writeOffset + axis_sums[3][AxisToSum::xindex];
edgeIds[3] = writeOffset + axis_sums[2][AxisToSum::xindex];
edgeIds[4] = writeOffset + axis_sums[0][AxisToSum::yindex]; // y-edges
edgeIds[5] = edgeIds[4] + edgeUses[4];
edgeIds[6] = writeOffset + axis_sums[3][AxisToSum::yindex];
edgeIds[7] = edgeIds[6] + edgeUses[6];
edgeIds[8] = writeOffset + axis_sums[0][AxisToSum::zindex]; // z-edges
edgeIds[9] = edgeIds[8] + edgeUses[8];
edgeIds[10] = writeOffset + axis_sums[1][AxisToSum::zindex];
edgeIds[11] = edgeIds[10] + edgeUses[10];
}
// Helper function to advance the point ids along voxel rows.
//----------------------------------------------------------------------------
VTKM_EXEC inline void advance_voxelIds(vtkm::UInt8 const* const edgeUses, vtkm::Id* edgeIds)
{
edgeIds[0] += edgeUses[0]; // x-edges
edgeIds[1] += edgeUses[1];
edgeIds[2] += edgeUses[2];
edgeIds[3] += edgeUses[3];
edgeIds[4] += edgeUses[4]; // y-edges
edgeIds[5] = edgeIds[4] + edgeUses[5];
edgeIds[6] += edgeUses[6];
edgeIds[7] = edgeIds[6] + edgeUses[7];
edgeIds[8] += edgeUses[8]; // z-edges
edgeIds[9] = edgeIds[8] + edgeUses[9];
edgeIds[10] += edgeUses[10];
edgeIds[11] = edgeIds[10] + edgeUses[11];
}
//----------------------------------------------------------------------------
struct Pass4TrimState
{
vtkm::Id left, right;
vtkm::Id3 ijk;
vtkm::Id4 startPos;
vtkm::Id axis_inc;
vtkm::UInt8 yzLoc;
bool valid = true;
template <typename AxisToSum,
typename ThreadIndices,
typename FieldInPointId,
typename WholeEdgeField>
VTKM_EXEC Pass4TrimState(AxisToSum,
const vtkm::Id3& pdims,
const ThreadIndices& threadIndices,
const FieldInPointId& axis_mins,
const FieldInPointId& axis_maxs,
const WholeEdgeField& edges)
{
// find adjusted trim values.
left = vtkm::Min(axis_mins[0], axis_mins[1]);
left = vtkm::Min(left, axis_mins[2]);
left = vtkm::Min(left, axis_mins[3]);
right = vtkm::Max(axis_maxs[0], axis_maxs[1]);
right = vtkm::Max(right, axis_maxs[2]);
right = vtkm::Max(right, axis_maxs[3]);
ijk = compute_ijk(AxisToSum{}, threadIndices.GetInputIndex3D());
startPos = compute_neighbor_starts(AxisToSum{}, ijk, pdims);
axis_inc = compute_inc(AxisToSum{}, pdims);
if (left == pdims[AxisToSum::xindex] && right == 0)
{
//verify that we have nothing to generate and early terminate.
bool mins_same = (axis_mins[0] == axis_mins[1] && axis_mins[0] == axis_mins[2] &&
axis_mins[0] == axis_mins[3]);
bool maxs_same = (axis_maxs[0] == axis_maxs[1] && axis_maxs[0] == axis_maxs[2] &&
axis_maxs[0] == axis_maxs[3]);
if (mins_same && maxs_same)
{
valid = false;
return;
}
else
{
left = 0;
right = pdims[AxisToSum::xindex] - 1;
}
}
// The trim edges may need adjustment if the contour travels between rows
// of edges (without intersecting these edges). This means checking
// whether the trim faces at (left,right) made up of the edges intersect
// the contour.
adjustTrimBounds(pdims[AxisToSum::xindex] - 1, edges, startPos, axis_inc, left, right);
if (left == right)
{
valid = false;
return;
}
const vtkm::UInt8 yLoc =
(ijk[AxisToSum::yindex] < 1
? FlyingEdges3D::MinBoundary
: (ijk[AxisToSum::yindex] >= (pdims[AxisToSum::yindex] - 2) ? FlyingEdges3D::MaxBoundary
: FlyingEdges3D::Interior));
const vtkm::UInt8 zLoc =
(ijk[AxisToSum::zindex] < 1
? FlyingEdges3D::MinBoundary
: (ijk[AxisToSum::zindex] >= (pdims[AxisToSum::zindex] - 2) ? FlyingEdges3D::MaxBoundary
: FlyingEdges3D::Interior));
yzLoc = static_cast<vtkm::UInt8>((yLoc << 2) | (zLoc << 4));
}
};
template <typename T, typename AxisToSum>
struct ComputePass4 : public vtkm::worklet::WorkletVisitCellsWithPoints
{
vtkm::Id3 PointDims;
vtkm::Vec3f Origin;
vtkm::Vec3f Spacing;
T IsoValue;
vtkm::Id CellWriteOffset;
vtkm::Id PointWriteOffset;
ComputePass4() {}
ComputePass4(T value,
const vtkm::Id3& pdims,
const vtkm::Vec3f& origin,
const vtkm::Vec3f& spacing,
vtkm::Id multiContourCellOffset,
vtkm::Id multiContourPointOffset)
launchComputePass4(const vtkm::Id3& pdims,
const vtkm::Vec3f& origin,
const vtkm::Vec3f& spacing,
vtkm::Id multiContourCellOffset,
vtkm::Id multiContourPointOffset)
: PointDims(pdims)
, Origin(origin)
, Spacing(spacing)
, IsoValue(value)
, CellWriteOffset(multiContourCellOffset)
, PointWriteOffset(multiContourPointOffset)
{
}
using ControlSignature = void(CellSetIn,
FieldInPoint axis_sums,
FieldInPoint axis_mins,
FieldInPoint axis_maxs,
WholeArrayIn cell_tri_count,
WholeArrayIn edgeData,
WholeArrayIn data,
WholeArrayOut connectivity,
WholeArrayOut edgeIds,
WholeArrayOut weights,
WholeArrayOut points,
WholeArrayOut inputCellIds);
using ExecutionSignature =
void(ThreadIndices, _2, _3, _4, _5, _6, _7, _8, _9, _10, _11, _12, WorkIndex);
template <typename ThreadIndices,
typename FieldInPointId3,
typename FieldInPointId,
typename WholeTriField,
typename WholeEdgeField,
typename WholeDataField,
typename WholeConnField,
typename WholeEdgeIdField,
typename WholeWeightField,
typename WholeCellIdField,
typename WholePointField>
VTKM_EXEC void operator()(const ThreadIndices& threadIndices,
const FieldInPointId3& axis_sums,
const FieldInPointId& axis_mins,
const FieldInPointId& axis_maxs,
const WholeTriField& cellTriCount,
const WholeEdgeField& edges,
const WholeDataField& field,
const WholeConnField& conn,
const WholeEdgeIdField& interpolatedEdgeIds,
const WholeWeightField& weights,
const WholeCellIdField& inputCellIds,
const WholePointField& points,
vtkm::Id oidx) const
template <typename DeviceAdapterTag,
typename T,
typename StorageTagField,
typename MeshSums,
typename PointType,
typename NormalType>
VTKM_CONT bool operator()(DeviceAdapterTag device,
vtkm::Id vtkmNotUsed(newPointSize),
T isoval,
const vtkm::cont::ArrayHandle<T, StorageTagField>& inputField,
vtkm::cont::ArrayHandle<vtkm::UInt8> edgeCases,
vtkm::cont::CellSetStructured<2>& metaDataMesh2D,
const MeshSums& metaDataSums,
const vtkm::cont::ArrayHandle<vtkm::Id>& metaDataMin,
const vtkm::cont::ArrayHandle<vtkm::Id>& metaDataMax,
const vtkm::cont::ArrayHandle<vtkm::Int32>& metaDataNumTris,
vtkm::worklet::contour::CommonState& sharedState,
vtkm::cont::ArrayHandle<vtkm::Id>& triangle_topology,
PointType& points,
NormalType& normals) const
{
//This works as cellTriCount was computed with ScanExtended
//and therefore has one more entry than the number of cells
vtkm::Id cell_tri_offset = cellTriCount.Get(oidx);
vtkm::Id next_tri_offset = cellTriCount.Get(oidx + 1);
if (cell_tri_offset == next_tri_offset)
{ //we produce nothing
return;
}
cell_tri_offset += this->CellWriteOffset;
const Pass4TrimState state(
AxisToSum{}, this->PointDims, threadIndices, axis_mins, axis_maxs, edges);
if (!state.valid)
vtkm::cont::Invoker invoke(device);
if (sharedState.GenerateNormals)
{
return;
ComputePass4XWithNormals<T> worklet4(isoval,
this->PointDims,
this->Origin,
this->Spacing,
this->CellWriteOffset,
this->PointWriteOffset);
invoke(worklet4,
metaDataMesh2D,
metaDataSums,
metaDataMin,
metaDataMax,
metaDataNumTris,
edgeCases,
inputField,
triangle_topology,
sharedState.InterpolationEdgeIds,
sharedState.InterpolationWeights,
sharedState.CellIdMap,
points,
normals);
}
const vtkm::Id3 pdims = this->PointDims;
const vtkm::Id3 increments = compute_incs3d(pdims);
vtkm::Id edgeIds[12];
auto edgeCase = getEdgeCase(edges, state.startPos, (state.axis_inc * state.left));
init_voxelIds(AxisToSum{}, this->PointWriteOffset, edgeCase, axis_sums, edgeIds);
for (vtkm::Id i = state.left; i < state.right; ++i) // run along the trimmed voxels
else
{
auto ijk = state.ijk;
ijk[AxisToSum::xindex] = i;
edgeCase = getEdgeCase(edges, state.startPos, (state.axis_inc * i));
vtkm::UInt8 numTris = data::GetNumberOfPrimitives(edgeCase);
if (numTris > 0)
{
//compute what the current cellId is
vtkm::Id cellId = compute_start(AxisToSum{}, ijk, pdims - vtkm::Id3{ 1, 1, 1 });
// Start by generating triangles for this case
generate_tris(cellId, edgeCase, numTris, edgeIds, cell_tri_offset, conn, inputCellIds);
// Now generate edgeIds and weights along voxel axes if needed. Remember to take
// boundary into account.
vtkm::UInt8 loc = static_cast<vtkm::UInt8>(
state.yzLoc | (i < 1 ? FlyingEdges3D::MinBoundary
: (i >= (pdims[AxisToSum::xindex] - 2) ? FlyingEdges3D::MaxBoundary
: FlyingEdges3D::Interior)));
auto* edgeUses = data::GetEdgeUses(edgeCase);
if (loc != FlyingEdges3D::Interior || case_includes_axes(edgeUses))
{
this->Generate(loc,
ijk,
field,
interpolatedEdgeIds,
weights,
points,
state.startPos,
increments,
(state.axis_inc * i),
edgeUses,
edgeIds);
}
advance_voxelIds(edgeUses, edgeIds);
}
ComputePass4X<T> worklet4(isoval,
this->PointDims,
this->Origin,
this->Spacing,
this->CellWriteOffset,
this->PointWriteOffset);
invoke(worklet4,
metaDataMesh2D,
metaDataSums,
metaDataMin,
metaDataMax,
metaDataNumTris,
edgeCases,
inputField,
triangle_topology,
sharedState.InterpolationEdgeIds,
sharedState.InterpolationWeights,
sharedState.CellIdMap,
points);
}
return true;
}
//----------------------------------------------------------------------------
template <typename WholeDataField,
typename WholeIEdgeField,
typename WholeWeightField,
typename WholePointField>
VTKM_EXEC inline void Generate(vtkm::UInt8 loc,
const vtkm::Id3& ijk,
const WholeDataField& field,
const WholeIEdgeField& interpolatedEdgeIds,
const WholeWeightField& weights,
const WholePointField& points,
const vtkm::Id4& startPos,
const vtkm::Id3& incs,
vtkm::Id offset,
vtkm::UInt8 const* const edgeUses,
vtkm::Id* edgeIds) const
template <typename T,
typename StorageTagField,
typename MeshSums,
typename PointType,
typename NormalType>
VTKM_CONT bool operator()(vtkm::cont::DeviceAdapterTagCuda device,
vtkm::Id newPointSize,
T isoval,
const vtkm::cont::ArrayHandle<T, StorageTagField>& inputField,
vtkm::cont::ArrayHandle<vtkm::UInt8> edgeCases,
vtkm::cont::CellSetStructured<2>& metaDataMesh2D,
const MeshSums& metaDataSums,
const vtkm::cont::ArrayHandle<vtkm::Id>& metaDataMin,
const vtkm::cont::ArrayHandle<vtkm::Id>& metaDataMax,
const vtkm::cont::ArrayHandle<vtkm::Int32>& metaDataNumTris,
vtkm::worklet::contour::CommonState& sharedState,
vtkm::cont::ArrayHandle<vtkm::Id>& triangle_topology,
PointType& points,
NormalType& normals) const
{
vtkm::Id2 pos(startPos[0] + offset, 0);
{
auto s0 = field.Get(pos[0]);
vtkm::cont::Invoker invoke(device);
//EdgesUses 0,4,8 work for Y axis
if (edgeUses[0])
{ // edgesUses[0] == i axes edge
auto writeIndex = edgeIds[0];
pos[1] = startPos[0] + offset + incs[AxisToSum::xindex];
auto s1 = field.Get(pos[1]);
auto t = (this->IsoValue - s0) / (s1 - s0);
ComputePass4Y<T> worklet4(
isoval, this->PointDims, this->CellWriteOffset, this->PointWriteOffset);
invoke(worklet4,
metaDataMesh2D,
metaDataSums,
metaDataMin,
metaDataMax,
metaDataNumTris,
edgeCases,
inputField,
triangle_topology,
sharedState.InterpolationEdgeIds,
sharedState.InterpolationWeights,
sharedState.CellIdMap);
interpolatedEdgeIds.Set(writeIndex, pos);
weights.Set(writeIndex, static_cast<vtkm::FloatDefault>(t));
//This needs to be done on array handle view ( start = this->PointWriteOffset, len = newPointSize)
ComputePass5Y<T> worklet5(this->PointDims,
this->Origin,
this->Spacing,
this->PointWriteOffset,
sharedState.GenerateNormals);
invoke(worklet5,
vtkm::cont::make_ArrayHandleView(
sharedState.InterpolationEdgeIds, this->PointWriteOffset, newPointSize),
vtkm::cont::make_ArrayHandleView(
sharedState.InterpolationWeights, this->PointWriteOffset, newPointSize),
vtkm::cont::make_ArrayHandleView(points, this->PointWriteOffset, newPointSize),
inputField,
normals);
auto coord = this->InterpolateCoordinate(t, ijk, ijk + vtkm::Id3{ 1, 0, 0 });
points.Set(writeIndex, coord);
}
if (edgeUses[4])
{ // edgesUses[4] == j axes edge
auto writeIndex = edgeIds[4];
pos[1] = startPos[1] + offset;
auto s1 = field.Get(pos[1]);
auto t = (this->IsoValue - s0) / (s1 - s0);
interpolatedEdgeIds.Set(writeIndex, pos);
weights.Set(writeIndex, static_cast<vtkm::FloatDefault>(t));
auto coord = this->InterpolateCoordinate(t, ijk, ijk + vtkm::Id3{ 0, 1, 0 });
points.Set(writeIndex, coord);
}
if (edgeUses[8])
{ // edgesUses[8] == k axes edge
auto writeIndex = edgeIds[8];
pos[1] = startPos[2] + offset;
auto s1 = field.Get(pos[1]);
auto t = (this->IsoValue - s0) / (s1 - s0);
interpolatedEdgeIds.Set(writeIndex, pos);
weights.Set(writeIndex, static_cast<vtkm::FloatDefault>(t));
auto coord = this->InterpolateCoordinate(t, ijk, ijk + vtkm::Id3{ 0, 0, 1 });
points.Set(writeIndex, coord);
}
}
// On the boundary cells special work has to be done to cover the partial
// cell axes. These are boundary situations where the voxel axes is not
// fully formed. These situations occur on the +x,+y,+z volume
// boundaries. (The other cases fall through the default: case which is
// expected.)
//
// Note that loc is one of 27 regions in the volume, with (0,1,2)
// indicating (interior, min, max) along coordinate axes.
switch (loc)
{
case 2:
case 6:
case 18:
case 22: //+x
this->InterpolateEdge(
ijk, pos[0], incs, 5, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
this->InterpolateEdge(
ijk, pos[0], incs, 9, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
break;
case 8:
case 9:
case 24:
case 25: //+y
this->InterpolateEdge(
ijk, pos[0], incs, 1, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
this->InterpolateEdge(
ijk, pos[0], incs, 10, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
break;
case 32:
case 33:
case 36:
case 37: //+z
this->InterpolateEdge(
ijk, pos[0], incs, 2, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
this->InterpolateEdge(
ijk, pos[0], incs, 6, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
break;
case 10:
case 26: //+x +y
this->InterpolateEdge(
ijk, pos[0], incs, 1, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
this->InterpolateEdge(
ijk, pos[0], incs, 5, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
this->InterpolateEdge(
ijk, pos[0], incs, 9, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
this->InterpolateEdge(
ijk, pos[0], incs, 10, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
this->InterpolateEdge(
ijk, pos[0], incs, 11, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
break;
case 34:
case 38: //+x +z
this->InterpolateEdge(
ijk, pos[0], incs, 2, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
this->InterpolateEdge(
ijk, pos[0], incs, 5, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
this->InterpolateEdge(
ijk, pos[0], incs, 9, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
this->InterpolateEdge(
ijk, pos[0], incs, 6, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
this->InterpolateEdge(
ijk, pos[0], incs, 7, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
break;
case 40:
case 41: //+y +z
this->InterpolateEdge(
ijk, pos[0], incs, 1, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
this->InterpolateEdge(
ijk, pos[0], incs, 2, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
this->InterpolateEdge(
ijk, pos[0], incs, 3, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
this->InterpolateEdge(
ijk, pos[0], incs, 6, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
this->InterpolateEdge(
ijk, pos[0], incs, 10, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
break;
case 42: //+x +y +z happens no more than once per volume
this->InterpolateEdge(
ijk, pos[0], incs, 1, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
this->InterpolateEdge(
ijk, pos[0], incs, 2, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
this->InterpolateEdge(
ijk, pos[0], incs, 3, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
this->InterpolateEdge(
ijk, pos[0], incs, 5, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
this->InterpolateEdge(
ijk, pos[0], incs, 9, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
this->InterpolateEdge(
ijk, pos[0], incs, 10, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
this->InterpolateEdge(
ijk, pos[0], incs, 11, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
this->InterpolateEdge(
ijk, pos[0], incs, 6, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
this->InterpolateEdge(
ijk, pos[0], incs, 7, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
break;
default: // interior, or -x,-y,-z boundaries
return;
}
}
// Indicate whether voxel axes need processing for this case.
//----------------------------------------------------------------------------
template <typename WholeField,
typename WholeIEdgeField,
typename WholeWeightField,
typename WholePointField>
VTKM_EXEC inline void InterpolateEdge(const vtkm::Id3& ijk,
vtkm::Id currentIdx,
const vtkm::Id3& incs,
vtkm::Id edgeNum,
vtkm::UInt8 const* const edgeUses,
vtkm::Id* edgeIds,
const WholeField& field,
const WholeIEdgeField& interpolatedEdgeIds,
const WholeWeightField& weights,
const WholePointField& points) const
{
// if this edge is not used then get out
if (!edgeUses[edgeNum])
{
return;
}
const vtkm::Id writeIndex = edgeIds[edgeNum];
// build the edge information
vtkm::Vec<vtkm::UInt8, 2> verts = data::GetVertMap(edgeNum);
vtkm::Id3 offsets1 = data::GetVertOffsets(AxisToSum{}, verts[0]);
vtkm::Id3 offsets2 = data::GetVertOffsets(AxisToSum{}, verts[1]);
vtkm::Id2 iEdge(currentIdx + vtkm::Dot(offsets1, incs), currentIdx + vtkm::Dot(offsets2, incs));
interpolatedEdgeIds.Set(writeIndex, iEdge);
auto s0 = field.Get(iEdge[0]);
auto s1 = field.Get(iEdge[1]);
auto t = (this->IsoValue - s0) / (s1 - s0);
weights.Set(writeIndex, static_cast<vtkm::FloatDefault>(t));
auto coord = this->InterpolateCoordinate(t, ijk + offsets1, ijk + offsets2);
points.Set(writeIndex, coord);
}
//----------------------------------------------------------------------------
inline VTKM_EXEC vtkm::Vec3f InterpolateCoordinate(T t,
const vtkm::Id3& ijk0,
const vtkm::Id3& ijk1) const
{
return vtkm::Vec3f(
this->Origin[0] +
this->Spacing[0] * static_cast<vtkm::FloatDefault>(ijk0[0] + t * (ijk1[0] - ijk0[0])),
this->Origin[1] +
this->Spacing[1] * static_cast<vtkm::FloatDefault>(ijk0[1] + t * (ijk1[1] - ijk0[1])),
this->Origin[2] +
this->Spacing[2] * static_cast<vtkm::FloatDefault>(ijk0[2] + t * (ijk1[2] - ijk0[2])));
return true;
}
};
}

190
vtkm/worklet/contour/FlyingEdgesPass4Common.h Normal file
View File

@ -0,0 +1,190 @@
//============================================================================
// 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.
//============================================================================
#ifndef vtk_m_worklet_contour_flyingedges_pass4_common_h
#define vtk_m_worklet_contour_flyingedges_pass4_common_h
#include <vtkm/worklet/contour/FlyingEdgesHelpers.h>
#include <vtkm/worklet/contour/FlyingEdgesTables.h>
namespace vtkm
{
namespace worklet
{
namespace flying_edges
{
VTKM_EXEC inline vtkm::Id3 compute_incs3d(const vtkm::Id3& dims)
{
return vtkm::Id3{ 1, dims[0], (dims[0] * dims[1]) };
}
VTKM_EXEC inline bool case_includes_axes(vtkm::UInt8 const* const edgeUses)
{
return (edgeUses[0] != 0 || edgeUses[4] != 0 || edgeUses[8] != 0);
}
template <typename WholeConnField, typename WholeCellIdField>
VTKM_EXEC inline void generate_tris(vtkm::Id inputCellId,
vtkm::UInt8 edgeCase,
vtkm::UInt8 numTris,
vtkm::Id* edgeIds,
vtkm::Id& triId,
const WholeConnField& conn,
const WholeCellIdField& cellIds)
{
auto* edges = data::GetTriEdgeCases(edgeCase);
vtkm::Id edgeIndex = 1;
vtkm::Id index = static_cast<vtkm::Id>(triId) * 3;
for (vtkm::UInt8 i = 0; i < numTris; ++i)
{
cellIds.Set(triId + i, inputCellId);
//This keeps the same winding for the triangles that marching cells
//produced. By keeping the winding the same we make sure
//that 'fast' normals are consistent with the marching
//cells version
conn.Set(index, edgeIds[edges[edgeIndex]]);
conn.Set(index + 1, edgeIds[edges[edgeIndex + 2]]);
conn.Set(index + 2, edgeIds[edges[edgeIndex + 1]]);
index += 3;
edgeIndex += 3;
}
triId += numTris;
}
// Helper function to set up the point ids on voxel edges.
//----------------------------------------------------------------------------
template <typename AxisToSum, typename FieldInPointId3>
VTKM_EXEC inline void init_voxelIds(AxisToSum,
vtkm::Id writeOffset,
vtkm::UInt8 edgeCase,
const FieldInPointId3& axis_sums,
vtkm::Id* edgeIds)
{
auto* edgeUses = data::GetEdgeUses(edgeCase);
edgeIds[0] = writeOffset + axis_sums[0][AxisToSum::xindex]; // x-edges
edgeIds[1] = writeOffset + axis_sums[1][AxisToSum::xindex];
edgeIds[2] = writeOffset + axis_sums[3][AxisToSum::xindex];
edgeIds[3] = writeOffset + axis_sums[2][AxisToSum::xindex];
edgeIds[4] = writeOffset + axis_sums[0][AxisToSum::yindex]; // y-edges
edgeIds[5] = edgeIds[4] + edgeUses[4];
edgeIds[6] = writeOffset + axis_sums[3][AxisToSum::yindex];
edgeIds[7] = edgeIds[6] + edgeUses[6];
edgeIds[8] = writeOffset + axis_sums[0][AxisToSum::zindex]; // z-edges
edgeIds[9] = edgeIds[8] + edgeUses[8];
edgeIds[10] = writeOffset + axis_sums[1][AxisToSum::zindex];
edgeIds[11] = edgeIds[10] + edgeUses[10];
}
// Helper function to advance the point ids along voxel rows.
//----------------------------------------------------------------------------
VTKM_EXEC inline void advance_voxelIds(vtkm::UInt8 const* const edgeUses, vtkm::Id* edgeIds)
{
edgeIds[0] += edgeUses[0]; // x-edges
edgeIds[1] += edgeUses[1];
edgeIds[2] += edgeUses[2];
edgeIds[3] += edgeUses[3];
edgeIds[4] += edgeUses[4]; // y-edges
edgeIds[5] = edgeIds[4] + edgeUses[5];
edgeIds[6] += edgeUses[6];
edgeIds[7] = edgeIds[6] + edgeUses[7];
edgeIds[8] += edgeUses[8]; // z-edges
edgeIds[9] = edgeIds[8] + edgeUses[9];
edgeIds[10] += edgeUses[10];
edgeIds[11] = edgeIds[10] + edgeUses[11];
}
//----------------------------------------------------------------------------
struct Pass4TrimState
{
vtkm::Id left, right;
vtkm::Id3 ijk;
vtkm::Id4 startPos;
vtkm::Id axis_inc;
vtkm::UInt8 yzLoc;
bool valid = true;
template <typename AxisToSum,
typename ThreadIndices,
typename FieldInPointId,
typename WholeEdgeField>
VTKM_EXEC Pass4TrimState(AxisToSum,
const vtkm::Id3& pdims,
const ThreadIndices& threadIndices,
const FieldInPointId& axis_mins,
const FieldInPointId& axis_maxs,
const WholeEdgeField& edges)
{
// find adjusted trim values.
left = vtkm::Min(axis_mins[0], axis_mins[1]);
left = vtkm::Min(left, axis_mins[2]);
left = vtkm::Min(left, axis_mins[3]);
right = vtkm::Max(axis_maxs[0], axis_maxs[1]);
right = vtkm::Max(right, axis_maxs[2]);
right = vtkm::Max(right, axis_maxs[3]);
ijk = compute_ijk(AxisToSum{}, threadIndices.GetInputIndex3D());
startPos = compute_neighbor_starts(AxisToSum{}, ijk, pdims);
axis_inc = compute_inc(AxisToSum{}, pdims);
if (left == pdims[AxisToSum::xindex] && right == 0)
{
//verify that we have nothing to generate and early terminate.
bool mins_same = (axis_mins[0] == axis_mins[1] && axis_mins[0] == axis_mins[2] &&
axis_mins[0] == axis_mins[3]);
bool maxs_same = (axis_maxs[0] == axis_maxs[1] && axis_maxs[0] == axis_maxs[2] &&
axis_maxs[0] == axis_maxs[3]);
if (mins_same && maxs_same)
{
valid = false;
return;
}
else
{
left = 0;
right = pdims[AxisToSum::xindex] - 1;
}
}
// The trim edges may need adjustment if the contour travels between rows
// of edges (without intersecting these edges). This means checking
// whether the trim faces at (left,right) made up of the edges intersect
// the contour.
adjustTrimBounds(pdims[AxisToSum::xindex] - 1, edges, startPos, axis_inc, left, right);
if (left == right)
{
valid = false;
return;
}
const vtkm::UInt8 yLoc =
(ijk[AxisToSum::yindex] < 1
? FlyingEdges3D::MinBoundary
: (ijk[AxisToSum::yindex] >= (pdims[AxisToSum::yindex] - 2) ? FlyingEdges3D::MaxBoundary
: FlyingEdges3D::Interior));
const vtkm::UInt8 zLoc =
(ijk[AxisToSum::zindex] < 1
? FlyingEdges3D::MinBoundary
: (ijk[AxisToSum::zindex] >= (pdims[AxisToSum::zindex] - 2) ? FlyingEdges3D::MaxBoundary
: FlyingEdges3D::Interior));
yzLoc = static_cast<vtkm::UInt8>((yLoc << 2) | (zLoc << 4));
}
};
}
}
}
#endif

388
vtkm/worklet/contour/FlyingEdgesPass4X.h Normal file
View File

@ -0,0 +1,388 @@
//============================================================================
// 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.
//============================================================================
#ifndef vtk_m_worklet_contour_flyingedges_pass4x_h
#define vtk_m_worklet_contour_flyingedges_pass4x_h
#include <vtkm/worklet/contour/FlyingEdgesHelpers.h>
#include <vtkm/worklet/contour/FlyingEdgesTables.h>
namespace vtkm
{
namespace worklet
{
namespace flying_edges
{
template <typename T>
struct ComputePass4X : public vtkm::worklet::WorkletVisitCellsWithPoints
{
vtkm::Id3 PointDims;
vtkm::Vec3f Origin;
vtkm::Vec3f Spacing;
T IsoValue;
vtkm::Id CellWriteOffset;
vtkm::Id PointWriteOffset;
ComputePass4X() {}
ComputePass4X(T value,
const vtkm::Id3& pdims,
const vtkm::Vec3f& origin,
const vtkm::Vec3f& spacing,
vtkm::Id multiContourCellOffset,
vtkm::Id multiContourPointOffset)
: PointDims(pdims)
, Origin(origin)
, Spacing(spacing)
, IsoValue(value)
, CellWriteOffset(multiContourCellOffset)
, PointWriteOffset(multiContourPointOffset)
{
}
using ControlSignature = void(CellSetIn,
FieldInPoint axis_sums,
FieldInPoint axis_mins,
FieldInPoint axis_maxs,
WholeArrayIn cell_tri_count,
WholeArrayIn edgeData,
WholeArrayIn data,
WholeArrayOut connectivity,
WholeArrayOut edgeIds,
WholeArrayOut weights,
WholeArrayOut inputCellIds,
WholeArrayOut points);
using ExecutionSignature =
void(ThreadIndices, _2, _3, _4, _5, _6, _7, _8, _9, _10, _11, _12, WorkIndex);
template <typename ThreadIndices,
typename FieldInPointId3,
typename FieldInPointId,
typename WholeTriField,
typename WholeEdgeField,
typename WholeDataField,
typename WholeConnField,
typename WholeEdgeIdField,
typename WholeWeightField,
typename WholeCellIdField,
typename WholePointField>
VTKM_EXEC void operator()(const ThreadIndices& threadIndices,
const FieldInPointId3& axis_sums,
const FieldInPointId& axis_mins,
const FieldInPointId& axis_maxs,
const WholeTriField& cellTriCount,
const WholeEdgeField& edges,
const WholeDataField& field,
const WholeConnField& conn,
const WholeEdgeIdField& interpolatedEdgeIds,
const WholeWeightField& weights,
const WholeCellIdField& inputCellIds,
const WholePointField& points,
vtkm::Id oidx) const
{
using AxisToSum = SumXAxis;
//This works as cellTriCount was computed with ScanExtended
//and therefore has one more entry than the number of cells
vtkm::Id cell_tri_offset = cellTriCount.Get(oidx);
vtkm::Id next_tri_offset = cellTriCount.Get(oidx + 1);
if (cell_tri_offset == next_tri_offset)
{ //we produce nothing
return;
}
cell_tri_offset += this->CellWriteOffset;
const Pass4TrimState state(
AxisToSum{}, this->PointDims, threadIndices, axis_mins, axis_maxs, edges);
if (!state.valid)
{
return;
}
const vtkm::Id3 pdims = this->PointDims;
const vtkm::Id3 increments = compute_incs3d(pdims);
vtkm::Id edgeIds[12];
auto edgeCase = getEdgeCase(edges, state.startPos, (state.axis_inc * state.left));
init_voxelIds(AxisToSum{}, this->PointWriteOffset, edgeCase, axis_sums, edgeIds);
for (vtkm::Id i = state.left; i < state.right; ++i) // run along the trimmed voxels
{
auto ijk = state.ijk;
ijk[AxisToSum::xindex] = i;
edgeCase = getEdgeCase(edges, state.startPos, (state.axis_inc * i));
vtkm::UInt8 numTris = data::GetNumberOfPrimitives(edgeCase);
if (numTris > 0)
{
//compute what the current cellId is
vtkm::Id cellId = compute_start(AxisToSum{}, ijk, pdims - vtkm::Id3{ 1, 1, 1 });
// Start by generating triangles for this case
generate_tris(cellId, edgeCase, numTris, edgeIds, cell_tri_offset, conn, inputCellIds);
// Now generate edgeIds and weights along voxel axes if needed. Remember to take
// boundary into account.
vtkm::UInt8 loc = static_cast<vtkm::UInt8>(
state.yzLoc | (i < 1 ? FlyingEdges3D::MinBoundary
: (i >= (pdims[AxisToSum::xindex] - 2) ? FlyingEdges3D::MaxBoundary
: FlyingEdges3D::Interior)));
auto* edgeUses = data::GetEdgeUses(edgeCase);
if (loc != FlyingEdges3D::Interior || case_includes_axes(edgeUses))
{
this->Generate(loc,
ijk,
field,
interpolatedEdgeIds,
weights,
points,
state.startPos,
increments,
(state.axis_inc * i),
edgeUses,
edgeIds);
}
advance_voxelIds(edgeUses, edgeIds);
}
}
}
//----------------------------------------------------------------------------
template <typename WholeDataField,
typename WholeIEdgeField,
typename WholeWeightField,
typename WholePointField>
VTKM_EXEC inline void Generate(vtkm::UInt8 loc,
const vtkm::Id3& ijk,
const WholeDataField& field,
const WholeIEdgeField& interpolatedEdgeIds,
const WholeWeightField& weights,
const WholePointField& points,
const vtkm::Id4& startPos,
const vtkm::Id3& incs,
vtkm::Id offset,
vtkm::UInt8 const* const edgeUses,
vtkm::Id* edgeIds) const
{
using AxisToSum = SumXAxis;
vtkm::Id2 pos(startPos[0] + offset, 0);
{
auto s0 = field.Get(pos[0]);
//EdgesUses 0,4,8 work for Y axis
if (edgeUses[0])
{ // edgesUses[0] == i axes edge
auto writeIndex = edgeIds[0];
pos[1] = startPos[0] + offset + incs[AxisToSum::xindex];
auto s1 = field.Get(pos[1]);
T t = static_cast<T>((this->IsoValue - s0) / (s1 - s0));
interpolatedEdgeIds.Set(writeIndex, pos);
weights.Set(writeIndex, static_cast<vtkm::FloatDefault>(t));
auto coord = this->InterpolateCoordinate(t, ijk, ijk + vtkm::Id3{ 1, 0, 0 });
points.Set(writeIndex, coord);
}
if (edgeUses[4])
{ // edgesUses[4] == j axes edge
auto writeIndex = edgeIds[4];
pos[1] = startPos[1] + offset;
auto s1 = field.Get(pos[1]);
T t = static_cast<T>((this->IsoValue - s0) / (s1 - s0));
interpolatedEdgeIds.Set(writeIndex, pos);
weights.Set(writeIndex, static_cast<vtkm::FloatDefault>(t));
auto coord = this->InterpolateCoordinate(t, ijk, ijk + vtkm::Id3{ 0, 1, 0 });
points.Set(writeIndex, coord);
}
if (edgeUses[8])
{ // edgesUses[8] == k axes edge
auto writeIndex = edgeIds[8];
pos[1] = startPos[2] + offset;
auto s1 = field.Get(pos[1]);
T t = static_cast<T>((this->IsoValue - s0) / (s1 - s0));
interpolatedEdgeIds.Set(writeIndex, pos);
weights.Set(writeIndex, static_cast<vtkm::FloatDefault>(t));
auto coord = this->InterpolateCoordinate(t, ijk, ijk + vtkm::Id3{ 0, 0, 1 });
points.Set(writeIndex, coord);
}
}
// On the boundary cells special work has to be done to cover the partial
// cell axes. These are boundary situations where the voxel axes is not
// fully formed. These situations occur on the +x,+y,+z volume
// boundaries. (The other cases fall through the default: case which is
// expected.)
//
// Note that loc is one of 27 regions in the volume, with (0,1,2)
// indicating (interior, min, max) along coordinate axes.
switch (loc)
{
case 2:
case 6:
case 18:
case 22: //+x
this->InterpolateEdge(
ijk, pos[0], incs, 5, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
this->InterpolateEdge(
ijk, pos[0], incs, 9, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
break;
case 8:
case 9:
case 24:
case 25: //+y
this->InterpolateEdge(
ijk, pos[0], incs, 1, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
this->InterpolateEdge(
ijk, pos[0], incs, 10, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
break;
case 32:
case 33:
case 36:
case 37: //+z
this->InterpolateEdge(
ijk, pos[0], incs, 2, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
this->InterpolateEdge(
ijk, pos[0], incs, 6, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
break;
case 10:
case 26: //+x +y
this->InterpolateEdge(
ijk, pos[0], incs, 1, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
this->InterpolateEdge(
ijk, pos[0], incs, 5, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
this->InterpolateEdge(
ijk, pos[0], incs, 9, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
this->InterpolateEdge(
ijk, pos[0], incs, 10, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
this->InterpolateEdge(
ijk, pos[0], incs, 11, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
break;
case 34:
case 38: //+x +z
this->InterpolateEdge(
ijk, pos[0], incs, 2, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
this->InterpolateEdge(
ijk, pos[0], incs, 5, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
this->InterpolateEdge(
ijk, pos[0], incs, 9, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
this->InterpolateEdge(
ijk, pos[0], incs, 6, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
this->InterpolateEdge(
ijk, pos[0], incs, 7, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
break;
case 40:
case 41: //+y +z
this->InterpolateEdge(
ijk, pos[0], incs, 1, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
this->InterpolateEdge(
ijk, pos[0], incs, 2, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
this->InterpolateEdge(
ijk, pos[0], incs, 3, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
this->InterpolateEdge(
ijk, pos[0], incs, 6, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
this->InterpolateEdge(
ijk, pos[0], incs, 10, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
break;
case 42: //+x +y +z happens no more than once per volume
this->InterpolateEdge(
ijk, pos[0], incs, 1, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
this->InterpolateEdge(
ijk, pos[0], incs, 2, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
this->InterpolateEdge(
ijk, pos[0], incs, 3, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
this->InterpolateEdge(
ijk, pos[0], incs, 5, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
this->InterpolateEdge(
ijk, pos[0], incs, 9, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
this->InterpolateEdge(
ijk, pos[0], incs, 10, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
this->InterpolateEdge(
ijk, pos[0], incs, 11, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
this->InterpolateEdge(
ijk, pos[0], incs, 6, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
this->InterpolateEdge(
ijk, pos[0], incs, 7, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points);
break;
default: // interior, or -x,-y,-z boundaries
return;
}
}
// Indicate whether voxel axes need processing for this case.
//----------------------------------------------------------------------------
template <typename WholeField,
typename WholeIEdgeField,
typename WholeWeightField,
typename WholePointField>
VTKM_EXEC inline void InterpolateEdge(const vtkm::Id3& ijk,
vtkm::Id currentIdx,
const vtkm::Id3& incs,
vtkm::Id edgeNum,
vtkm::UInt8 const* const edgeUses,
vtkm::Id* edgeIds,
const WholeField& field,
const WholeIEdgeField& interpolatedEdgeIds,
const WholeWeightField& weights,
const WholePointField& points) const
{
using AxisToSum = SumXAxis;
// if this edge is not used then get out
if (!edgeUses[edgeNum])
{
return;
}
const vtkm::Id writeIndex = edgeIds[edgeNum];
// build the edge information
vtkm::Vec<vtkm::UInt8, 2> verts = data::GetVertMap(edgeNum);
vtkm::Id3 offsets1 = data::GetVertOffsets(AxisToSum{}, verts[0]);
vtkm::Id3 offsets2 = data::GetVertOffsets(AxisToSum{}, verts[1]);
vtkm::Id2 iEdge(currentIdx + vtkm::Dot(offsets1, incs), currentIdx + vtkm::Dot(offsets2, incs));
interpolatedEdgeIds.Set(writeIndex, iEdge);
auto s0 = field.Get(iEdge[0]);
auto s1 = field.Get(iEdge[1]);
T t = static_cast<T>((this->IsoValue - s0) / (s1 - s0));
weights.Set(writeIndex, static_cast<vtkm::FloatDefault>(t));
auto coord = this->InterpolateCoordinate(t, ijk + offsets1, ijk + offsets2);
points.Set(writeIndex, coord);
}
//----------------------------------------------------------------------------
inline VTKM_EXEC vtkm::Vec3f InterpolateCoordinate(T t,
const vtkm::Id3& ijk0,
const vtkm::Id3& ijk1) const
{
return vtkm::Vec3f(
this->Origin[0] +
this->Spacing[0] * static_cast<vtkm::FloatDefault>(ijk0[0] + t * (ijk1[0] - ijk0[0])),
this->Origin[1] +
this->Spacing[1] * static_cast<vtkm::FloatDefault>(ijk0[1] + t * (ijk1[1] - ijk0[1])),
this->Origin[2] +
this->Spacing[2] * static_cast<vtkm::FloatDefault>(ijk0[2] + t * (ijk1[2] - ijk0[2])));
}
};
}
}
}
#endif

36
vtkm/worklet/contour/FlyingEdgesPass4WithNormals.h → vtkm/worklet/contour/FlyingEdgesPass4XWithNormals.h
View File

@ -10,8 +10,8 @@
//============================================================================
#ifndef vtk_m_worklet_contour_flyingedges_pass4_with_norms_h
#define vtk_m_worklet_contour_flyingedges_pass4_with_norms_h
#ifndef vtk_m_worklet_contour_flyingedges_pass4x_with_norms_h
#define vtk_m_worklet_contour_flyingedges_pass4x_with_norms_h
#include <vtkm/worklet/contour/FlyingEdgesHelpers.h>
@ -26,8 +26,8 @@ namespace worklet
namespace flying_edges
{
template <typename T, typename AxisToSum>
struct ComputePass4WithNormals : public vtkm::worklet::WorkletVisitCellsWithPoints
template <typename T>
struct ComputePass4XWithNormals : public vtkm::worklet::WorkletVisitCellsWithPoints
{
vtkm::Id3 PointDims;
@ -39,13 +39,13 @@ struct ComputePass4WithNormals : public vtkm::worklet::WorkletVisitCellsWithPoin
vtkm::Id CellWriteOffset;
vtkm::Id PointWriteOffset;
ComputePass4WithNormals() {}
ComputePass4WithNormals(T value,
const vtkm::Id3& pdims,
const vtkm::Vec3f& origin,
const vtkm::Vec3f& spacing,
vtkm::Id multiContourCellOffset,
vtkm::Id multiContourPointOffset)
ComputePass4XWithNormals() {}
ComputePass4XWithNormals(T value,
const vtkm::Id3& pdims,
const vtkm::Vec3f& origin,
const vtkm::Vec3f& spacing,
vtkm::Id multiContourCellOffset,
vtkm::Id multiContourPointOffset)
: PointDims(pdims)
, Origin(origin)
, Spacing(spacing)
@ -98,6 +98,8 @@ struct ComputePass4WithNormals : public vtkm::worklet::WorkletVisitCellsWithPoin
const WholeNormalsField& normals,
vtkm::Id oidx) const
{
using AxisToSum = SumXAxis;
//This works as cellTriCount was computed with ScanExtended
//and therefore has one more entry than the number of cells
vtkm::Id cell_tri_offset = cellTriCount.Get(oidx);
@ -181,6 +183,8 @@ struct ComputePass4WithNormals : public vtkm::worklet::WorkletVisitCellsWithPoin
vtkm::UInt8 const* const edgeUses,
vtkm::Id* edgeIds) const
{
using AxisToSum = SumXAxis;
vtkm::Id2 pos(startPos[0] + offset, 0);
{
auto s0 = field.Get(pos[0]);
@ -192,7 +196,7 @@ struct ComputePass4WithNormals : public vtkm::worklet::WorkletVisitCellsWithPoin
auto writeIndex = edgeIds[0];
pos[1] = startPos[0] + offset + incs[AxisToSum::xindex];
auto s1 = field.Get(pos[1]);
auto t = (this->IsoValue - s0) / (s1 - s0);
T t = static_cast<T>((this->IsoValue - s0) / (s1 - s0));
interpolatedEdgeIds.Set(writeIndex, pos);
weights.Set(writeIndex, static_cast<vtkm::FloatDefault>(t));
@ -211,7 +215,7 @@ struct ComputePass4WithNormals : public vtkm::worklet::WorkletVisitCellsWithPoin
auto writeIndex = edgeIds[4];
pos[1] = startPos[1] + offset;
auto s1 = field.Get(pos[1]);
auto t = (this->IsoValue - s0) / (s1 - s0);
T t = static_cast<T>((this->IsoValue - s0) / (s1 - s0));
interpolatedEdgeIds.Set(writeIndex, pos);
weights.Set(writeIndex, static_cast<vtkm::FloatDefault>(t));
@ -230,7 +234,7 @@ struct ComputePass4WithNormals : public vtkm::worklet::WorkletVisitCellsWithPoin
auto writeIndex = edgeIds[8];
pos[1] = startPos[2] + offset;
auto s1 = field.Get(pos[1]);
auto t = (this->IsoValue - s0) / (s1 - s0);
T t = static_cast<T>((this->IsoValue - s0) / (s1 - s0));
interpolatedEdgeIds.Set(writeIndex, pos);
weights.Set(writeIndex, static_cast<vtkm::FloatDefault>(t));
@ -369,6 +373,8 @@ struct ComputePass4WithNormals : public vtkm::worklet::WorkletVisitCellsWithPoin
const WholePointField& points,
const WholeNormalField& normals) const
{
using AxisToSum = SumXAxis;
// if this edge is not used then get out
if (!edgeUses[edgeNum])
{
@ -388,7 +394,7 @@ struct ComputePass4WithNormals : public vtkm::worklet::WorkletVisitCellsWithPoin
auto s0 = field.Get(iEdge[0]);
auto s1 = field.Get(iEdge[1]);
auto t = (this->IsoValue - s0) / (s1 - s0);
T t = static_cast<T>((this->IsoValue - s0) / (s1 - s0));
weights.Set(writeIndex, static_cast<vtkm::FloatDefault>(t));
auto coord = this->InterpolateCoordinate(t, ijk + offsets1, ijk + offsets2);

425
vtkm/worklet/contour/FlyingEdgesPass4Y.h Normal file
View File

@ -0,0 +1,425 @@
//============================================================================
// 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.
//============================================================================
#ifndef vtk_m_worklet_contour_flyingedges_pass4y_h
#define vtk_m_worklet_contour_flyingedges_pass4y_h
#include <vtkm/worklet/contour/FlyingEdgesHelpers.h>
#include <vtkm/worklet/contour/FlyingEdgesTables.h>
#include <vtkm/VectorAnalysis.h>
#include <vtkm/worklet/gradient/StructuredPointGradient.h>
namespace vtkm
{
namespace worklet
{
namespace flying_edges
{
template <typename T>
struct ComputePass4Y : public vtkm::worklet::WorkletVisitCellsWithPoints
{
vtkm::Id3 PointDims;
T IsoValue;
vtkm::Id CellWriteOffset;
vtkm::Id PointWriteOffset;
ComputePass4Y() {}
ComputePass4Y(T value,
const vtkm::Id3& pdims,
vtkm::Id multiContourCellOffset,
vtkm::Id multiContourPointOffset)
: PointDims(pdims)
, IsoValue(value)
, CellWriteOffset(multiContourCellOffset)
, PointWriteOffset(multiContourPointOffset)
{
}
using ControlSignature = void(CellSetIn,
FieldInPoint axis_sums,
FieldInPoint axis_mins,
FieldInPoint axis_maxs,
WholeArrayIn cell_tri_count,
WholeArrayIn edgeData,
WholeArrayIn data,
WholeArrayOut connectivity,
WholeArrayOut edgeIds,
WholeArrayOut weights,
WholeArrayOut inputCellIds);
using ExecutionSignature =
void(ThreadIndices, _2, _3, _4, _5, _6, _7, _8, _9, _10, _11, WorkIndex);
template <typename ThreadIndices,
typename FieldInPointId3,
typename FieldInPointId,
typename WholeTriField,
typename WholeEdgeField,
typename WholeDataField,
typename WholeConnField,
typename WholeEdgeIdField,
typename WholeWeightField,
typename WholeCellIdField>
VTKM_EXEC void operator()(const ThreadIndices& threadIndices,
const FieldInPointId3& axis_sums,
const FieldInPointId& axis_mins,
const FieldInPointId& axis_maxs,
const WholeTriField& cellTriCount,
const WholeEdgeField& edges,
const WholeDataField& field,
const WholeConnField& conn,
const WholeEdgeIdField& interpolatedEdgeIds,
const WholeWeightField& weights,
const WholeCellIdField& inputCellIds,
vtkm::Id oidx) const
{
using AxisToSum = SumYAxis;
//This works as cellTriCount was computed with ScanExtended
//and therefore has one more entry than the number of cells
vtkm::Id cell_tri_offset = cellTriCount.Get(oidx);
vtkm::Id next_tri_offset = cellTriCount.Get(oidx + 1);
if (cell_tri_offset == next_tri_offset)
{ //we produce nothing
return;
}
cell_tri_offset += this->CellWriteOffset;
const Pass4TrimState state(
AxisToSum{}, this->PointDims, threadIndices, axis_mins, axis_maxs, edges);
if (!state.valid)
{
return;
}
const vtkm::Id3 pdims = this->PointDims;
const vtkm::Id3 increments = compute_incs3d(pdims);
vtkm::Id edgeIds[12];
auto edgeCase = getEdgeCase(edges, state.startPos, (state.axis_inc * state.left));
init_voxelIds(AxisToSum{}, this->PointWriteOffset, edgeCase, axis_sums, edgeIds);
for (vtkm::Id i = state.left; i < state.right; ++i) // run along the trimmed voxels
{
auto ijk = state.ijk;
ijk[AxisToSum::xindex] = i;
edgeCase = getEdgeCase(edges, state.startPos, (state.axis_inc * i));
vtkm::UInt8 numTris = data::GetNumberOfPrimitives(edgeCase);
if (numTris > 0)
{
//compute what the current cellId is
vtkm::Id cellId = compute_start(AxisToSum{}, ijk, pdims - vtkm::Id3{ 1, 1, 1 });
// Start by generating triangles for this case
generate_tris(cellId, edgeCase, numTris, edgeIds, cell_tri_offset, conn, inputCellIds);
// Now generate edgeIds and weights along voxel axes if needed. Remember to take
// boundary into account.
vtkm::UInt8 loc = static_cast<vtkm::UInt8>(
state.yzLoc | (i < 1 ? FlyingEdges3D::MinBoundary
: (i >= (pdims[AxisToSum::xindex] - 2) ? FlyingEdges3D::MaxBoundary
: FlyingEdges3D::Interior)));
auto* edgeUses = data::GetEdgeUses(edgeCase);
if (loc != FlyingEdges3D::Interior || case_includes_axes(edgeUses))
{
this->Generate(loc,
field,
interpolatedEdgeIds,
weights,
state.startPos,
increments,
(state.axis_inc * i),
edgeUses,
edgeIds);
}
advance_voxelIds(edgeUses, edgeIds);
}
}
}
//----------------------------------------------------------------------------
template <typename WholeDataField, typename WholeIEdgeField, typename WholeWeightField>
VTKM_EXEC inline void Generate(vtkm::UInt8 loc,
const WholeDataField& field,
const WholeIEdgeField& interpolatedEdgeIds,
const WholeWeightField& weights,
const vtkm::Id4& startPos,
const vtkm::Id3& incs,
vtkm::Id offset,
vtkm::UInt8 const* const edgeUses,
vtkm::Id* edgeIds) const
{
using AxisToSum = SumYAxis;
vtkm::Id2 pos(startPos[0] + offset, 0);
{
auto s0 = field.Get(pos[0]);
//EdgesUses 0,4,8 work for Y axis
if (edgeUses[0])
{ // edgesUses[0] == i axes edge
auto writeIndex = edgeIds[0];
pos[1] = startPos[0] + offset + incs[AxisToSum::xindex];
auto s1 = field.Get(pos[1]);
T t = static_cast<T>((this->IsoValue - s0) / (s1 - s0));
interpolatedEdgeIds.Set(writeIndex, pos);
weights.Set(writeIndex, static_cast<vtkm::FloatDefault>(t));
}
if (edgeUses[4])
{ // edgesUses[4] == j axes edge
auto writeIndex = edgeIds[4];
pos[1] = startPos[1] + offset;
auto s1 = field.Get(pos[1]);
T t = static_cast<T>((this->IsoValue - s0) / (s1 - s0));
interpolatedEdgeIds.Set(writeIndex, pos);
weights.Set(writeIndex, static_cast<vtkm::FloatDefault>(t));
}
if (edgeUses[8])
{ // edgesUses[8] == k axes edge
auto writeIndex = edgeIds[8];
pos[1] = startPos[2] + offset;
auto s1 = field.Get(pos[1]);
T t = static_cast<T>((this->IsoValue - s0) / (s1 - s0));
interpolatedEdgeIds.Set(writeIndex, pos);
weights.Set(writeIndex, static_cast<vtkm::FloatDefault>(t));
}
}
// On the boundary cells special work has to be done to cover the partial
// cell axes. These are boundary situations where the voxel axes is not
// fully formed. These situations occur on the +x,+y,+z volume
// boundaries. (The other cases fall through the default: case which is
// expected.)
//
// Note that loc is one of 27 regions in the volume, with (0,1,2)
// indicating (interior, min, max) along coordinate axes.
switch (loc)
{
case 2:
case 6:
case 18:
case 22: //+x
this->InterpolateEdge(
pos[0], incs, 5, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
this->InterpolateEdge(
pos[0], incs, 9, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
break;
case 8:
case 9:
case 24:
case 25: //+y
this->InterpolateEdge(
pos[0], incs, 1, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
this->InterpolateEdge(
pos[0], incs, 10, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
break;
case 32:
case 33:
case 36:
case 37: //+z
this->InterpolateEdge(
pos[0], incs, 2, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
this->InterpolateEdge(
pos[0], incs, 6, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
break;
case 10:
case 26: //+x +y
this->InterpolateEdge(
pos[0], incs, 1, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
this->InterpolateEdge(
pos[0], incs, 5, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
this->InterpolateEdge(
pos[0], incs, 9, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
this->InterpolateEdge(
pos[0], incs, 10, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
this->InterpolateEdge(
pos[0], incs, 11, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
break;
case 34:
case 38: //+x +z
this->InterpolateEdge(
pos[0], incs, 2, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
this->InterpolateEdge(
pos[0], incs, 5, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
this->InterpolateEdge(
pos[0], incs, 9, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
this->InterpolateEdge(
pos[0], incs, 6, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
this->InterpolateEdge(
pos[0], incs, 7, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
break;
case 40:
case 41: //+y +z
this->InterpolateEdge(
pos[0], incs, 1, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
this->InterpolateEdge(
pos[0], incs, 2, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
this->InterpolateEdge(
pos[0], incs, 3, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
this->InterpolateEdge(
pos[0], incs, 6, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
this->InterpolateEdge(
pos[0], incs, 10, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
break;
case 42: //+x +y +z happens no more than once per volume
this->InterpolateEdge(
pos[0], incs, 1, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
this->InterpolateEdge(
pos[0], incs, 2, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
this->InterpolateEdge(
pos[0], incs, 3, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
this->InterpolateEdge(
pos[0], incs, 5, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
this->InterpolateEdge(
pos[0], incs, 9, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
this->InterpolateEdge(
pos[0], incs, 10, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
this->InterpolateEdge(
pos[0], incs, 11, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
this->InterpolateEdge(
pos[0], incs, 6, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
this->InterpolateEdge(
pos[0], incs, 7, edgeUses, edgeIds, field, interpolatedEdgeIds, weights);
break;
default: // interior, or -x,-y,-z boundaries
return;
}
}
// Indicate whether voxel axes need processing for this case.
//----------------------------------------------------------------------------
template <typename WholeField, typename WholeIEdgeField, typename WholeWeightField>
VTKM_EXEC inline void InterpolateEdge(vtkm::Id currentIdx,
const vtkm::Id3& incs,
vtkm::Id edgeNum,
vtkm::UInt8 const* const edgeUses,
vtkm::Id* edgeIds,
const WholeField& field,
const WholeIEdgeField& interpolatedEdgeIds,
const WholeWeightField& weights) const
{
using AxisToSum = SumYAxis;
// if this edge is not used then get out
if (!edgeUses[edgeNum])
{
return;
}
const vtkm::Id writeIndex = edgeIds[edgeNum];
// build the edge information
vtkm::Vec<vtkm::UInt8, 2> verts = data::GetVertMap(edgeNum);
vtkm::Id3 offsets1 = data::GetVertOffsets(AxisToSum{}, verts[0]);
vtkm::Id3 offsets2 = data::GetVertOffsets(AxisToSum{}, verts[1]);
vtkm::Id2 iEdge(currentIdx + vtkm::Dot(offsets1, incs), currentIdx + vtkm::Dot(offsets2, incs));
interpolatedEdgeIds.Set(writeIndex, iEdge);
auto s0 = field.Get(iEdge[0]);
auto s1 = field.Get(iEdge[1]);
T t = static_cast<T>((this->IsoValue - s0) / (s1 - s0));
weights.Set(writeIndex, static_cast<vtkm::FloatDefault>(t));
}
};
template <typename T>
struct ComputePass5Y : public vtkm::worklet::WorkletMapField
{
vtkm::internal::ArrayPortalUniformPointCoordinates Coordinates;
vtkm::Id NormalWriteOffset;
ComputePass5Y() {}
ComputePass5Y(const vtkm::Id3& pdims,
const vtkm::Vec3f& origin,
const vtkm::Vec3f& spacing,
vtkm::Id normalWriteOffset,
bool generateNormals)
: Coordinates(pdims, origin, spacing)
, NormalWriteOffset(normalWriteOffset)
{
if (!generateNormals)
{
this->NormalWriteOffset = -1;
}
}
using ControlSignature = void(FieldIn interpEdgeIds,
FieldIn interpWeight,
FieldOut points,
WholeArrayIn field,
WholeArrayOut normals);
using ExecutionSignature = void(_1, _2, _3, _4, _5, WorkIndex);
template <typename PT, typename WholeInputField, typename WholeNormalField>
VTKM_EXEC void operator()(const vtkm::Id2& interpEdgeIds,
vtkm::FloatDefault weight,
vtkm::Vec<PT, 3>& outPoint,
const WholeInputField& field,
WholeNormalField& normals,
vtkm::Id oidx) const
{
{
vtkm::Vec3f point1 = this->Coordinates.Get(interpEdgeIds[0]);
vtkm::Vec3f point2 = this->Coordinates.Get(interpEdgeIds[1]);
outPoint = vtkm::Lerp(point1, point2, weight);
}
//NormalWriteOffset of -1 means no normals
if (this->NormalWriteOffset >= 0)
{
vtkm::Vec<T, 3> g0, g1;
const vtkm::Id3& dims = this->Coordinates.GetDimensions();
vtkm::Id3 ijk{ interpEdgeIds[0] % dims[0],
(interpEdgeIds[0] / dims[0]) % dims[1],
interpEdgeIds[0] / (dims[0] * dims[1]) };
vtkm::worklet::gradient::StructuredPointGradient gradient;
vtkm::exec::BoundaryState boundary(ijk, dims);
vtkm::exec::FieldNeighborhood<vtkm::internal::ArrayPortalUniformPointCoordinates>
coord_neighborhood(this->Coordinates, boundary);
vtkm::exec::FieldNeighborhood<WholeInputField> field_neighborhood(field, boundary);
//compute the gradient at point 1
gradient(boundary, coord_neighborhood, field_neighborhood, g0);
//compute the gradient at point 2. This optimization can be optimized
boundary.IJK = vtkm::Id3{ interpEdgeIds[1] % dims[0],
(interpEdgeIds[1] / dims[0]) % dims[1],
interpEdgeIds[1] / (dims[0] * dims[1]) };
gradient(boundary, coord_neighborhood, field_neighborhood, g1);
vtkm::Vec3f n = vtkm::Lerp(g0, g1, weight);
const auto mag2 = vtkm::MagnitudeSquared(n);
if (mag2 > 0.)
{
n = n * vtkm::RSqrt(mag2);
}
normals.Set(this->NormalWriteOffset + oidx, n);
}
}
};
}
}
}
#endif

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

@ -33,6 +33,7 @@
#include <vtkm/cont/CoordinateSystem.h>
#include <vtkm/cont/DataSet.h>
#include <vtkm/filter/CleanGrid.h>
#include <vtkm/filter/Contour.h>
#include <vtkm/filter/PolicyBase.h>
#include <vtkm/filter/SurfaceNormals.h>
@ -56,14 +57,16 @@ vtkm::cont::DataSet CreateDataSet(bool pointNormals, bool cellNormals)
wavelet.SetMagnitude({ 5 });
auto dataSet = wavelet.Execute();
// Cut a contour
vtkm::filter::CleanGrid toGrid;
// unstructured grid contour
vtkm::filter::Contour contour;
contour.SetActiveField("scalars", vtkm::cont::Field::Association::POINTS);
contour.SetNumberOfIsoValues(1);
contour.SetIsoValue(192);
contour.SetMergeDuplicatePoints(true);
contour.SetGenerateNormals(false);
dataSet = contour.Execute(dataSet);
dataSet = contour.Execute(toGrid.Execute(dataSet));
vtkm::filter::SurfaceNormals normals;
normals.SetGeneratePointNormals(pointNormals);