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

Add point merge capabilities to CleanGrid filter

Browse Source
This commit is contained in:
Kenneth Moreland 2019-01-28 09:24:48 -05:00
parent 5e2e0d9ff1
commit 5688375c99
16 changed files with 808 additions and 59 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

26
docs/changelog/point-merge.md Normal file
View File

@ -0,0 +1,26 @@
# Add point merge capabilities to CleanGrid filter
We have added a `PointMerge` worklet that uses a virtual grid approach to
identify nearby points. The worklet works by creating a very fine but
sparsely represented locator grid. It then groups points by grid bins and
finds those within a specified radius.
This functionality has been integrated into the `CleanGrid` filter. The
following flags have been added to `CleanGrid` to modify the behavior of
point merging.
* `Set`/`GetMergePoints` - a flag to turn on/off the merging of
duplicated coincident points. This extra operation will find points
spatially located near each other and merge them together.
* `Set`/`GetTolerance` - Defines the tolerance used when determining
whether two points are considered coincident. If the
`ToleranceIsAbsolute` flag is false (the default), then this tolerance
is scaled by the diagonal of the points. This parameter is only used
when merge points is on.
* `Set`/`GetToleranceIsAbsolute` - When ToleranceIsAbsolute is false (the
default) then the tolerance is scaled by the diagonal of the bounds of
the dataset. If true, then the tolerance is taken as the actual
distance to use. This parameter is only used when merge points is on.
* `Set`/`GetFastMerge` - When FastMerge is true (the default), some
corners are cut when computing coincident points. The point merge will
go faster but the tolerance will not be strictly followed.

40
vtkm/filter/CleanGrid.h
View File

@ -22,6 +22,7 @@
#include <vtkm/filter/FilterDataSet.h>
#include <vtkm/worklet/PointMerge.h>
#include <vtkm/worklet/RemoveUnusedPoints.h>
namespace vtkm
@ -52,10 +53,36 @@ public:
/// When the CompactPointFields flag is true, the filter will identify any
/// points that are not used by the topology. This is on by default.
///
VTKM_CONT
bool GetCompactPointFields() const { return this->CompactPointFields; }
VTKM_CONT
void SetCompactPointFields(bool flag) { this->CompactPointFields = flag; }
VTKM_CONT bool GetCompactPointFields() const { return this->CompactPointFields; }
VTKM_CONT void SetCompactPointFields(bool flag) { this->CompactPointFields = flag; }
/// When the MergePoints flag is true, the filter will identify any coincident
/// points and merge them together. The distance two points can be to considered
/// coincident is set with the tolerance flags. This is on by default.
///
VTKM_CONT bool GetMergePoints() const { return this->MergePoints; }
VTKM_CONT void SetMergePoints(bool flag) { this->MergePoints = flag; }
/// Defines the tolerance used when determining whether two points are considered
/// coincident. If the ToleranceIsAbsolute flag is false (the default), then this
/// tolerance is scaled by the diagonal of the points.
///
VTKM_CONT vtkm::Float64 GetTolerance() const { return this->Tolerance; }
VTKM_CONT void SetTolerance(vtkm::Float64 tolerance) { this->Tolerance = tolerance; }
/// When ToleranceIsAbsolute is false (the default) then the tolerance is scaled
/// by the diagonal of the bounds of the dataset. If true, then the tolerance is
/// taken as the actual distance to use.
///
VTKM_CONT bool GetToleranceIsAbsolute() const { return this->ToleranceIsAbsolute; }
VTKM_CONT void SetToleranceIsAbsolute(bool flag) { this->ToleranceIsAbsolute = flag; }
/// When FastMerge is true (the default), some corners are cut when computing
/// coincident points. The point merge will go faster but the tolerance will not
/// be strictly followed.
///
VTKM_CONT bool GetFastMerge() const { return this->FastMerge; }
VTKM_CONT void SetFastMerge(bool flag) { this->FastMerge = flag; }
template <typename Policy>
VTKM_CONT vtkm::cont::DataSet DoExecute(const vtkm::cont::DataSet& inData,
@ -73,8 +100,13 @@ public:
private:
bool CompactPointFields;
bool MergePoints;
vtkm::Float64 Tolerance;
bool ToleranceIsAbsolute;
bool FastMerge;
vtkm::worklet::RemoveUnusedPoints PointCompactor;
vtkm::worklet::PointMerge PointMerger;
};
}
} // namespace vtkm::filter

121
vtkm/filter/CleanGrid.hxx
View File

@ -18,6 +18,11 @@
// this software.
//============================================================================
#ifndef vtk_m_filter_CleanGrid_hxx
#define vtk_m_filter_CleanGrid_hxx
#include <vtkm/filter/CleanGrid.h>
#include <vtkm/worklet/CellDeepCopy.h>
#include <vtkm/worklet/RemoveUnusedPoints.h>
@ -30,6 +35,10 @@ namespace filter
inline VTKM_CONT CleanGrid::CleanGrid()
: CompactPointFields(true)
, MergePoints(true)
, Tolerance(1.0e-6)
, ToleranceIsAbsolute(false)
, FastMerge(true)
{
}
@ -41,11 +50,10 @@ inline VTKM_CONT vtkm::cont::DataSet CleanGrid::DoExecute(const vtkm::cont::Data
using VecId = std::vector<CellSetType>::size_type;
VecId numCellSets = static_cast<VecId>(inData.GetNumberOfCellSets());
std::vector<CellSetType> outputCellSets(numCellSets);
// Do a deep copy of the cells to new CellSetExplicit structures
for (VecId cellSetIndex = 0; cellSetIndex < numCellSets; cellSetIndex++)
for (VecId cellSetIndex = 0; cellSetIndex < numCellSets; ++cellSetIndex)
{
vtkm::cont::DynamicCellSet inCellSet =
inData.GetCellSet(static_cast<vtkm::IdComponent>(cellSetIndex));
@ -54,6 +62,16 @@ inline VTKM_CONT vtkm::cont::DataSet CleanGrid::DoExecute(const vtkm::cont::Data
outputCellSets[cellSetIndex]);
}
VecId numCoordSystems = static_cast<VecId>(inData.GetNumberOfCoordinateSystems());
std::vector<vtkm::cont::CoordinateSystem> outputCoordinateSystems(numCoordSystems);
// Start with a shallow copy of the coordinate systems
for (VecId coordSystemIndex = 0; coordSystemIndex < numCoordSystems; ++coordSystemIndex)
{
outputCoordinateSystems[coordSystemIndex] =
inData.GetCoordinateSystem(static_cast<vtkm::IdComponent>(coordSystemIndex));
}
// Optionally adjust the cell set indices to remove all unused points
if (this->GetCompactPointFields())
{
@ -64,10 +82,56 @@ inline VTKM_CONT vtkm::cont::DataSet CleanGrid::DoExecute(const vtkm::cont::Data
}
this->PointCompactor.FindPointsEnd();
for (VecId cellSetIndex = 0; cellSetIndex < numCellSets; cellSetIndex++)
for (VecId cellSetIndex = 0; cellSetIndex < numCellSets; ++cellSetIndex)
{
outputCellSets[cellSetIndex] = this->PointCompactor.MapCellSet(outputCellSets[cellSetIndex]);
}
for (VecId coordSystemIndex = 0; coordSystemIndex < numCoordSystems; ++coordSystemIndex)
{
outputCoordinateSystems[coordSystemIndex] =
vtkm::cont::CoordinateSystem(outputCoordinateSystems[coordSystemIndex].GetName(),
this->PointCompactor.MapPointFieldDeep(
outputCoordinateSystems[coordSystemIndex].GetData()));
}
}
// Optionally find and merge coincident points
if (this->GetMergePoints())
{
vtkm::cont::CoordinateSystem activeCoordSystem =
outputCoordinateSystems[static_cast<VecId>(this->GetActiveCoordinateSystemIndex())];
vtkm::Bounds bounds = activeCoordSystem.GetBounds();
vtkm::Float64 delta = this->GetTolerance();
if (!this->GetToleranceIsAbsolute())
{
delta *=
vtkm::Magnitude(vtkm::make_Vec(bounds.X.Length(), bounds.Y.Length(), bounds.Z.Length()));
}
auto coordArray = activeCoordSystem.GetData();
this->PointMerger.Run(delta, this->GetFastMerge(), bounds, coordArray);
activeCoordSystem = vtkm::cont::CoordinateSystem(activeCoordSystem.GetName(), coordArray);
for (VecId coordSystemIndex = 0; coordSystemIndex < numCoordSystems; ++coordSystemIndex)
{
if (coordSystemIndex == static_cast<VecId>(this->GetActiveCoordinateSystemIndex()))
{
outputCoordinateSystems[coordSystemIndex] = activeCoordSystem;
}
else
{
outputCoordinateSystems[coordSystemIndex] = vtkm::cont::CoordinateSystem(
outputCoordinateSystems[coordSystemIndex].GetName(),
this->PointMerger.MapPointField(outputCoordinateSystems[coordSystemIndex].GetData()));
}
}
for (VecId cellSetIndex = 0; cellSetIndex < numCellSets; ++cellSetIndex)
{
outputCellSets[cellSetIndex] = this->PointMerger.MapCellSet(outputCellSets[cellSetIndex]);
}
}
// Construct resulting data set with new cell sets
@ -78,27 +142,9 @@ inline VTKM_CONT vtkm::cont::DataSet CleanGrid::DoExecute(const vtkm::cont::Data
}
// Pass the coordinate systems
// TODO: This is very awkward. First of all, there is no support for dealing
// with coordinate systems at all. That is fine if you are computing a new
// coordinate system, but a pain if you are deriving the coordinate system
// array. Second, why is it that coordinate systems are automatically mapped
// but other fields are not? Why shouldn't the Execute of a filter also set
// up all the fields of the output data set?
for (vtkm::IdComponent coordSystemIndex = 0;
coordSystemIndex < inData.GetNumberOfCoordinateSystems();
coordSystemIndex++)
for (VecId coordSystemIndex = 0; coordSystemIndex < numCoordSystems; ++coordSystemIndex)
{
vtkm::cont::CoordinateSystem coordSystem = inData.GetCoordinateSystem(coordSystemIndex);
if (this->GetCompactPointFields())
{
auto outArray = this->MapPointField(coordSystem.GetData());
outData.AddCoordinateSystem(vtkm::cont::CoordinateSystem(coordSystem.GetName(), outArray));
}
else
{
outData.AddCoordinateSystem(coordSystem);
}
outData.AddCoordinateSystem(outputCoordinateSystems[coordSystemIndex]);
}
return outData;
@ -111,9 +157,21 @@ inline VTKM_CONT bool CleanGrid::DoMapField(
const vtkm::filter::FieldMetadata& fieldMeta,
vtkm::filter::PolicyBase<Policy>)
{
if (this->GetCompactPointFields() && fieldMeta.IsPointField())
if (fieldMeta.IsPointField() && (this->GetCompactPointFields() || this->GetMergePoints()))
{
vtkm::cont::ArrayHandle<ValueType> compactedArray = this->MapPointField(input);
vtkm::cont::ArrayHandle<ValueType> compactedArray;
if (this->GetCompactPointFields())
{
compactedArray = this->PointCompactor.MapPointFieldDeep(input);
if (this->GetMergePoints())
{
compactedArray = this->PointMerger.MapPointField(compactedArray);
}
}
else if (this->GetMergePoints())
{
compactedArray = this->PointMerger.MapPointField(input);
}
result.AddField(fieldMeta.AsField(compactedArray));
}
else
@ -123,14 +181,7 @@ inline VTKM_CONT bool CleanGrid::DoMapField(
return true;
}
}
}
template <typename ValueType, typename Storage>
inline VTKM_CONT vtkm::cont::ArrayHandle<ValueType> CleanGrid::MapPointField(
const vtkm::cont::ArrayHandle<ValueType, Storage>& inArray) const
{
VTKM_ASSERT(this->GetCompactPointFields());
return this->PointCompactor.MapPointFieldDeep(inArray);
}
}
}
#endif //vtk_m_filter_CleanGrid_hxx

1
vtkm/filter/ExternalFaces.hxx
View File

@ -100,6 +100,7 @@ inline VTKM_CONT vtkm::cont::DataSet ExternalFaces::DoExecute(
if (this->CompactPoints)
{
this->Compactor.SetCompactPointFields(true);
this->Compactor.SetMergePoints(false);
return this->Compactor.DoExecute(output, GetCellSetExplicitPolicy(policy));
}
else

1
vtkm/filter/ExtractPoints.hxx
View File

@ -80,6 +80,7 @@ inline vtkm::cont::DataSet ExtractPoints::DoExecute(const vtkm::cont::DataSet& i
if (this->CompactPoints)
{
this->Compactor.SetCompactPointFields(true);
this->Compactor.SetMergePoints(false);
return this->Compactor.DoExecute(output, GetCellSetSingleTypePolicy(policy));
}
else

1
vtkm/filter/Mask.h
View File

@ -65,7 +65,6 @@ public:
private:
vtkm::Id Stride;
bool CompactPoints;
vtkm::filter::CleanGrid Compactor;
vtkm::worklet::Mask Worklet;
};
}

1
vtkm/filter/MaskPoints.hxx
View File

@ -73,6 +73,7 @@ inline VTKM_CONT vtkm::cont::DataSet MaskPoints::DoExecute(
if (this->CompactPoints)
{
this->Compactor.SetCompactPointFields(true);
this->Compactor.SetMergePoints(false);
return this->Compactor.DoExecute(output, GetCellSetSingleTypePolicy(policy));
}
else

1
vtkm/filter/ThresholdPoints.hxx
View File

@ -198,6 +198,7 @@ inline VTKM_CONT vtkm::cont::DataSet ThresholdPoints::DoExecute(
if (this->CompactPoints)
{
this->Compactor.SetCompactPointFields(true);
this->Compactor.SetMergePoints(true);
return this->Compactor.DoExecute(output, GetCellSetSingleTypePolicy(policy));
}
else

124
vtkm/filter/testing/UnitTestCleanGrid.cxx
View File

@ -20,6 +20,8 @@
#include <vtkm/filter/CleanGrid.h>
#include <vtkm/filter/MarchingCubes.h>
#include <vtkm/cont/testing/MakeTestDataSet.h>
#include <vtkm/cont/testing/Testing.h>
@ -37,46 +39,144 @@ void TestUniformGrid(vtkm::filter::CleanGrid clean)
clean.SetFieldsToPass({ "pointvar", "cellvar" });
vtkm::cont::DataSet outData = clean.Execute(inData);
VTKM_TEST_ASSERT(outData.HasField("pointvar"), "Failed to map point field");
VTKM_TEST_ASSERT(outData.HasField("cellvar"), "Failed to map point field");
VTKM_TEST_ASSERT(outData.HasField("cellvar"), "Failed to map cell field");
vtkm::cont::CellSetExplicit<> outCellSet;
outData.GetCellSet().CopyTo(outCellSet);
VTKM_TEST_ASSERT(outCellSet.GetNumberOfPoints() == 6, "Wrong number of points");
VTKM_TEST_ASSERT(outCellSet.GetNumberOfCells() == 2, "Wrong number of cells");
VTKM_TEST_ASSERT(outCellSet.GetNumberOfPoints() == 6,
"Wrong number of points: ",
outCellSet.GetNumberOfPoints());
VTKM_TEST_ASSERT(
outCellSet.GetNumberOfCells() == 2, "Wrong number of cells: ", outCellSet.GetNumberOfCells());
vtkm::Vec<vtkm::Id, 4> cellIds;
outCellSet.GetIndices(0, cellIds);
VTKM_TEST_ASSERT((cellIds == vtkm::Vec<vtkm::Id, 4>(0, 1, 4, 3)), "Bad cell ids");
VTKM_TEST_ASSERT((cellIds == vtkm::Vec<vtkm::Id, 4>(0, 1, 4, 3)), "Bad cell ids: ", cellIds);
outCellSet.GetIndices(1, cellIds);
VTKM_TEST_ASSERT((cellIds == vtkm::Vec<vtkm::Id, 4>(1, 2, 5, 4)), "Bad cell ids");
VTKM_TEST_ASSERT((cellIds == vtkm::Vec<vtkm::Id, 4>(1, 2, 5, 4)), "Bad cell ids: ", cellIds);
vtkm::cont::ArrayHandle<vtkm::Float32> outPointField;
outData.GetField("pointvar").GetData().CopyTo(outPointField);
VTKM_TEST_ASSERT(outPointField.GetNumberOfValues() == 6, "Wrong point field size.");
VTKM_TEST_ASSERT(outPointField.GetNumberOfValues() == 6,
"Wrong point field size: ",
outPointField.GetNumberOfValues());
VTKM_TEST_ASSERT(test_equal(outPointField.GetPortalConstControl().Get(1), 20.1),
"Bad point field value");
"Bad point field value: ",
outPointField.GetPortalConstControl().Get(1));
VTKM_TEST_ASSERT(test_equal(outPointField.GetPortalConstControl().Get(4), 50.1),
"Bad point field value");
"Bad point field value: ",
outPointField.GetPortalConstControl().Get(1));
vtkm::cont::ArrayHandle<vtkm::Float32> outCellField;
outData.GetField("cellvar").GetData().CopyTo(outCellField);
VTKM_TEST_ASSERT(outCellField.GetNumberOfValues() == 2, "Wrong cell field size.");
VTKM_TEST_ASSERT(test_equal(outCellField.GetPortalConstControl().Get(0), 100.1),
"Bad cell field value");
"Bad cell field value",
outCellField.GetPortalConstControl().Get(0));
VTKM_TEST_ASSERT(test_equal(outCellField.GetPortalConstControl().Get(1), 200.1),
"Bad cell field value");
"Bad cell field value",
outCellField.GetPortalConstControl().Get(0));
}
void TestPointMerging()
{
vtkm::cont::testing::MakeTestDataSet makeDataSet;
vtkm::cont::DataSet baseData = makeDataSet.Make3DUniformDataSet3(vtkm::Id3(4, 4, 4));
vtkm::filter::MarchingCubes marchingCubes;
marchingCubes.SetIsoValue(0.05);
marchingCubes.SetMergeDuplicatePoints(false);
marchingCubes.SetActiveField("pointvar");
vtkm::cont::DataSet inData = marchingCubes.Execute(baseData);
constexpr vtkm::Id originalNumPoints = 228;
constexpr vtkm::Id originalNumCells = 76;
VTKM_TEST_ASSERT(inData.GetCellSet().GetNumberOfPoints() == originalNumPoints);
VTKM_TEST_ASSERT(inData.GetCellSet().GetNumberOfCells() == originalNumCells);
vtkm::filter::CleanGrid cleanGrid;
std::cout << "Clean grid without any merging" << std::endl;
cleanGrid.SetCompactPointFields(false);
cleanGrid.SetMergePoints(false);
vtkm::cont::DataSet noMerging = cleanGrid.Execute(inData);
VTKM_TEST_ASSERT(noMerging.GetCellSet().GetNumberOfCells() == originalNumCells);
VTKM_TEST_ASSERT(noMerging.GetCellSet().GetNumberOfPoints() == originalNumPoints);
VTKM_TEST_ASSERT(noMerging.GetCoordinateSystem().GetData().GetNumberOfValues() ==
originalNumPoints);
VTKM_TEST_ASSERT(noMerging.GetField("pointvar").GetData().GetNumberOfValues() ==
originalNumPoints);
std::cout << "Clean grid by merging very close points" << std::endl;
cleanGrid.SetMergePoints(true);
cleanGrid.SetFastMerge(false);
vtkm::cont::DataSet closeMerge = cleanGrid.Execute(inData);
constexpr vtkm::Id closeMergeNumPoints = 62;
VTKM_TEST_ASSERT(closeMerge.GetCellSet().GetNumberOfCells() == originalNumCells);
VTKM_TEST_ASSERT(closeMerge.GetCellSet().GetNumberOfPoints() == closeMergeNumPoints);
VTKM_TEST_ASSERT(closeMerge.GetCoordinateSystem().GetData().GetNumberOfValues() ==
closeMergeNumPoints);
VTKM_TEST_ASSERT(closeMerge.GetField("pointvar").GetData().GetNumberOfValues() ==
closeMergeNumPoints);
std::cout << "Clean grid by merging very close points with fast merge" << std::endl;
cleanGrid.SetFastMerge(true);
vtkm::cont::DataSet closeFastMerge = cleanGrid.Execute(inData);
VTKM_TEST_ASSERT(closeFastMerge.GetCellSet().GetNumberOfCells() == originalNumCells);
VTKM_TEST_ASSERT(closeFastMerge.GetCellSet().GetNumberOfPoints() == closeMergeNumPoints);
VTKM_TEST_ASSERT(closeFastMerge.GetCoordinateSystem().GetData().GetNumberOfValues() ==
closeMergeNumPoints);
VTKM_TEST_ASSERT(closeFastMerge.GetField("pointvar").GetData().GetNumberOfValues() ==
closeMergeNumPoints);
std::cout << "Clean grid with largely separated points" << std::endl;
cleanGrid.SetFastMerge(false);
cleanGrid.SetTolerance(0.1);
vtkm::cont::DataSet farMerge = cleanGrid.Execute(inData);
constexpr vtkm::Id farMergeNumPoints = 36;
VTKM_TEST_ASSERT(farMerge.GetCellSet().GetNumberOfCells() == originalNumCells);
VTKM_TEST_ASSERT(farMerge.GetCellSet().GetNumberOfPoints() == farMergeNumPoints);
VTKM_TEST_ASSERT(farMerge.GetCoordinateSystem().GetData().GetNumberOfValues() ==
farMergeNumPoints);
VTKM_TEST_ASSERT(farMerge.GetField("pointvar").GetData().GetNumberOfValues() ==
farMergeNumPoints);
std::cout << "Clean grid with largerly separated points quickly" << std::endl;
cleanGrid.SetFastMerge(true);
vtkm::cont::DataSet farFastMerge = cleanGrid.Execute(inData);
constexpr vtkm::Id farFastMergeNumPoints = 19;
VTKM_TEST_ASSERT(farFastMerge.GetCellSet().GetNumberOfCells() == originalNumCells);
VTKM_TEST_ASSERT(farFastMerge.GetCellSet().GetNumberOfPoints() == farFastMergeNumPoints);
VTKM_TEST_ASSERT(farFastMerge.GetCoordinateSystem().GetData().GetNumberOfValues() ==
farFastMergeNumPoints);
VTKM_TEST_ASSERT(farFastMerge.GetField("pointvar").GetData().GetNumberOfValues() ==
farFastMergeNumPoints);
}
void RunTest()
{
vtkm::filter::CleanGrid clean;
std::cout << "*** Test wqith compact point fields on" << std::endl;
std::cout << "*** Test with compact point fields on merge points off" << std::endl;
clean.SetCompactPointFields(true);
clean.SetMergePoints(false);
TestUniformGrid(clean);
std::cout << "*** Test wqith compact point fields off" << std::endl;
std::cout << "*** Test with compact point fields off merge points off" << std::endl;
clean.SetCompactPointFields(false);
clean.SetMergePoints(false);
TestUniformGrid(clean);
std::cout << "*** Test with compact point fields on merge points on" << std::endl;
clean.SetCompactPointFields(true);
clean.SetMergePoints(true);
TestUniformGrid(clean);
std::cout << "*** Test with compact point fields off merge points on" << std::endl;
clean.SetCompactPointFields(false);
clean.SetMergePoints(true);
TestUniformGrid(clean);
std::cout << "*** Test point merging" << std::endl;
TestPointMerging();
}
} // anonymous namespace

2
vtkm/filter/testing/UnitTestExternalFacesFilter.cxx
View File

@ -35,6 +35,8 @@ vtkm::cont::DataSet MakeDataTestSet1()
vtkm::cont::DataSet ds = MakeTestDataSet().Make3DUniformDataSet1();
vtkm::filter::CleanGrid clean;
clean.SetCompactPointFields(false);
clean.SetMergePoints(false);
return clean.Execute(ds);
}

1
vtkm/filter/testing/UnitTestMarchingCubesFilter.cxx
View File

@ -503,6 +503,7 @@ void TestMarchingCubesNormals()
std::cout << "\tUnstructured dataset\n";
vtkm::filter::CleanGrid makeUnstructured;
makeUnstructured.SetCompactPointFields(false);
makeUnstructured.SetMergePoints(false);
makeUnstructured.SetFieldsToPass("pointvar");
auto result = makeUnstructured.Execute(dataset);
TestNormals(result, false);

1
vtkm/filter/testing/UnitTestZFP.cxx
View File

@ -25,7 +25,6 @@
#include <vtkm/cont/DataSetBuilderUniform.h>
#include <vtkm/cont/DataSetFieldAdd.h>
#include <vtkm/cont/testing/Testing.h>
#include <vtkm/filter/CleanGrid.h>
#include <vtkm/filter/ZFPCompressor1D.h>
#include <vtkm/filter/ZFPCompressor2D.h>

2
vtkm/internal/ArrayPortalValueReference.h
View File

@ -63,7 +63,7 @@ struct ArrayPortalValueReference
VTKM_SUPPRESS_EXEC_WARNINGS
VTKM_EXEC_CONT
operator ValueType(void) const { return this->Portal.Get(this->Index); }
operator ValueType(void) const { return this->Get(); }
// Declaring Set as const seems a little weird because we are changing the value. But remember
// that ArrayPortalReference is only a reference class. The reference itself does not change,

1
vtkm/worklet/CMakeLists.txt
View File

@ -62,6 +62,7 @@ set(headers
ParticleAdvection.h
PointAverage.h
PointElevation.h
PointMerge.h
PointTransform.h
Probe.h
RemoveUnusedPoints.h

502
vtkm/worklet/PointMerge.h Normal file
View File

@ -0,0 +1,502 @@
//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//
// Copyright 2014 National Technology & Engineering Solutions of Sandia, LLC (NTESS).
// Copyright 2014 UT-Battelle, LLC.
// Copyright 2014 Los Alamos National Security.
//
// Under the terms of Contract DE-NA0003525 with NTESS,
// the U.S. Government retains certain rights in this software.
//
// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
#ifndef vtk_m_worklet_PointMerge_h
#define vtk_m_worklet_PointMerge_h
#include <vtkm/worklet/AverageByKey.h>
#include <vtkm/worklet/DispatcherMapField.h>
#include <vtkm/worklet/DispatcherReduceByKey.h>
#include <vtkm/worklet/Invoker.h>
#include <vtkm/worklet/Keys.h>
#include <vtkm/worklet/RemoveUnusedPoints.h>
#include <vtkm/worklet/WorkletMapField.h>
#include <vtkm/worklet/WorkletReduceByKey.h>
#include <vtkm/cont/ArrayCopy.h>
#include <vtkm/cont/ArrayHandleIndex.h>
#include <vtkm/cont/ArrayHandlePermutation.h>
#include <vtkm/cont/ArrayHandleVirtual.h>
#include <vtkm/cont/CellSetExplicit.h>
#include <vtkm/cont/ExecutionAndControlObjectBase.h>
#include <vtkm/Bounds.h>
#include <vtkm/Hash.h>
#include <vtkm/Math.h>
#include <vtkm/VectorAnalysis.h>
namespace vtkm
{
namespace worklet
{
class PointMerge
{
public:
// This class can take point worldCoords as inputs and return the bin index of the enclosing bin.
class BinLocator : public vtkm::cont::ExecutionAndControlObjectBase
{
vtkm::Vec<vtkm::Float64, 3> Offset;
vtkm::Vec<vtkm::Float64, 3> Scale;
#ifdef VTKM_USE_64BIT_IDS
// IEEE double precision floating point as 53 bits for the significand, so it would not be
// possible to represent a number with more precision than that. We also back off a few bits to
// avoid potential issues with numerical imprecision in the scaling.
static constexpr vtkm::IdComponent BitsPerDimension = 50;
#else
static constexpr vtkm::IdComponent BitsPerDimension = 31;
#endif
static constexpr vtkm::Id MaxBinsPerDimension =
static_cast<vtkm::Id>((1LL << BitsPerDimension) - 1);
public:
VTKM_CONT BinLocator()
: Offset(0.0)
, Scale(0.0)
{
}
VTKM_CONT
static vtkm::Vec<vtkm::Float64, 3> ComputeBinWidths(const vtkm::Bounds& bounds,
vtkm::Float64 delta)
{
const vtkm::Vec<vtkm::Float64, 3> boundLengths(
bounds.X.Length() + delta, bounds.Y.Length() + delta, bounds.Z.Length() + delta);
vtkm::Vec<vtkm::Float64, 3> binWidths;
for (vtkm::IdComponent dimIndex = 0; dimIndex < 3; ++dimIndex)
{
if (boundLengths[dimIndex] > vtkm::Epsilon64())
{
vtkm::Float64 minBinWidth = boundLengths[dimIndex] / (MaxBinsPerDimension - 1);
if (minBinWidth < (2 * delta))
{
// We can accurately represent delta with the precision of the bin indices. The bin
// size is 2*delta, which means we scale the (offset) point coordinates by 1/delta to
// get the bin index.
binWidths[dimIndex] = 2.0 * delta;
}
else
{
// Scale the (offset) point coordinates by 1/minBinWidth, which will give us bin
// indices between 0 and MaxBinsPerDimension - 1.
binWidths[dimIndex] = minBinWidth;
}
}
else
{
// Bounds are essentially 0 in this dimension. The scale does not matter so much.
binWidths[dimIndex] = 1.0;
}
}
return binWidths;
}
// Constructs a BinLocator such that all bins are at least 2*delta large. The bins might be
// made larger than that if there would be too many bins for the precision of vtkm::Id.
VTKM_CONT
BinLocator(const vtkm::Bounds& bounds, vtkm::Float64 delta = 0.0)
: Offset(bounds.X.Min, bounds.Y.Min, bounds.Z.Min)
{
const vtkm::Vec<vtkm::Float64, 3> binWidths = ComputeBinWidths(bounds, delta);
this->Scale = vtkm::Vec<vtkm::Float64, 3>(1.0) / binWidths;
}
// Shifts the grid by delta in the specified directions. This will allow the bins to cover
// neighbors that straddled the boundaries of the original.
VTKM_CONT
BinLocator ShiftBins(const vtkm::Bounds& bounds,
vtkm::Float64 delta,
const vtkm::Vec<bool, 3>& directions)
{
const vtkm::Vec<vtkm::Float64, 3> binWidths = ComputeBinWidths(bounds, delta);
BinLocator shiftedLocator(*this);
for (vtkm::IdComponent dimIndex = 0; dimIndex < 3; ++dimIndex)
{
if (directions[dimIndex])
{
shiftedLocator.Offset[dimIndex] -= (0.5 * binWidths[dimIndex]);
}
}
return shiftedLocator;
}
template <typename T>
VTKM_EXEC_CONT vtkm::Id3 FindBin(const vtkm::Vec<T, 3>& worldCoords) const
{
vtkm::Vec<vtkm::Float64, 3> relativeCoords = (worldCoords - this->Offset) * this->Scale;
return vtkm::Id3(vtkm::Floor(relativeCoords));
}
// Because this class is a POD, we can reuse it in both control and execution environments.
template <typename Device>
BinLocator PrepareForExecution(Device) const
{
return *this;
}
BinLocator PrepareForControl() const { return *this; }
};
// Converts point coordinates to a hash that represents the bin.
struct CoordsToHash : public vtkm::worklet::WorkletMapField
{
using ControlSignature = void(FieldIn pointCoordinates,
ExecObject binLocator,
FieldOut hashesOut);
using ExecutionSignature = void(_1, _2, _3);
template <typename T>
VTKM_EXEC void operator()(const vtkm::Vec<T, 3>& coordiantes,
const BinLocator binLocator,
vtkm::HashType& hashOut) const
{
vtkm::Id3 binId = binLocator.FindBin(coordiantes);
hashOut = vtkm::Hash(binId);
}
};
class FindNeighbors : public vtkm::worklet::WorkletReduceByKey
{
vtkm::Float64 DeltaSquared;
bool FastCheck;
public:
VTKM_CONT
FindNeighbors(bool fastCheck = true, vtkm::Float64 delta = vtkm::Epsilon64())
: DeltaSquared(delta * delta)
, FastCheck(fastCheck)
{
}
using ControlSignature = void(KeysIn keys,
ValuesInOut pointIndices,
ValuesInOut pointCoordinates,
ExecObject binLocator,
ValuesOut neighborIndices);
using ExecutionSignature = void(_2, _3, _4, _5);
template <typename IndexVecInType, typename CoordinateVecInType, typename IndexVecOutType>
VTKM_EXEC void operator()(IndexVecInType& pointIndices,
CoordinateVecInType& pointCoordinates,
const BinLocator& binLocator,
IndexVecOutType& neighborIndices) const
{
// For each point we are going to find all points close enough to be considered neighbors. We
// record the neighbors by filling in the same index into neighborIndices. That is, if two
// items in neighborIndices have the same value, they should be considered neighbors.
// Otherwise, they should not. We will use the "local" index, which refers to index in the
// vec-like objects passed into this worklet. This allows us to quickly identify the local
// point without sorting through the global indices.
using CoordType = typename CoordinateVecInType::ComponentType;
vtkm::IdComponent numPoints = pointIndices.GetNumberOfComponents();
VTKM_ASSERT(numPoints == pointCoordinates.GetNumberOfComponents());
VTKM_ASSERT(numPoints == neighborIndices.GetNumberOfComponents());
// Initially, set every point to be its own neighbor.
for (vtkm::IdComponent i = 0; i < numPoints; ++i)
{
neighborIndices[i] = i;
}
// Iterate over every point and look for neighbors. Only need to look to numPoints-1 since we
// only need to check points after the current index (earlier points are already checked).
for (vtkm::IdComponent i = 0; i < (numPoints - 1); ++i)
{
CoordType p0 = pointCoordinates[i];
vtkm::Id3 bin0 = binLocator.FindBin(p0);
// Check all points after this one. (All those before already checked themselves to this.)
for (vtkm::IdComponent j = i + 1; j < numPoints; ++j)
{
if (neighborIndices[i] == neighborIndices[j])
{
// We have already identified these points as neighbors. Can skip the check.
continue;
}
CoordType p1 = pointCoordinates[j];
vtkm::Id3 bin1 = binLocator.FindBin(p1);
// Check to see if these points should be considered neighbors. First, check to make sure
// that they are in the same bin. If they are not, then they cannot be neighbors. Next,
// check the FastCheck flag. If fast checking is on, then all points in the same bin are
// considered neighbors. Otherwise, check that the distance is within the specified
// delta. If so, mark them as neighbors.
if ((bin0 == bin1) &&
(this->FastCheck || (this->DeltaSquared >= vtkm::MagnitudeSquared(p0 - p1))))
{
// The two points should be merged. But we also might need to merge larger
// neighborhoods.
if (neighborIndices[j] == j)
{
// Second point not yet merged into another neighborhood. We can just take it.
neighborIndices[j] = neighborIndices[i];
}
else
{
// The second point is already part of a neighborhood. Merge the neighborhood with
// the largest index into the neighborhood with the smaller index.
vtkm::IdComponent neighborhoodToGrow;
vtkm::IdComponent neighborhoodToAbsorb;
if (neighborIndices[i] < neighborIndices[j])
{
neighborhoodToGrow = neighborIndices[i];
neighborhoodToAbsorb = neighborIndices[j];
}
else
{
neighborhoodToGrow = neighborIndices[j];
neighborhoodToAbsorb = neighborIndices[i];
}
// Change all neighborhoodToAbsorb indices to neighborhoodToGrow.
for (vtkm::IdComponent k = neighborhoodToAbsorb; k < numPoints; ++k)
{
if (neighborIndices[k] == neighborhoodToAbsorb)
{
neighborIndices[k] = neighborhoodToGrow;
}
}
}
} // if merge points
} // for each p1
} // for each p0
// We have finished grouping neighbors. neighborIndices contains a unique local index for
// each neighbor group. Now find the average (centroid) point coordinates for each group and
// write those coordinates back into the coordinates array. Also modify the point indices
// so that all indices of a group are the same. (This forms a map from old point indices to
// merged point indices.)
for (vtkm::IdComponent i = 0; i < numPoints; ++i)
{
vtkm::IdComponent neighborhood = neighborIndices[i];
if (i == neighborhood)
{
// Found a new group. Find the centroid.
CoordType centroid = pointCoordinates[i];
vtkm::IdComponent numInGroup = 1;
for (vtkm::IdComponent j = i + 1; j < numPoints; ++j)
{
if (neighborhood == neighborIndices[j])
{
centroid = centroid + pointCoordinates[j];
++numInGroup;
}
}
centroid = centroid / numInGroup;
// Now that we have the centroid, write new point coordinates and index.
vtkm::Id groupIndex = pointIndices[i];
pointCoordinates[i] = centroid;
for (vtkm::IdComponent j = i + 1; j < numPoints; ++j)
{
if (neighborhood == neighborIndices[j])
{
pointCoordinates[j] = centroid;
pointIndices[j] = groupIndex;
}
}
}
}
}
};
struct BuildPointInputToOutputMap : vtkm::worklet::WorkletReduceByKey
{
using ControlSignature = void(KeysIn, ValuesOut PointInputToOutputMap);
using ExecutionSignature = void(InputIndex, _2);
template <typename MapPortalType>
VTKM_EXEC void operator()(vtkm::Id newIndex, MapPortalType outputIndices) const
{
const vtkm::IdComponent numIndices = outputIndices.GetNumberOfComponents();
for (vtkm::IdComponent i = 0; i < numIndices; ++i)
{
outputIndices[i] = newIndex;
}
}
};
private:
template <typename T>
VTKM_CONT static void RunOneIteration(
vtkm::Float64 delta, // Distance to consider two points coincident
bool fastCheck, // If true, approximate distances are used
const BinLocator& binLocator, // Used to find nearby points
vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>>& points, // coordinates, modified to merge close
vtkm::cont::ArrayHandle<vtkm::Id> indexNeighborMap) // identifies each neighbor group, updated
{
Invoker invoker;
vtkm::cont::ArrayHandle<vtkm::HashType> hashes;
invoker(CoordsToHash(), points, binLocator, hashes);
vtkm::worklet::Keys<HashType> keys(hashes);
// Really just scratch space
vtkm::cont::ArrayHandle<vtkm::IdComponent> neighborIndices;
invoker(
FindNeighbors(fastCheck, delta), keys, indexNeighborMap, points, binLocator, neighborIndices);
}
public:
template <typename T>
VTKM_CONT void Run(
vtkm::Float64 delta, // Distance to consider two points coincident
bool fastCheck, // If true, approximate distances are used
const vtkm::Bounds& bounds, // Bounds of points
vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>>& points) // coordinates, modified to merge close
{
Invoker invoker;
BinLocator binLocator(bounds, delta);
vtkm::cont::ArrayHandle<vtkm::Id> indexNeighborMap;
vtkm::cont::ArrayCopy(vtkm::cont::ArrayHandleIndex(points.GetNumberOfValues()),
indexNeighborMap);
this->RunOneIteration(delta, fastCheck, binLocator, points, indexNeighborMap);
if (!fastCheck)
{
// Run the algorithm again after shifting the bins to capture nearby points that straddled
// the previous bins.
this->RunOneIteration(delta,
fastCheck,
binLocator.ShiftBins(bounds, delta, vtkm::make_Vec(true, false, false)),
points,
indexNeighborMap);
this->RunOneIteration(delta,
fastCheck,
binLocator.ShiftBins(bounds, delta, vtkm::make_Vec(false, true, false)),
points,
indexNeighborMap);
this->RunOneIteration(delta,
fastCheck,
binLocator.ShiftBins(bounds, delta, vtkm::make_Vec(false, false, true)),
points,
indexNeighborMap);
this->RunOneIteration(delta,
fastCheck,
binLocator.ShiftBins(bounds, delta, vtkm::make_Vec(true, true, false)),
points,
indexNeighborMap);
this->RunOneIteration(delta,
fastCheck,
binLocator.ShiftBins(bounds, delta, vtkm::make_Vec(true, false, true)),
points,
indexNeighborMap);
this->RunOneIteration(delta,
fastCheck,
binLocator.ShiftBins(bounds, delta, vtkm::make_Vec(false, true, true)),
points,
indexNeighborMap);
this->RunOneIteration(delta,
fastCheck,
binLocator.ShiftBins(bounds, delta, vtkm::make_Vec(true, true, true)),
points,
indexNeighborMap);
}
this->MergeKeys = vtkm::worklet::Keys<vtkm::Id>(indexNeighborMap);
invoker(BuildPointInputToOutputMap(), this->MergeKeys, this->PointInputToOutputMap);
// Need to pull out the unique point coordiantes
vtkm::cont::ArrayHandle<vtkm::Vec<T, 3>> uniquePointCoordinates;
vtkm::cont::ArrayCopy(
vtkm::cont::make_ArrayHandlePermutation(this->MergeKeys.GetUniqueKeys(), points),
uniquePointCoordinates);
points = uniquePointCoordinates;
}
VTKM_CONT void Run(
vtkm::Float64 delta, // Distance to consider two points coincident
bool fastCheck, // If true, approximate distances are used
const vtkm::Bounds& bounds, // Bounds of points
vtkm::cont::ArrayHandleVirtualCoordinates& points) // coordinates, modified to merge close
{
// Get a cast to a concrete set of point coordiantes so that it can be modified in place
vtkm::cont::ArrayHandle<vtkm::Vec<vtkm::FloatDefault, 3>> concretePoints;
if (points.IsType<decltype(concretePoints)>())
{
concretePoints = points.Cast<decltype(concretePoints)>();
}
else
{
vtkm::cont::ArrayCopy(points, concretePoints);
}
Run(delta, fastCheck, bounds, concretePoints);
// Make sure that the modified points are reflected back in the virtual array.
points = vtkm::cont::ArrayHandleVirtualCoordinates(concretePoints);
}
template <typename ShapeStorage,
typename NumIndicesStorage,
typename ConnectivityStorage,
typename OffsetsStorage>
VTKM_CONT vtkm::cont::CellSetExplicit<ShapeStorage,
NumIndicesStorage,
VTKM_DEFAULT_CONNECTIVITY_STORAGE_TAG,
OffsetsStorage>
MapCellSet(const vtkm::cont::CellSetExplicit<ShapeStorage,
NumIndicesStorage,
ConnectivityStorage,
OffsetsStorage>& inCellSet) const
{
return vtkm::worklet::RemoveUnusedPoints::MapCellSet(
inCellSet, this->PointInputToOutputMap, this->MergeKeys.GetInputRange());
}
template <typename InArrayHandle, typename OutArrayHandle>
VTKM_CONT void MapPointField(const InArrayHandle& inArray, OutArrayHandle& outArray) const
{
VTKM_IS_ARRAY_HANDLE(InArrayHandle);
VTKM_IS_ARRAY_HANDLE(OutArrayHandle);
vtkm::worklet::AverageByKey::Run(this->MergeKeys, inArray, outArray);
}
template <typename InArrayHandle>
VTKM_CONT vtkm::cont::ArrayHandle<typename InArrayHandle::ValueType> MapPointField(
const InArrayHandle& inArray) const
{
VTKM_IS_ARRAY_HANDLE(InArrayHandle);
vtkm::cont::ArrayHandle<typename InArrayHandle::ValueType> outArray;
this->MapPointField(inArray, outArray);
return outArray;
}
private:
vtkm::worklet::Keys<vtkm::Id> MergeKeys;
vtkm::cont::ArrayHandle<vtkm::Id> PointInputToOutputMap;
};
}
} // namespace vtkm::worklet
#endif //vtk_m_worklet_PointMerge_h

42
vtkm/worklet/RemoveUnusedPoints.h
View File

@ -141,7 +141,7 @@ public:
/// \brief Map cell indices
///
/// Given a cell set (typically the same one passed to the constructor) and
/// Given a cell set (typically the same one passed to the constructor)
/// returns a new cell set with cell points transformed to use the indices of
/// the new reduced point arrays.
///
@ -157,22 +157,54 @@ public:
NumIndicesStorage,
ConnectivityStorage,
OffsetsStorage>& inCellSet) const
{
VTKM_ASSERT(this->PointScatter);
return MapCellSet(inCellSet,
this->PointScatter->GetInputToOutputMap(),
this->PointScatter->GetOutputToInputMap().GetNumberOfValues());
}
/// \brief Map cell indices
///
/// Given a cell set (typically the same one passed to the constructor) and
/// an array that maps point indices from an old set of indices to a new set,
/// returns a new cell set with cell points transformed to use the indices of
/// the new reduced point arrays.
///
/// This helper method can be used by external items that do similar operations
/// that remove points or otherwise rearange points in a cell set. If points
/// were removed by calling \c FindPoints, then you should use the other form
/// of \c MapCellSet.
///
template <typename ShapeStorage,
typename NumIndicesStorage,
typename ConnectivityStorage,
typename OffsetsStorage,
typename MapStorage>
VTKM_CONT static vtkm::cont::CellSetExplicit<ShapeStorage,
NumIndicesStorage,
VTKM_DEFAULT_CONNECTIVITY_STORAGE_TAG,
OffsetsStorage>
MapCellSet(const vtkm::cont::CellSetExplicit<ShapeStorage,
NumIndicesStorage,
ConnectivityStorage,
OffsetsStorage>& inCellSet,
const vtkm::cont::ArrayHandle<vtkm::Id, MapStorage>& inputToOutputPointMap,
vtkm::Id numberOfPoints)
{
using FromTopology = vtkm::TopologyElementTagPoint;
using ToTopology = vtkm::TopologyElementTagCell;
using NewConnectivityStorage = VTKM_DEFAULT_CONNECTIVITY_STORAGE_TAG;
VTKM_ASSERT(this->PointScatter);
vtkm::cont::ArrayHandle<vtkm::Id, NewConnectivityStorage> newConnectivityArray;
vtkm::worklet::DispatcherMapField<TransformPointIndices> dispatcher;
dispatcher.Invoke(inCellSet.GetConnectivityArray(FromTopology(), ToTopology()),
this->PointScatter->GetInputToOutputMap(),
inputToOutputPointMap,
newConnectivityArray);
vtkm::Id numberOfPoints = this->PointScatter->GetOutputToInputMap().GetNumberOfValues();
vtkm::cont::
CellSetExplicit<ShapeStorage, NumIndicesStorage, NewConnectivityStorage, OffsetsStorage>
outCellSet(inCellSet.GetName());