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);