Flying Edges now works with multiple contours

This commit is contained in:
Robert Maynard 2020-04-10 16:59:47 -04:00
parent 6b46bfc54b
commit 2434254218
3 changed files with 71 additions and 33 deletions

@ -128,6 +128,8 @@ vtkm::cont::CellSetSingleType<> execute(
vtkm::cont::ArrayHandle<vtkm::Id> triangle_topology;
for (std::size_t i = 0; i < isovalues.size(); ++i)
{
auto multiContourCellOffset = sharedState.CellIdMap.GetNumberOfValues();
auto multiContourPointOffset = sharedState.InterpolationWeights.GetNumberOfValues();
ValueType isoval = isovalues[i];
//----------------------------------------------------------------------------
@ -161,16 +163,16 @@ vtkm::cont::CellSetSingleType<> execute(
detail::extend_by(triangle_topology, 3 * sumTris);
detail::extend_by(sharedState.CellIdMap, sumTris);
//By using the current length of points array as the starting value we
//handle multiple contours, by propagating the correct write offset
auto newPointSize = vtkm::cont::Algorithm::ScanExclusive(
metaDataLinearSums, metaDataLinearSums, vtkm::Sum{}, points.GetNumberOfValues());
auto 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
ComputePass4<ValueType, AxisToSum> worklet4(isoval, pdims);
ComputePass4<ValueType, AxisToSum> worklet4(
isoval, pdims, multiContourCellOffset, multiContourPointOffset);
invoke(worklet4,
metaDataMesh2D,
metaDataSums,

@ -55,7 +55,7 @@ VTKM_EXEC inline void generate_tris(vtkm::Id inputCellId,
vtkm::UInt8 edgeCase,
vtkm::UInt8 numTris,
vtkm::Id* edgeIds,
vtkm::Int32& triId,
vtkm::Id& triId,
const WholeConnField& conn,
const WholeCellIdField& cellIds)
{
@ -84,22 +84,23 @@ VTKM_EXEC inline void generate_tris(vtkm::Id inputCellId,
//----------------------------------------------------------------------------
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] = axis_sums[0][AxisToSum::xindex]; // x-edges
edgeIds[1] = axis_sums[1][AxisToSum::xindex];
edgeIds[2] = axis_sums[3][AxisToSum::xindex];
edgeIds[3] = axis_sums[2][AxisToSum::xindex];
edgeIds[4] = axis_sums[0][AxisToSum::yindex]; // y-edges
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] = axis_sums[3][AxisToSum::yindex];
edgeIds[6] = writeOffset + axis_sums[3][AxisToSum::yindex];
edgeIds[7] = edgeIds[6] + edgeUses[6];
edgeIds[8] = axis_sums[0][AxisToSum::zindex]; // z-edges
edgeIds[8] = writeOffset + axis_sums[0][AxisToSum::zindex]; // z-edges
edgeIds[9] = edgeIds[8] + edgeUses[8];
edgeIds[10] = axis_sums[1][AxisToSum::zindex];
edgeIds[10] = writeOffset + axis_sums[1][AxisToSum::zindex];
edgeIds[11] = edgeIds[10] + edgeUses[10];
}
@ -121,19 +122,24 @@ VTKM_EXEC inline void advance_voxelIds(vtkm::UInt8 const* const edgeUses, vtkm::
edgeIds[11] = edgeIds[10] + edgeUses[11];
}
template <typename T, typename AxisToSum>
struct ComputePass4 : public vtkm::worklet::WorkletVisitCellsWithPoints
{
vtkm::Id3 PointDims;
T IsoValue;
vtkm::Id CellWriteOffset;
vtkm::Id PointWriteOffset;
ComputePass4() {}
ComputePass4(T value, const vtkm::Id3& pdims)
ComputePass4(T value,
const vtkm::Id3& pdims,
vtkm::Id multiContourCellOffset,
vtkm::Id multiContourPointOffset)
: PointDims(pdims)
, IsoValue(value)
, CellWriteOffset(multiContourCellOffset)
, PointWriteOffset(multiContourPointOffset)
{
}
@ -176,12 +182,13 @@ struct ComputePass4 : public vtkm::worklet::WorkletVisitCellsWithPoints
{
//This works as cellTriCount was computed with ScanExtended
//and therefore has one more entry than the number of cells
vtkm::Int32 cell_tri_offset = cellTriCount.Get(oidx);
vtkm::Int32 next_tri_offset = cellTriCount.Get(oidx + 1);
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;
// find adjusted trim values.
vtkm::Id left = vtkm::Min(axis_mins[0], axis_mins[1]);
@ -242,7 +249,7 @@ struct ComputePass4 : public vtkm::worklet::WorkletVisitCellsWithPoints
const vtkm::Id3 increments = compute_incs3d(pdims);
vtkm::Id edgeIds[12];
auto edgeCase = getEdgeCase(edges, startPos, (axis_inc * left));
init_voxelIds(AxisToSum{}, edgeCase, axis_sums, edgeIds);
init_voxelIds(AxisToSum{}, this->PointWriteOffset, edgeCase, axis_sums, edgeIds);
for (vtkm::Id i = left; i < right; ++i) // run along the trimmed voxels
{
ijk[AxisToSum::xindex] = i;

@ -196,25 +196,26 @@ void TestContourUniformGrid()
vtkm::cont::ArrayHandle<vtkm::FloatDefault> cellFieldArray;
dataSet.GetField("cellvar").GetData().CopyTo(cellFieldArray);
vtkm::worklet::Contour isosurfaceFilter;
isosurfaceFilter.SetMergeDuplicatePoints(false);
vtkm::worklet::Contour contour;
contour.SetMergeDuplicatePoints(false);
std::vector<vtkm::Float32> contourValue{ 0.5f };
std::vector<vtkm::Float32> contourValue{ 0.5f, 0.5f };
const vtkm::Id numContours = static_cast<vtkm::Id>(contourValue.size());
vtkm::cont::ArrayHandle<vtkm::Vec3f_32> verticesArray;
vtkm::cont::ArrayHandle<vtkm::Vec3f_32> normalsArray;
vtkm::cont::ArrayHandle<vtkm::Float32> scalarsArray;
auto result = isosurfaceFilter.Run(contourValue,
cellSet,
dataSet.GetCoordinateSystem(),
pointFieldArray,
verticesArray,
normalsArray);
auto result = contour.Run(contourValue,
cellSet,
dataSet.GetCoordinateSystem(),
pointFieldArray,
verticesArray,
normalsArray);
scalarsArray = isosurfaceFilter.ProcessPointField(pointFieldArray);
scalarsArray = contour.ProcessPointField(pointFieldArray);
vtkm::cont::ArrayHandle<vtkm::FloatDefault> cellFieldArrayOut;
cellFieldArrayOut = isosurfaceFilter.ProcessCellField(cellFieldArray);
cellFieldArrayOut = contour.ProcessCellField(cellFieldArray);
std::cout << "vertices: ";
vtkm::cont::printSummary_ArrayHandle(verticesArray, std::cout);
@ -231,9 +232,37 @@ void TestContourUniformGrid()
VTKM_TEST_ASSERT(result.GetNumberOfCells() == cellFieldArrayOut.GetNumberOfValues());
VTKM_TEST_ASSERT(result.GetNumberOfCells() == 160);
VTKM_TEST_ASSERT(result.GetNumberOfCells() == (160 * numContours));
VTKM_TEST_ASSERT(verticesArray.GetNumberOfValues() == 72);
VTKM_TEST_ASSERT(verticesArray.GetNumberOfValues() == (72 * numContours));
// Verify that multiple contours of the same iso value are identical
{
auto normal_portal = normalsArray.ReadPortal();
for (vtkm::Id i = 0; i < 72; ++i)
{
for (vtkm::Id j = 1; j < numContours; ++j)
{
vtkm::Id jIndex = i + (72 * j);
VTKM_TEST_ASSERT(test_equal(normal_portal.Get(i), normal_portal.Get(jIndex)),
"multi contour failed");
}
}
auto outCellPortal =
result.GetConnectivityArray(vtkm::TopologyElementTagCell(), vtkm::TopologyElementTagPoint())
.ReadPortal();
for (vtkm::Id i = 0; i < 480; ++i)
{ //(3*160) as we are iterating triangle soup so the length is numOfCells*3
for (vtkm::Id j = 1; j < numContours; ++j)
{
vtkm::Id jIndex = i + (480 * j);
vtkm::Id expectedValue = (72 * j) + outCellPortal.Get(i);
VTKM_TEST_ASSERT(test_equal(outCellPortal.Get(jIndex), expectedValue),
"multi contour failed");
}
}
}
}
void TestContourExplicit()