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

Merge topic 'fe-structured'

Browse Source 
1538fc02f Implement Flying Edges for all structured CellSets

Acked-by: Kitware Robot <kwrobot@kitware.com>
Acked-by: Kenneth Moreland <morelandkd@ornl.gov>
Merge-request: !3056
This commit is contained in:
Louis Gombert 2023-05-12 12:06:20 +00:00 committed by Kitware Robot
commit 57abd08684
11 changed files with 212 additions and 164 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
docs/changelog/flying-edges-structured.md Normal file
View File

@ -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.

10
vtkm/filter/clean_grid/testing/UnitTestCleanGrid.cxx
View File

@ -13,7 +13,7 @@
#include <vtkm/cont/ArrayCopy.h>
#include <vtkm/cont/testing/MakeTestDataSet.h>
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/filter/contour/Contour.h>
#include <vtkm/filter/contour/ContourMarchingCells.h>
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<vtkm::Vec3f> 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");

7
vtkm/filter/contour/Contour.cxx
View File

@ -34,11 +34,8 @@ vtkm::cont::DataSet Contour::DoExecute(const vtkm::cont::DataSet& inDataSet)
auto inCoords = inDataSet.GetCoordinateSystem(this->GetActiveCoordinateSystemIndex()).GetData();
std::unique_ptr<vtkm::filter::contour::AbstractContour> implementation;
// For now, Flying Edges is only used for 3D Structured CellSets,
// using Uniform coordinates.
if (inCellSet.template IsType<vtkm::cont::CellSetStructured<3>>() &&
inCoords.template IsType<
vtkm::cont::ArrayHandle<vtkm::Vec3f, vtkm::cont::StorageTagUniformPoints>>())
// Flying Edges is only used for 3D Structured CellSets
if (inCellSet.template IsType<vtkm::cont::CellSetStructured<3>>())
{
VTKM_LOG_S(vtkm::cont::LogLevel::Info, "Using flying edges");
implementation.reset(new vtkm::filter::contour::ContourFlyingEdges);

22
vtkm/filter/contour/ContourFlyingEdges.cxx
View File

@ -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<vtkm::cont::CellSetStructured<3>>() ||
!inCoords.GetData()
.template IsType<
vtkm::cont::ArrayHandle<vtkm::Vec3f, vtkm::cont::StorageTagUniformPoints>>())
if (!inCellSet.template IsType<vtkm::cont::CellSetStructured<3>>())
{
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<vtkm::cont::ArrayHandleUniformPointCoordinates>(),
concrete,
vertices,
normals);
outputCells = worklet.Run(ivalues, inputCells, inCoords, concrete, vertices, normals);
}
else
{
outputCells = worklet.Run(
ivalues,
inputCells,
inCoords.GetData().AsArrayHandle<vtkm::cont::ArrayHandleUniformPointCoordinates>(),
concrete,
vertices);
outputCells = worklet.Run(ivalues, inputCells, inCoords, concrete, vertices);
}
};

68
vtkm/filter/contour/testing/UnitTestContourFilter.cxx
View File

@ -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 <typename ContourFilterType>
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::cont::ArrayHandle<vtkm::Vec3f>>();
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>({ 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::cont::ArrayHandle<vtkm::Vec3f>>();
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<vtkm::filter::contour::Contour>(72);
@ -220,6 +260,10 @@ public:
this->TestContourWedges<vtkm::filter::contour::Contour>();
this->TestContourWedges<vtkm::filter::contour::ContourMarchingCells>();
this->TestNonUniformStructured<vtkm::filter::contour::Contour>();
this->TestNonUniformStructured<vtkm::filter::contour::ContourFlyingEdges>();
this->TestNonUniformStructured<vtkm::filter::contour::ContourMarchingCells>();
this->TestUnsupportedFlyingEdges();
}

6
vtkm/filter/contour/worklet/ContourFlyingEdges.h
View File

@ -67,13 +67,14 @@ public:
// Filter called without normals generation
template <typename ValueType,
typename CoordsType,
typename StorageTagField,
typename CoordinateType,
typename StorageTagVertices>
vtkm::cont::CellSetSingleType<> Run(
const std::vector<ValueType>& isovalues,
const vtkm::cont::CellSetStructured<3>& cells,
const vtkm::cont::ArrayHandleUniformPointCoordinates& coordinateSystem,
const CoordsType& coordinateSystem,
const vtkm::cont::ArrayHandle<ValueType, StorageTagField>& input,
vtkm::cont::ArrayHandle<vtkm::Vec<CoordinateType, 3>, StorageTagVertices>& vertices)
{
@ -87,6 +88,7 @@ public:
// Filter called with normals generation
template <typename ValueType,
typename CoordsType,
typename StorageTagField,
typename CoordinateType,
typename StorageTagVertices,
@ -94,7 +96,7 @@ public:
vtkm::cont::CellSetSingleType<> Run(
const std::vector<ValueType>& isovalues,
const vtkm::cont::CellSetStructured<3>& cells,
const vtkm::cont::ArrayHandleUniformPointCoordinates& coordinateSystem,
const CoordsType& coordinateSystem,
const vtkm::cont::ArrayHandle<ValueType, StorageTagField>& input,
vtkm::cont::ArrayHandle<vtkm::Vec<CoordinateType, 3>, StorageTagVertices>& vertices,
vtkm::cont::ArrayHandle<vtkm::Vec<CoordinateType, 3>, StorageTagNormals>& normals)

17
vtkm/filter/contour/worklet/contour/FlyingEdges.h
View File

@ -41,6 +41,7 @@ vtkm::Id extend_by(vtkm::cont::ArrayHandle<T, S>& handle, vtkm::Id size)
//----------------------------------------------------------------------------
template <typename ValueType,
typename CoordsType,
typename StorageTagField,
typename StorageTagVertices,
typename StorageTagNormals,
@ -48,7 +49,7 @@ template <typename ValueType,
typename NormalType>
vtkm::cont::CellSetSingleType<> execute(
const vtkm::cont::CellSetStructured<3>& cells,
const vtkm::cont::ArrayHandleUniformPointCoordinates& coordinateSystem,
const CoordsType coordinateSystem,
const std::vector<ValueType>& isovalues,
const vtkm::cont::ArrayHandle<ValueType, StorageTagField>& inputField,
vtkm::cont::ArrayHandle<vtkm::Vec<CoordinateType, 3>, 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<vtkm::UInt8> edgeCases;
edgeCases.Allocate(coordinateSystem.GetNumberOfValues());
edgeCases.Allocate(coordinateSystem.GetData().GetNumberOfValues());
vtkm::cont::CellSetStructured<2> metaDataMesh2D;
vtkm::cont::ArrayHandle<vtkm::Id> 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,

36
vtkm/filter/contour/worklet/contour/FlyingEdgesPass4.h
View File

@ -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 <typename DeviceAdapterTag,
typename T,
typename CoordsType,
typename StorageTagField,
typename MeshSums,
typename PointType,
@ -56,6 +51,7 @@ struct launchComputePass4
VTKM_CONT bool LaunchXAxis(DeviceAdapterTag device,
vtkm::Id vtkmNotUsed(newPointSize),
T isoval,
CoordsType coordinateSystem,
const vtkm::cont::ArrayHandle<T, StorageTagField>& inputField,
vtkm::cont::ArrayHandle<vtkm::UInt8> edgeCases,
vtkm::cont::CellSetStructured<2>& metaDataMesh2D,
@ -71,12 +67,8 @@ struct launchComputePass4
vtkm::cont::Invoker invoke(device);
if (sharedState.GenerateNormals)
{
ComputePass4XWithNormals<T> worklet4(isoval,
this->PointDims,
this->Origin,
this->Spacing,
this->CellWriteOffset,
this->PointWriteOffset);
ComputePass4XWithNormals<T> 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<T> worklet4(isoval,
this->PointDims,
this->Origin,
this->Spacing,
this->CellWriteOffset,
this->PointWriteOffset);
ComputePass4X<T> 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 <typename DeviceAdapterTag,
typename T,
typename CoordsType,
typename StorageTagField,
typename MeshSums,
typename PointType,
@ -127,6 +118,7 @@ struct launchComputePass4
VTKM_CONT bool LaunchYAxis(DeviceAdapterTag device,
vtkm::Id newPointSize,
T isoval,
CoordsType coordinateSystem,
const vtkm::cont::ArrayHandle<T, StorageTagField>& inputField,
vtkm::cont::ArrayHandle<vtkm::UInt8> 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<T> worklet5(this->PointDims,
this->Origin,
this->Spacing,
this->PointWriteOffset,
sharedState.GenerateNormals);
ComputePass5Y<T> 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;

80
vtkm/filter/contour/worklet/contour/FlyingEdgesPass4X.h
View File

@ -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 <typename ThreadIndices,
typename FieldInPointId3,
typename FieldInPointId,
typename WholeTriField,
typename WholeEdgeField,
typename WholeCoordsField,
typename WholeDataField,
typename WholeConnField,
typename WholeEdgeIdField,
@ -88,6 +83,7 @@ struct ComputePass4X : public vtkm::worklet::WorkletVisitCellsWithPoints
const FieldInPointId& axis_maxs,
const WholeTriField& cellTriCount,
const WholeEdgeField& edges,
const WholeCoordsField& coords,
const WholeDataField& field,
const WholeConnField& conn,
const WholeEdgeIdField& interpolatedEdgeIds,
@ -141,6 +137,7 @@ struct ComputePass4X : public vtkm::worklet::WorkletVisitCellsWithPoints
field,
interpolatedEdgeIds,
weights,
coords,
points,
state.startPos,
increments,
@ -158,12 +155,14 @@ struct ComputePass4X : public vtkm::worklet::WorkletVisitCellsWithPoints
template <typename WholeDataField,
typename WholeIEdgeField,
typename WholeWeightField,
typename WholeCoordsField,
typename WholePointField>
VTKM_EXEC inline void Generate(const vtkm::Vec<vtkm::UInt8, 3>& 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<vtkm::FloatDefault>(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<vtkm::FloatDefault>(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<vtkm::FloatDefault>(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 <typename WholeField,
typename WholeIEdgeField,
typename WholeWeightField,
typename WholeCoordsField,
typename WholePointField>
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<T>((this->IsoValue - s0) / (s1 - s0));
weights.Set(writeIndex, static_cast<vtkm::FloatDefault>(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<vtkm::FloatDefault>(ijk0[0]) +
static_cast<vtkm::FloatDefault>(t) * static_cast<vtkm::FloatDefault>(ijk1[0] - ijk0[0])),
this->Origin[1] +
this->Spacing[1] *
coords.GetOrigin()[1] +
coords.GetSpacing()[1] *
(static_cast<vtkm::FloatDefault>(ijk0[1]) +
static_cast<vtkm::FloatDefault>(t) * static_cast<vtkm::FloatDefault>(ijk1[1] - ijk0[1])),
this->Origin[2] +
this->Spacing[2] *
coords.GetOrigin()[2] +
coords.GetSpacing()[2] *
(static_cast<vtkm::FloatDefault>(ijk0[2]) +
static_cast<vtkm::FloatDefault>(t) *
static_cast<vtkm::FloatDefault>(ijk1[2] - ijk0[2])));
}
// Interpolation for explicit coordinates
//----------------------------------------------------------------------------
template <typename CoordsPortal>
inline VTKM_EXEC vtkm::Vec3f InterpolateCoordinate(CoordsPortal coords,
T t,
const vtkm::Id3& ijk0,
const vtkm::Id3& ijk1) const
{
return (1.0f - static_cast<vtkm::FloatDefault>(t)) *
coords.Get(ijk0[0] + this->PointDims[0] * ijk0[1] +
this->PointDims[0] * this->PointDims[1] * ijk0[2]) +
static_cast<vtkm::FloatDefault>(t) *
coords.Get(ijk1[0] + this->PointDims[0] * ijk1[1] +
this->PointDims[0] * this->PointDims[1] * ijk1[2]);
}
};
}
}

82
vtkm/filter/contour/worklet/contour/FlyingEdgesPass4XWithNormals.h
View File

@ -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 <typename ThreadIndices,
typename FieldInPointId3,
typename FieldInPointId,
typename WholeTriField,
typename WholeEdgeField,
typename WholeCoordsField,
typename WholeDataField,
typename WholeConnField,
typename WholeEdgeIdField,
@ -89,6 +84,7 @@ struct ComputePass4XWithNormals : public vtkm::worklet::WorkletVisitCellsWithPoi
const FieldInPointId& axis_maxs,
const WholeTriField& cellTriCount,
const WholeEdgeField& edges,
const WholeCoordsField& coords,
const WholeDataField& field,
const WholeConnField& conn,
const WholeEdgeIdField& interpolatedEdgeIds,
@ -144,6 +140,7 @@ struct ComputePass4XWithNormals : public vtkm::worklet::WorkletVisitCellsWithPoi
field,
interpolatedEdgeIds,
weights,
coords,
points,
normals,
state.startPos,
@ -162,6 +159,7 @@ struct ComputePass4XWithNormals : public vtkm::worklet::WorkletVisitCellsWithPoi
template <typename WholeDataField,
typename WholeIEdgeField,
typename WholeWeightField,
typename WholeCoordsField,
typename WholePointField,
typename WholeNormalField>
VTKM_EXEC inline void Generate(const vtkm::Vec<vtkm::UInt8, 3>& 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<vtkm::FloatDefault>(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<vtkm::FloatDefault>(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<vtkm::FloatDefault>(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<T>((this->IsoValue - s0) / (s1 - s0));
weights.Set(writeIndex, static_cast<vtkm::FloatDefault>(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<vtkm::FloatDefault>(ijk0[0]) +
static_cast<vtkm::FloatDefault>(t) * static_cast<vtkm::FloatDefault>(ijk1[0] - ijk0[0])),
this->Origin[1] +
this->Spacing[1] *
coords.GetOrigin()[1] +
coords.GetSpacing()[1] *
(static_cast<vtkm::FloatDefault>(ijk0[1]) +
static_cast<vtkm::FloatDefault>(t) * static_cast<vtkm::FloatDefault>(ijk1[1] - ijk0[1])),
this->Origin[2] +
this->Spacing[2] *
coords.GetOrigin()[2] +
coords.GetSpacing()[2] *
(static_cast<vtkm::FloatDefault>(ijk0[2]) +
static_cast<vtkm::FloatDefault>(t) *
static_cast<vtkm::FloatDefault>(ijk1[2] - ijk0[2])));
}
// Interpolation for explicit coordinates
//----------------------------------------------------------------------------
template <typename CoordsPortal>
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<vtkm::FloatDefault>(t)) *
coords.Get(ijk0[0] + this->PointDims[0] * ijk0[1] +
this->PointDims[0] * this->PointDims[1] * ijk0[2]) +
static_cast<vtkm::FloatDefault>(t) *
coords.Get(ijk1[0] + this->PointDims[0] * ijk1[1] +
this->PointDims[0] * this->PointDims[1] * ijk1[2]);
}
//----------------------------------------------------------------------------
template <typename WholeDataField>
VTKM_EXEC vtkm::Vec3f ComputeGradient(bool fullyInterior,

43
vtkm/filter/contour/worklet/contour/FlyingEdgesPass4Y.h
View File

@ -276,17 +276,11 @@ struct ComputePass4Y : public vtkm::worklet::WorkletVisitCellsWithPoints
template <typename T>
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 <typename PT, typename WholeInputField, typename WholeNormalField>
template <typename PT,
typename WholeInputField,
typename WholeNormalField,
typename WholeCoordsField>
VTKM_EXEC void operator()(const vtkm::Id2& interpEdgeIds,
vtkm::FloatDefault weight,
vtkm::Vec<PT, 3>& 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<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::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<vtkm::internal::ArrayPortalUniformPointCoordinates>
coord_neighborhood(this->Coordinates, boundary);
vtkm::exec::BoundaryState boundary(ijk, this->PointDims);
vtkm::exec::FieldNeighborhood<WholeCoordsField> coord_neighborhood(coords, boundary);
vtkm::exec::FieldNeighborhood<WholeInputField> 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);