From 1538fc02fbc5b9cbe41bc29bcf3b151a3a3dfb2b Mon Sep 17 00:00:00 2001 From: Louis Gombert Date: Fri, 28 Apr 2023 16:32:11 +0200 Subject: [PATCH] Implement Flying Edges for all structured CellSets In addition to using uniform coordinates, the ContourFlyingEdges filter can now process any type of coordinate system, making the filter use Flying Edges in more cases --- docs/changelog/flying-edges-structured.md | 5 ++ .../clean_grid/testing/UnitTestCleanGrid.cxx | 10 +-- vtkm/filter/contour/Contour.cxx | 7 +- vtkm/filter/contour/ContourFlyingEdges.cxx | 22 +---- .../contour/testing/UnitTestContourFilter.cxx | 68 ++++++++++++--- .../contour/worklet/ContourFlyingEdges.h | 6 +- .../contour/worklet/contour/FlyingEdges.h | 17 ++-- .../worklet/contour/FlyingEdgesPass4.h | 36 +++----- .../worklet/contour/FlyingEdgesPass4X.h | 80 +++++++++++------- .../contour/FlyingEdgesPass4XWithNormals.h | 82 ++++++++++++------- .../worklet/contour/FlyingEdgesPass4Y.h | 43 +++++----- 11 files changed, 212 insertions(+), 164 deletions(-) create mode 100644 docs/changelog/flying-edges-structured.md diff --git a/docs/changelog/flying-edges-structured.md b/docs/changelog/flying-edges-structured.md new file mode 100644 index 000000000..a961afb70 --- /dev/null +++ b/docs/changelog/flying-edges-structured.md @@ -0,0 +1,5 @@ +# Implement Flying Edges for structured cellsets with rectilinear and curvilinear coordinates + +When Flying Edges was introduced to compute contours of a 3D structured cellset, it could only process uniform coordinates. This limitation is now lifted : an alternative interpolation function can be used in the fourth pass of the algorithm in order to support rectilinear and curvilinear coordinate systems. + +Accordingly, the `Contour` filter now calls `ContourFlyingEdges` instead of `ContourMarchingCells` for these newly supported cases. diff --git a/vtkm/filter/clean_grid/testing/UnitTestCleanGrid.cxx b/vtkm/filter/clean_grid/testing/UnitTestCleanGrid.cxx index 740e6d9dc..c5eb5cc3f 100644 --- a/vtkm/filter/clean_grid/testing/UnitTestCleanGrid.cxx +++ b/vtkm/filter/clean_grid/testing/UnitTestCleanGrid.cxx @@ -13,7 +13,7 @@ #include #include #include -#include +#include namespace { @@ -72,13 +72,7 @@ void TestPointMerging() vtkm::cont::testing::MakeTestDataSet makeDataSet; vtkm::cont::DataSet baseData = makeDataSet.Make3DUniformDataSet3(vtkm::Id3(4, 4, 4)); - //Convert the baseData implicit points to explicit points, since the contour - //filter for uniform data always does point merging - vtkm::cont::ArrayHandle newcoords; - vtkm::cont::ArrayCopy(baseData.GetCoordinateSystem().GetData(), newcoords); - baseData.AddPointField(baseData.GetCoordinateSystemName(), newcoords); - - vtkm::filter::contour::Contour marchingCubes; + vtkm::filter::contour::ContourMarchingCells marchingCubes; marchingCubes.SetIsoValue(0.05); marchingCubes.SetMergeDuplicatePoints(false); marchingCubes.SetActiveField("pointvar"); diff --git a/vtkm/filter/contour/Contour.cxx b/vtkm/filter/contour/Contour.cxx index ebbd896c0..466a3657d 100644 --- a/vtkm/filter/contour/Contour.cxx +++ b/vtkm/filter/contour/Contour.cxx @@ -34,11 +34,8 @@ vtkm::cont::DataSet Contour::DoExecute(const vtkm::cont::DataSet& inDataSet) auto inCoords = inDataSet.GetCoordinateSystem(this->GetActiveCoordinateSystemIndex()).GetData(); std::unique_ptr implementation; - // For now, Flying Edges is only used for 3D Structured CellSets, - // using Uniform coordinates. - if (inCellSet.template IsType>() && - inCoords.template IsType< - vtkm::cont::ArrayHandle>()) + // Flying Edges is only used for 3D Structured CellSets + if (inCellSet.template IsType>()) { VTKM_LOG_S(vtkm::cont::LogLevel::Info, "Using flying edges"); implementation.reset(new vtkm::filter::contour::ContourFlyingEdges); diff --git a/vtkm/filter/contour/ContourFlyingEdges.cxx b/vtkm/filter/contour/ContourFlyingEdges.cxx index f457a0be8..457565e55 100644 --- a/vtkm/filter/contour/ContourFlyingEdges.cxx +++ b/vtkm/filter/contour/ContourFlyingEdges.cxx @@ -44,13 +44,10 @@ vtkm::cont::DataSet ContourFlyingEdges::DoExecute(const vtkm::cont::DataSet& inD const vtkm::cont::CoordinateSystem& inCoords = inDataSet.GetCoordinateSystem(this->GetActiveCoordinateSystemIndex()); - if (!inCellSet.template IsType>() || - !inCoords.GetData() - .template IsType< - vtkm::cont::ArrayHandle>()) + if (!inCellSet.template IsType>()) { throw vtkm::cont::ErrorFilterExecution("This filter is only available for 3-Dimensional " - "Structured Cell Sets using uniform point coordinates."); + "Structured Cell Sets"); } // Get the CellSet's known dynamic type @@ -74,22 +71,11 @@ vtkm::cont::DataSet ContourFlyingEdges::DoExecute(const vtkm::cont::DataSet& inD if (this->GenerateNormals && !this->GetComputeFastNormals()) { - outputCells = worklet.Run( - ivalues, - inputCells, - inCoords.GetData().AsArrayHandle(), - concrete, - vertices, - normals); + outputCells = worklet.Run(ivalues, inputCells, inCoords, concrete, vertices, normals); } else { - outputCells = worklet.Run( - ivalues, - inputCells, - inCoords.GetData().AsArrayHandle(), - concrete, - vertices); + outputCells = worklet.Run(ivalues, inputCells, inCoords, concrete, vertices); } }; diff --git a/vtkm/filter/contour/testing/UnitTestContourFilter.cxx b/vtkm/filter/contour/testing/UnitTestContourFilter.cxx index d1e159ce5..86aaa3aeb 100644 --- a/vtkm/filter/contour/testing/UnitTestContourFilter.cxx +++ b/vtkm/filter/contour/testing/UnitTestContourFilter.cxx @@ -176,24 +176,12 @@ public: void TestUnsupportedFlyingEdges() const { vtkm::cont::testing::MakeTestDataSet maker; - - vtkm::cont::DataSet rectilinearDataset = maker.Make3DRectilinearDataSet0(); vtkm::cont::DataSet explicitDataSet = maker.Make3DExplicitDataSet0(); vtkm::filter::contour::ContourFlyingEdges filter; filter.SetIsoValue(2.0); filter.SetActiveField("pointvar"); - try - { - filter.Execute(rectilinearDataset); - VTKM_TEST_FAIL("Flying Edges filter should not run on datasets with rectilinear coordinates"); - } - catch (vtkm::cont::ErrorFilterExecution&) - { - std::cout << "Execution successfully aborted" << std::endl; - } - try { filter.Execute(explicitDataSet); @@ -205,6 +193,58 @@ public: } } + template + void TestNonUniformStructured() const + { + auto pathname = + vtkm::cont::testing::Testing::DataPath("rectilinear/simple_rectilinear1_ascii.vtk"); + vtkm::io::VTKDataSetReader reader(pathname); + vtkm::cont::DataSet rectilinearDataset = reader.ReadDataSet(); + + // Single-cell contour + ContourFilterType filter; + filter.SetActiveField("var"); + filter.SetIsoValue(2.0); + vtkm::cont::DataSet outputSingleCell = filter.Execute(rectilinearDataset); + auto coordinates = outputSingleCell.GetCoordinateSystem() + .GetData() + .AsArrayHandle>(); + + VTKM_TEST_ASSERT(outputSingleCell.GetNumberOfPoints() == 3, + "Wrong number of points in rectilinear contour"); + VTKM_TEST_ASSERT(outputSingleCell.GetNumberOfCells() == 1, + "Wrong number of cells in rectilinear contour"); + VTKM_TEST_ASSERT(outputSingleCell.GetCellSet().GetCellShape(0) == vtkm::CELL_SHAPE_TRIANGLE, + "Wrong contour cell shape"); + + auto expectedCoordinates = + vtkm::cont::make_ArrayHandle({ vtkm::Vec3f{ 10.0f, -10.0f, 9.66341f }, + vtkm::Vec3f{ 9.30578f, -10.0f, 10.0f }, + vtkm::Vec3f{ 10.0f, -9.78842f, 10.0f } }); + VTKM_TEST_ASSERT(test_equal_ArrayHandles(coordinates, expectedCoordinates), + "Wrong contour coordinates"); + + // Generating normals triggers a different worklet for Flying Edges pass 4, + // But it should not change anything on the contour itself. + filter.SetGenerateNormals(true); + vtkm::cont::DataSet outputNormals = filter.Execute(rectilinearDataset); + coordinates = outputSingleCell.GetCoordinateSystem() + .GetData() + .AsArrayHandle>(); + VTKM_TEST_ASSERT(test_equal_ArrayHandles(coordinates, expectedCoordinates), + "Wrong contour coordinates"); + + // Full contour + filter.SetIsoValue(3.0); + filter.SetGenerateNormals(false); + vtkm::cont::DataSet output = filter.Execute(rectilinearDataset); + + VTKM_TEST_ASSERT(output.GetNumberOfPoints() == 93, + "Wrong number of points in rectilinear contour"); + VTKM_TEST_ASSERT(output.GetNumberOfCells() == 144, + "Wrong number of cells in rectilinear contour"); + } + void operator()() const { this->TestContourUniformGrid(72); @@ -220,6 +260,10 @@ public: this->TestContourWedges(); this->TestContourWedges(); + this->TestNonUniformStructured(); + this->TestNonUniformStructured(); + this->TestNonUniformStructured(); + this->TestUnsupportedFlyingEdges(); } diff --git a/vtkm/filter/contour/worklet/ContourFlyingEdges.h b/vtkm/filter/contour/worklet/ContourFlyingEdges.h index 85dabcf78..1a77c014a 100644 --- a/vtkm/filter/contour/worklet/ContourFlyingEdges.h +++ b/vtkm/filter/contour/worklet/ContourFlyingEdges.h @@ -67,13 +67,14 @@ public: // Filter called without normals generation template vtkm::cont::CellSetSingleType<> Run( const std::vector& isovalues, const vtkm::cont::CellSetStructured<3>& cells, - const vtkm::cont::ArrayHandleUniformPointCoordinates& coordinateSystem, + const CoordsType& coordinateSystem, const vtkm::cont::ArrayHandle& input, vtkm::cont::ArrayHandle, StorageTagVertices>& vertices) { @@ -87,6 +88,7 @@ public: // Filter called with normals generation template Run( const std::vector& isovalues, const vtkm::cont::CellSetStructured<3>& cells, - const vtkm::cont::ArrayHandleUniformPointCoordinates& coordinateSystem, + const CoordsType& coordinateSystem, const vtkm::cont::ArrayHandle& input, vtkm::cont::ArrayHandle, StorageTagVertices>& vertices, vtkm::cont::ArrayHandle, StorageTagNormals>& normals) diff --git a/vtkm/filter/contour/worklet/contour/FlyingEdges.h b/vtkm/filter/contour/worklet/contour/FlyingEdges.h index b9f605124..e5bd379c8 100644 --- a/vtkm/filter/contour/worklet/contour/FlyingEdges.h +++ b/vtkm/filter/contour/worklet/contour/FlyingEdges.h @@ -41,6 +41,7 @@ vtkm::Id extend_by(vtkm::cont::ArrayHandle& handle, vtkm::Id size) //---------------------------------------------------------------------------- template vtkm::cont::CellSetSingleType<> execute( const vtkm::cont::CellSetStructured<3>& cells, - const vtkm::cont::ArrayHandleUniformPointCoordinates& coordinateSystem, + const CoordsType coordinateSystem, const std::vector& isovalues, const vtkm::cont::ArrayHandle& inputField, vtkm::cont::ArrayHandle, StorageTagVertices>& points, @@ -56,18 +57,10 @@ vtkm::cont::CellSetSingleType<> execute( vtkm::worklet::contour::CommonState& sharedState) { vtkm::cont::Invoker invoke; - - vtkm::Vec3f origin, spacing; - { //extract out the origin and spacing as these are needed for Pass4 to properly - //interpolate the new points - auto portal = coordinateSystem.ReadPortal(); - origin = portal.GetOrigin(); - spacing = portal.GetSpacing(); - } auto pdims = cells.GetPointDimensions(); vtkm::cont::ArrayHandle edgeCases; - edgeCases.Allocate(coordinateSystem.GetNumberOfValues()); + edgeCases.Allocate(coordinateSystem.GetData().GetNumberOfValues()); vtkm::cont::CellSetStructured<2> metaDataMesh2D; vtkm::cont::ArrayHandle metaDataLinearSums; //per point of metaDataMesh @@ -158,8 +151,7 @@ vtkm::cont::CellSetSingleType<> execute( { VTKM_LOG_SCOPE(vtkm::cont::LogLevel::Perf, "FlyingEdges Pass4"); - launchComputePass4 pass4( - pdims, origin, spacing, multiContourCellOffset, multiContourPointOffset); + auto pass4 = launchComputePass4(pdims, multiContourCellOffset, multiContourPointOffset); detail::extend_by(points, newPointSize); if (sharedState.GenerateNormals) @@ -171,6 +163,7 @@ vtkm::cont::CellSetSingleType<> execute( pass4, newPointSize, isoval, + coordinateSystem, inputField, edgeCases, metaDataMesh2D, diff --git a/vtkm/filter/contour/worklet/contour/FlyingEdgesPass4.h b/vtkm/filter/contour/worklet/contour/FlyingEdgesPass4.h index 9905e7bdb..fe66dac49 100644 --- a/vtkm/filter/contour/worklet/contour/FlyingEdgesPass4.h +++ b/vtkm/filter/contour/worklet/contour/FlyingEdgesPass4.h @@ -28,20 +28,14 @@ namespace flying_edges struct launchComputePass4 { vtkm::Id3 PointDims; - vtkm::Vec3f Origin; - vtkm::Vec3f Spacing; vtkm::Id CellWriteOffset; vtkm::Id PointWriteOffset; 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) , CellWriteOffset(multiContourCellOffset) , PointWriteOffset(multiContourPointOffset) { @@ -49,6 +43,7 @@ struct launchComputePass4 template & inputField, vtkm::cont::ArrayHandle edgeCases, vtkm::cont::CellSetStructured<2>& metaDataMesh2D, @@ -71,12 +67,8 @@ struct launchComputePass4 vtkm::cont::Invoker invoke(device); if (sharedState.GenerateNormals) { - ComputePass4XWithNormals worklet4(isoval, - this->PointDims, - this->Origin, - this->Spacing, - this->CellWriteOffset, - this->PointWriteOffset); + ComputePass4XWithNormals worklet4( + isoval, this->PointDims, this->CellWriteOffset, this->PointWriteOffset); invoke(worklet4, metaDataMesh2D, metaDataSums, @@ -84,6 +76,7 @@ struct launchComputePass4 metaDataMax, metaDataNumTris, edgeCases, + coordinateSystem, inputField, triangle_topology, sharedState.InterpolationEdgeIds, @@ -94,12 +87,8 @@ struct launchComputePass4 } else { - ComputePass4X worklet4(isoval, - this->PointDims, - this->Origin, - this->Spacing, - this->CellWriteOffset, - this->PointWriteOffset); + ComputePass4X worklet4( + isoval, this->PointDims, this->CellWriteOffset, this->PointWriteOffset); invoke(worklet4, metaDataMesh2D, metaDataSums, @@ -107,6 +96,7 @@ struct launchComputePass4 metaDataMax, metaDataNumTris, edgeCases, + coordinateSystem, inputField, triangle_topology, sharedState.InterpolationEdgeIds, @@ -120,6 +110,7 @@ struct launchComputePass4 template & inputField, vtkm::cont::ArrayHandle edgeCases, vtkm::cont::CellSetStructured<2>& metaDataMesh2D, @@ -157,11 +149,8 @@ struct launchComputePass4 sharedState.CellIdMap); //This needs to be done on array handle view ( start = this->PointWriteOffset, len = newPointSize) - ComputePass5Y worklet5(this->PointDims, - this->Origin, - this->Spacing, - this->PointWriteOffset, - sharedState.GenerateNormals); + ComputePass5Y worklet5(this->PointDims, this->PointWriteOffset, sharedState.GenerateNormals); + invoke(worklet5, vtkm::cont::make_ArrayHandleView( sharedState.InterpolationEdgeIds, this->PointWriteOffset, newPointSize), @@ -169,6 +158,7 @@ struct launchComputePass4 sharedState.InterpolationWeights, this->PointWriteOffset, newPointSize), vtkm::cont::make_ArrayHandleView(points, this->PointWriteOffset, newPointSize), inputField, + coordinateSystem, normals); return true; diff --git a/vtkm/filter/contour/worklet/contour/FlyingEdgesPass4X.h b/vtkm/filter/contour/worklet/contour/FlyingEdgesPass4X.h index 99cf6550a..f6cddc9b8 100644 --- a/vtkm/filter/contour/worklet/contour/FlyingEdgesPass4X.h +++ b/vtkm/filter/contour/worklet/contour/FlyingEdgesPass4X.h @@ -32,9 +32,6 @@ struct ComputePass4X : public vtkm::worklet::WorkletVisitCellsWithPoints { vtkm::Id3 PointDims; - vtkm::Vec3f Origin; - vtkm::Vec3f Spacing; - T IsoValue; vtkm::Id CellWriteOffset; @@ -43,13 +40,9 @@ struct ComputePass4X : public vtkm::worklet::WorkletVisitCellsWithPoints 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) @@ -62,6 +55,7 @@ struct ComputePass4X : public vtkm::worklet::WorkletVisitCellsWithPoints FieldInPoint axis_maxs, WholeArrayIn cell_tri_count, WholeArrayIn edgeData, + WholeArrayIn coords, WholeArrayIn data, WholeArrayOut connectivity, WholeArrayOut edgeIds, @@ -69,13 +63,14 @@ struct ComputePass4X : public vtkm::worklet::WorkletVisitCellsWithPoints WholeArrayOut inputCellIds, WholeArrayOut points); using ExecutionSignature = - void(ThreadIndices, _2, _3, _4, _5, _6, _7, _8, _9, _10, _11, _12, WorkIndex); + void(ThreadIndices, _2, _3, _4, _5, _6, _7, _8, _9, _10, _11, _12, _13, WorkIndex); template VTKM_EXEC inline void Generate(const vtkm::Vec& boundaryStatus, const vtkm::Id3& ijk, const WholeDataField& field, const WholeIEdgeField& interpolatedEdgeIds, const WholeWeightField& weights, + const WholeCoordsField coords, const WholePointField& points, const vtkm::Id4& startPos, const vtkm::Id3& incs, @@ -188,7 +187,7 @@ struct ComputePass4X : public vtkm::worklet::WorkletVisitCellsWithPoints interpolatedEdgeIds.Set(writeIndex, pos); weights.Set(writeIndex, static_cast(t)); - auto coord = this->InterpolateCoordinate(t, ijk, ijk + vtkm::Id3{ 1, 0, 0 }); + auto coord = this->InterpolateCoordinate(coords, t, ijk, ijk + vtkm::Id3{ 1, 0, 0 }); points.Set(writeIndex, coord); } if (edgeUses[4]) @@ -201,7 +200,7 @@ struct ComputePass4X : public vtkm::worklet::WorkletVisitCellsWithPoints interpolatedEdgeIds.Set(writeIndex, pos); weights.Set(writeIndex, static_cast(t)); - auto coord = this->InterpolateCoordinate(t, ijk, ijk + vtkm::Id3{ 0, 1, 0 }); + auto coord = this->InterpolateCoordinate(coords, t, ijk, ijk + vtkm::Id3{ 0, 1, 0 }); points.Set(writeIndex, coord); } if (edgeUses[8]) @@ -214,7 +213,7 @@ struct ComputePass4X : public vtkm::worklet::WorkletVisitCellsWithPoints interpolatedEdgeIds.Set(writeIndex, pos); weights.Set(writeIndex, static_cast(t)); - auto coord = this->InterpolateCoordinate(t, ijk, ijk + vtkm::Id3{ 0, 0, 1 }); + auto coord = this->InterpolateCoordinate(coords, t, ijk, ijk + vtkm::Id3{ 0, 0, 1 }); points.Set(writeIndex, coord); } } @@ -231,30 +230,30 @@ struct ComputePass4X : public vtkm::worklet::WorkletVisitCellsWithPoints const bool onZ = boundaryStatus[AxisToSum::zindex] & FlyingEdges3D::MaxBoundary; if (onX) //+x boundary { - 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, 5, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, coords, points); + this->InterpolateEdge(ijk, pos[0], incs, 9, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, coords, points); if (onY) //+x +y { - this->InterpolateEdge(ijk, pos[0], incs, 11, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points); + this->InterpolateEdge(ijk, pos[0], incs, 11, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, coords, points); } if (onZ) //+x +z { - this->InterpolateEdge(ijk, pos[0], incs, 7, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points); + this->InterpolateEdge(ijk, pos[0], incs, 7, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, coords, points); } } if (onY) //+y boundary { - 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); + this->InterpolateEdge(ijk, pos[0], incs, 1, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, coords, points); + this->InterpolateEdge(ijk, pos[0], incs, 10, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, coords, points); if (onZ) //+y +z boundary { - this->InterpolateEdge(ijk, pos[0], incs, 3, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points); + this->InterpolateEdge(ijk, pos[0], incs, 3, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, coords, points); } } if (onZ) //+z boundary { - 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); + this->InterpolateEdge(ijk, pos[0], incs, 2, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, coords, points); + this->InterpolateEdge(ijk, pos[0], incs, 6, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, coords, points); } // clang-format on } @@ -264,6 +263,7 @@ struct ComputePass4X : public vtkm::worklet::WorkletVisitCellsWithPoints template VTKM_EXEC inline void InterpolateEdge(const vtkm::Id3& ijk, vtkm::Id currentIdx, @@ -274,6 +274,7 @@ struct ComputePass4X : public vtkm::worklet::WorkletVisitCellsWithPoints const WholeField& field, const WholeIEdgeField& interpolatedEdgeIds, const WholeWeightField& weights, + const WholeCoordsField& coords, const WholePointField& points) const { using AxisToSum = SumXAxis; @@ -300,30 +301,49 @@ struct ComputePass4X : public vtkm::worklet::WorkletVisitCellsWithPoints T t = static_cast((this->IsoValue - s0) / (s1 - s0)); weights.Set(writeIndex, static_cast(t)); - auto coord = this->InterpolateCoordinate(t, ijk + offsets1, ijk + offsets2); + auto coord = this->InterpolateCoordinate(coords, t, ijk + offsets1, ijk + offsets2); points.Set(writeIndex, coord); } + // Fast interpolation method for uniform coordinates //---------------------------------------------------------------------------- - inline VTKM_EXEC vtkm::Vec3f InterpolateCoordinate(T t, - const vtkm::Id3& ijk0, - const vtkm::Id3& ijk1) const + inline VTKM_EXEC vtkm::Vec3f InterpolateCoordinate( + vtkm::internal::ArrayPortalUniformPointCoordinates coords, + T t, + const vtkm::Id3& ijk0, + const vtkm::Id3& ijk1) const { return vtkm::Vec3f( - this->Origin[0] + - this->Spacing[0] * + coords.GetOrigin()[0] + + coords.GetSpacing()[0] * (static_cast(ijk0[0]) + static_cast(t) * static_cast(ijk1[0] - ijk0[0])), - this->Origin[1] + - this->Spacing[1] * + coords.GetOrigin()[1] + + coords.GetSpacing()[1] * (static_cast(ijk0[1]) + static_cast(t) * static_cast(ijk1[1] - ijk0[1])), - this->Origin[2] + - this->Spacing[2] * + coords.GetOrigin()[2] + + coords.GetSpacing()[2] * (static_cast(ijk0[2]) + static_cast(t) * static_cast(ijk1[2] - ijk0[2]))); } + + // Interpolation for explicit coordinates + //---------------------------------------------------------------------------- + template + inline VTKM_EXEC vtkm::Vec3f InterpolateCoordinate(CoordsPortal coords, + T t, + const vtkm::Id3& ijk0, + const vtkm::Id3& ijk1) const + { + return (1.0f - static_cast(t)) * + coords.Get(ijk0[0] + this->PointDims[0] * ijk0[1] + + this->PointDims[0] * this->PointDims[1] * ijk0[2]) + + static_cast(t) * + coords.Get(ijk1[0] + this->PointDims[0] * ijk1[1] + + this->PointDims[0] * this->PointDims[1] * ijk1[2]); + } }; } } diff --git a/vtkm/filter/contour/worklet/contour/FlyingEdgesPass4XWithNormals.h b/vtkm/filter/contour/worklet/contour/FlyingEdgesPass4XWithNormals.h index aea23acc7..678052299 100644 --- a/vtkm/filter/contour/worklet/contour/FlyingEdgesPass4XWithNormals.h +++ b/vtkm/filter/contour/worklet/contour/FlyingEdgesPass4XWithNormals.h @@ -31,9 +31,6 @@ struct ComputePass4XWithNormals : public vtkm::worklet::WorkletVisitCellsWithPoi { vtkm::Id3 PointDims; - vtkm::Vec3f Origin; - vtkm::Vec3f Spacing; - T IsoValue; vtkm::Id CellWriteOffset; @@ -42,13 +39,9 @@ struct ComputePass4XWithNormals : public vtkm::worklet::WorkletVisitCellsWithPoi 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) , IsoValue(value) , CellWriteOffset(multiContourCellOffset) , PointWriteOffset(multiContourPointOffset) @@ -61,6 +54,7 @@ struct ComputePass4XWithNormals : public vtkm::worklet::WorkletVisitCellsWithPoi FieldInPoint axis_maxs, WholeArrayIn cell_tri_count, WholeArrayIn edgeData, + WholeArrayIn coords, WholeArrayIn data, WholeArrayOut connectivity, WholeArrayOut edgeIds, @@ -69,13 +63,14 @@ struct ComputePass4XWithNormals : public vtkm::worklet::WorkletVisitCellsWithPoi WholeArrayOut points, WholeArrayOut normals); using ExecutionSignature = - void(ThreadIndices, _2, _3, _4, _5, _6, _7, _8, _9, _10, _11, _12, _13, WorkIndex); + void(ThreadIndices, _2, _3, _4, _5, _6, _7, _8, _9, _10, _11, _12, _13, _14, WorkIndex); template VTKM_EXEC inline void Generate(const vtkm::Vec& boundaryStatus, @@ -169,6 +167,7 @@ struct ComputePass4XWithNormals : public vtkm::worklet::WorkletVisitCellsWithPoi const WholeDataField& field, const WholeIEdgeField& interpolatedEdgeIds, const WholeWeightField& weights, + const WholeCoordsField coords, const WholePointField& points, const WholeNormalField& normals, const vtkm::Id4& startPos, @@ -197,7 +196,7 @@ struct ComputePass4XWithNormals : public vtkm::worklet::WorkletVisitCellsWithPoi weights.Set(writeIndex, static_cast(t)); auto ijk1 = ijk + vtkm::Id3{ 1, 0, 0 }; - auto coord = this->InterpolateCoordinate(t, ijk, ijk1); + auto coord = this->InterpolateCoordinate(coords, t, ijk, ijk1); points.Set(writeIndex, coord); //gradient generation @@ -216,7 +215,7 @@ struct ComputePass4XWithNormals : public vtkm::worklet::WorkletVisitCellsWithPoi weights.Set(writeIndex, static_cast(t)); auto ijk1 = ijk + vtkm::Id3{ 0, 1, 0 }; - auto coord = this->InterpolateCoordinate(t, ijk, ijk1); + auto coord = this->InterpolateCoordinate(coords, t, ijk, ijk1); points.Set(writeIndex, coord); //gradient generation @@ -235,7 +234,7 @@ struct ComputePass4XWithNormals : public vtkm::worklet::WorkletVisitCellsWithPoi weights.Set(writeIndex, static_cast(t)); auto ijk1 = ijk + vtkm::Id3{ 0, 0, 1 }; - auto coord = this->InterpolateCoordinate(t, ijk, ijk1); + auto coord = this->InterpolateCoordinate(coords, t, ijk, ijk1); points.Set(writeIndex, coord); //gradient generation @@ -258,38 +257,38 @@ struct ComputePass4XWithNormals : public vtkm::worklet::WorkletVisitCellsWithPoi if (onX) //+x boundary { this->InterpolateEdge( - fullyInterior, ijk, pos[0], incs, 5, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points, normals); + fullyInterior, ijk, pos[0], incs, 5, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, coords, points, normals); this->InterpolateEdge( - fullyInterior, ijk, pos[0], incs, 9, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points, normals); + fullyInterior, ijk, pos[0], incs, 9, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, coords, points, normals); if (onY) //+x +y { this->InterpolateEdge( - fullyInterior, ijk, pos[0], incs, 11, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points, normals); + fullyInterior, ijk, pos[0], incs, 11, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, coords, points, normals); } if (onZ) //+x +z { this->InterpolateEdge( - fullyInterior, ijk, pos[0], incs, 7, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points, normals); + fullyInterior, ijk, pos[0], incs, 7, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, coords, points, normals); } } if (onY) //+y boundary { this->InterpolateEdge( - fullyInterior, ijk, pos[0], incs, 1, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points, normals); + fullyInterior, ijk, pos[0], incs, 1, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, coords, points, normals); this->InterpolateEdge( - fullyInterior, ijk, pos[0], incs, 10, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points, normals); + fullyInterior, ijk, pos[0], incs, 10, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, coords, points, normals); if (onZ) //+y +z boundary { this->InterpolateEdge( - fullyInterior, ijk, pos[0], incs, 3, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points, normals); + fullyInterior, ijk, pos[0], incs, 3, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, coords, points, normals); } } if (onZ) //+z boundary { this->InterpolateEdge( - fullyInterior, ijk, pos[0], incs, 2, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points, normals); + fullyInterior, ijk, pos[0], incs, 2, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, coords, points, normals); this->InterpolateEdge( - fullyInterior, ijk, pos[0], incs, 6, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, points, normals); + fullyInterior, ijk, pos[0], incs, 6, edgeUses, edgeIds, field, interpolatedEdgeIds, weights, coords, points, normals); } // clang-format on } @@ -300,6 +299,7 @@ struct ComputePass4XWithNormals : public vtkm::worklet::WorkletVisitCellsWithPoi typename WholeIEdgeField, typename WholeWeightField, typename WholePointField, + typename WholeCoordsField, typename WholeNormalField> VTKM_EXEC inline void InterpolateEdge(bool fullyInterior, const vtkm::Id3& ijk, @@ -311,6 +311,7 @@ struct ComputePass4XWithNormals : public vtkm::worklet::WorkletVisitCellsWithPoi const WholeField& field, const WholeIEdgeField& interpolatedEdgeIds, const WholeWeightField& weights, + const WholeCoordsField& coords, const WholePointField& points, const WholeNormalField& normals) const { @@ -338,36 +339,55 @@ struct ComputePass4XWithNormals : public vtkm::worklet::WorkletVisitCellsWithPoi T t = static_cast((this->IsoValue - s0) / (s1 - s0)); weights.Set(writeIndex, static_cast(t)); - auto coord = this->InterpolateCoordinate(t, ijk + offsets1, ijk + offsets2); + auto coord = this->InterpolateCoordinate(coords, t, ijk + offsets1, ijk + offsets2); points.Set(writeIndex, coord); - auto g0 = this->ComputeGradient(fullyInterior, ijk + offsets1, incs, iEdge[0], field); auto g1 = this->ComputeGradient(fullyInterior, ijk + offsets2, incs, iEdge[1], field); + auto g0 = this->ComputeGradient(fullyInterior, ijk + offsets1, incs, iEdge[0], field); g1 = g0 + (t * (g1 - g0)); normals.Set(writeIndex, vtkm::Normal(g1)); } + // Fast interpolation method for uniform coordinates //---------------------------------------------------------------------------- - inline VTKM_EXEC vtkm::Vec3f InterpolateCoordinate(T t, - const vtkm::Id3& ijk0, - const vtkm::Id3& ijk1) const + inline VTKM_EXEC vtkm::Vec3f InterpolateCoordinate( + const vtkm::internal::ArrayPortalUniformPointCoordinates& coords, + T t, + const vtkm::Id3& ijk0, + const vtkm::Id3& ijk1) const { return vtkm::Vec3f( - this->Origin[0] + - this->Spacing[0] * + coords.GetOrigin()[0] + + coords.GetSpacing()[0] * (static_cast(ijk0[0]) + static_cast(t) * static_cast(ijk1[0] - ijk0[0])), - this->Origin[1] + - this->Spacing[1] * + coords.GetOrigin()[1] + + coords.GetSpacing()[1] * (static_cast(ijk0[1]) + static_cast(t) * static_cast(ijk1[1] - ijk0[1])), - this->Origin[2] + - this->Spacing[2] * + coords.GetOrigin()[2] + + coords.GetSpacing()[2] * (static_cast(ijk0[2]) + static_cast(t) * static_cast(ijk1[2] - ijk0[2]))); } + // Interpolation for explicit coordinates + //---------------------------------------------------------------------------- + template + inline VTKM_EXEC vtkm::Vec3f InterpolateCoordinate(const CoordsPortal& coords, + T t, + const vtkm::Id3& ijk0, + const vtkm::Id3& ijk1) const + { + return (1.0f - static_cast(t)) * + coords.Get(ijk0[0] + this->PointDims[0] * ijk0[1] + + this->PointDims[0] * this->PointDims[1] * ijk0[2]) + + static_cast(t) * + coords.Get(ijk1[0] + this->PointDims[0] * ijk1[1] + + this->PointDims[0] * this->PointDims[1] * ijk1[2]); + } + //---------------------------------------------------------------------------- template VTKM_EXEC vtkm::Vec3f ComputeGradient(bool fullyInterior, diff --git a/vtkm/filter/contour/worklet/contour/FlyingEdgesPass4Y.h b/vtkm/filter/contour/worklet/contour/FlyingEdgesPass4Y.h index 05e73ba13..08a16f3b4 100644 --- a/vtkm/filter/contour/worklet/contour/FlyingEdgesPass4Y.h +++ b/vtkm/filter/contour/worklet/contour/FlyingEdgesPass4Y.h @@ -276,17 +276,11 @@ struct ComputePass4Y : public vtkm::worklet::WorkletVisitCellsWithPoints template struct ComputePass5Y : public vtkm::worklet::WorkletMapField { - - vtkm::internal::ArrayPortalUniformPointCoordinates Coordinates; + vtkm::Id3 PointDims; 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) + ComputePass5Y(const vtkm::Id3& pdims, vtkm::Id normalWriteOffset, bool generateNormals) + : PointDims(pdims) , NormalWriteOffset(normalWriteOffset) { if (!generateNormals) @@ -299,20 +293,25 @@ struct ComputePass5Y : public vtkm::worklet::WorkletMapField FieldIn interpWeight, FieldOut points, WholeArrayIn field, + WholeArrayIn coords, WholeArrayOut normals); - using ExecutionSignature = void(_1, _2, _3, _4, _5, WorkIndex); + using ExecutionSignature = void(_1, _2, _3, _4, _5, _6, WorkIndex); - template + template VTKM_EXEC void operator()(const vtkm::Id2& interpEdgeIds, vtkm::FloatDefault weight, vtkm::Vec& outPoint, const WholeInputField& field, + const WholeCoordsField& coords, WholeNormalField& normals, vtkm::Id oidx) const { { - vtkm::Vec3f point1 = this->Coordinates.Get(interpEdgeIds[0]); - vtkm::Vec3f point2 = this->Coordinates.Get(interpEdgeIds[1]); + vtkm::Vec3f point1 = coords.Get(interpEdgeIds[0]); + vtkm::Vec3f point2 = coords.Get(interpEdgeIds[1]); outPoint = vtkm::Lerp(point1, point2, weight); } @@ -320,15 +319,13 @@ struct ComputePass5Y : public vtkm::worklet::WorkletMapField if (this->NormalWriteOffset >= 0) { vtkm::Vec 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::Id3 ijk{ interpEdgeIds[0] % this->PointDims[0], + (interpEdgeIds[0] / this->PointDims[0]) % this->PointDims[1], + interpEdgeIds[0] / (this->PointDims[0] * this->PointDims[1]) }; vtkm::worklet::gradient::StructuredPointGradient gradient; - vtkm::exec::BoundaryState boundary(ijk, dims); - vtkm::exec::FieldNeighborhood - coord_neighborhood(this->Coordinates, boundary); + vtkm::exec::BoundaryState boundary(ijk, this->PointDims); + vtkm::exec::FieldNeighborhood coord_neighborhood(coords, boundary); vtkm::exec::FieldNeighborhood field_neighborhood(field, boundary); @@ -337,9 +334,9 @@ struct ComputePass5Y : public vtkm::worklet::WorkletMapField 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]) }; + boundary.IJK = vtkm::Id3{ interpEdgeIds[1] % this->PointDims[0], + (interpEdgeIds[1] / this->PointDims[0]) % this->PointDims[1], + interpEdgeIds[1] / (this->PointDims[0] * this->PointDims[1]) }; gradient(boundary, coord_neighborhood, field_neighborhood, g1); vtkm::Vec3f n = vtkm::Lerp(g0, g1, weight);