Add a new unstructured cell locator

This commit is contained in:
Dave Pugmire 2022-11-22 09:51:40 -05:00
parent a7fcbf4a62
commit 7d2d7fe1cc
8 changed files with 747 additions and 16 deletions

@ -149,6 +149,22 @@ struct Bounds
return vtkm::Vec3f_64(this->X.Center(), this->Y.Center(), this->Z.Center());
}
/// \b Returns the min point of the bounds
///
/// \c MinCorder returns the minium point of the bounds.If the bounds
/// are empty, the results are undefined.
///
VTKM_EXEC_CONT
vtkm::Vec3f_64 MinCorner() const { return vtkm::Vec3f_64(this->X.Min, this->Y.Min, this->Z.Min); }
/// \b Returns the max point of the bounds
///
/// \c MaxCorder returns the minium point of the bounds.If the bounds
/// are empty, the results are undefined.
///
VTKM_EXEC_CONT
vtkm::Vec3f_64 MaxCorner() const { return vtkm::Vec3f_64(this->X.Max, this->Y.Max, this->Z.Max); }
/// \b Expand bounds to include a point.
///
/// This version of \c Include expands the bounds just enough to include the

@ -62,6 +62,7 @@ set(headers
CellLocatorPartitioned.h
CellLocatorRectilinearGrid.h
CellLocatorTwoLevel.h
CellLocatorUniformBins.h
CellLocatorUniformGrid.h
CellSet.h
CellSetExplicit.h
@ -139,6 +140,7 @@ set(sources
CellLocatorGeneral.cxx
CellLocatorPartitioned.cxx
CellLocatorRectilinearGrid.cxx
CellLocatorUniformBins.cxx
CellLocatorUniformGrid.cxx
CellSet.cxx
CellSetStructured.cxx

@ -0,0 +1,340 @@
//============================================================================
// 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.
//============================================================================
#include <vtkm/cont/Algorithm.h>
#include <vtkm/cont/ArrayHandleGroupVecVariable.h>
#include <vtkm/cont/CellLocatorUniformBins.h>
#include <vtkm/cont/ConvertNumComponentsToOffsets.h>
#include <vtkm/cont/Invoker.h>
#include <vtkm/worklet/WorkletMapField.h>
#include <vtkm/worklet/WorkletMapTopology.h>
namespace vtkm
{
namespace cont
{
namespace
{
VTKM_EXEC
inline vtkm::Id ComputeFlatIndex(const vtkm::Id3& idx, const vtkm::Id3& dims)
{
return idx[0] + (dims[0] * (idx[1] + (dims[1] * idx[2])));
}
VTKM_EXEC
inline vtkm::Id3 GetFlatIndex(const vtkm::Vec3f& pt,
const vtkm::Vec3f& origin,
const vtkm::Vec3f& invSpacing,
const vtkm::Id3& maxCellIds)
{
auto temp = pt - origin;
temp = temp * invSpacing;
auto logicalIdx = vtkm::Min(vtkm::Id3(temp), maxCellIds);
return logicalIdx;
}
template <typename PointsVecType>
VTKM_EXEC inline void MinMaxIndicesForCellPoints(const PointsVecType& points,
const vtkm::Vec3f& origin,
const vtkm::Vec3f& invSpacing,
const vtkm::Id3& maxCellIds,
vtkm::Id3& minIdx,
vtkm::Id3& maxIdx)
{
auto numPoints = vtkm::VecTraits<PointsVecType>::GetNumberOfComponents(points);
vtkm::Bounds bounds;
for (vtkm::IdComponent i = 0; i < numPoints; ++i)
bounds.Include(points[i]);
//Get 8 corners of bbox.
minIdx = GetFlatIndex(bounds.MinCorner(), origin, invSpacing, maxCellIds);
maxIdx = GetFlatIndex(bounds.MaxCorner(), origin, invSpacing, maxCellIds);
}
class CountCellBins : public vtkm::worklet::WorkletVisitCellsWithPoints
{
public:
using ControlSignature = void(CellSetIn cellset, FieldInPoint coords, FieldOutCell bincount);
using ExecutionSignature = void(_2, _3);
using InputDomain = _1;
CountCellBins(const vtkm::Vec3f& origin,
const vtkm::Vec3f& invSpacing,
const vtkm::Id3& maxCellIds)
: InvSpacing(invSpacing)
, MaxCellIds(maxCellIds)
, Origin(origin)
{
}
template <typename PointsVecType>
VTKM_EXEC void operator()(const PointsVecType& points, vtkm::Id& numBins) const
{
vtkm::Id3 idx000, idx111;
MinMaxIndicesForCellPoints(
points, this->Origin, this->InvSpacing, this->MaxCellIds, idx000, idx111);
//Count the number of bins for this cell
numBins =
(idx111[0] - idx000[0] + 1) * (idx111[1] - idx000[1] + 1) * (idx111[2] - idx000[2] + 1);
}
private:
vtkm::Vec3f InvSpacing;
vtkm::Id3 MaxCellIds;
vtkm::Vec3f Origin;
};
class RecordBinsPerCell : public vtkm::worklet::WorkletVisitCellsWithPoints
{
public:
using ControlSignature = void(CellSetIn cellset,
FieldInPoint coords,
FieldInCell start,
WholeArrayInOut binsPerCell,
WholeArrayInOut cellIds,
AtomicArrayInOut cellCounts);
using ExecutionSignature = void(InputIndex, _2, _3, _4, _5, _6);
using InputDomain = _1;
RecordBinsPerCell(const vtkm::Vec3f& origin,
const vtkm::Vec3f& invSpacing,
const vtkm::Id3& dims,
const vtkm::Id3& maxCellIds)
: Dims(dims)
, InvSpacing(invSpacing)
, MaxCellIds(maxCellIds)
, Origin(origin)
{
}
template <typename PointsVecType, typename ResultArrayType, typename CellCountType>
VTKM_EXEC void operator()(const vtkm::Id& cellIdx,
const PointsVecType& points,
const vtkm::Id& start,
ResultArrayType& binsPerCell,
ResultArrayType& cellIds,
CellCountType cellCounts) const
{
vtkm::Id3 idx000, idx111;
MinMaxIndicesForCellPoints(
points, this->Origin, this->InvSpacing, this->MaxCellIds, idx000, idx111);
vtkm::Id cnt = 0;
vtkm::Id sliceStart = ComputeFlatIndex(idx000, this->Dims);
for (vtkm::Id k = idx000[2]; k <= idx111[2]; k++)
{
vtkm::Id shaftStart = sliceStart;
for (vtkm::Id j = idx000[1]; j <= idx111[1]; j++)
{
vtkm::Id flatIdx = shaftStart;
for (vtkm::Id i = idx000[0]; i <= idx111[0]; i++)
{
// set portals and increment cnt...
binsPerCell.Set(start + cnt, flatIdx);
cellIds.Set(start + cnt, cellIdx);
cellCounts.Add(flatIdx, 1);
++flatIdx;
++cnt;
}
shaftStart += this->Dims[0];
}
sliceStart += this->Dims[0] * this->Dims[1];
}
}
private:
vtkm::Id3 Dims;
vtkm::Vec3f InvSpacing;
vtkm::Id3 MaxCellIds;
vtkm::Vec3f Origin;
};
} //namespace detail
//----------------------------------------------------------------------------
/// Builds the cell locator lookup structure
///
VTKM_CONT void CellLocatorUniformBins::Build()
{
if (this->UniformDims[0] <= 0 || this->UniformDims[1] <= 0 || this->UniformDims[2] <= 0)
throw vtkm::cont::ErrorBadValue("Grid dimensions of CellLocatorUniformBins must be > 0");
VTKM_LOG_SCOPE(vtkm::cont::LogLevel::Perf, "CellLocatorUniformBins::Build");
this->MaxCellIds = (vtkm::Max(this->UniformDims - vtkm::Id3(1), vtkm::Id3(0)));
vtkm::Id totalNumBins = this->UniformDims[0] * this->UniformDims[1] * this->UniformDims[2];
vtkm::cont::Invoker invoker;
auto cellset = this->GetCellSet();
const auto& coords = this->GetCoordinates();
auto bounds = coords.GetBounds();
this->Origin = vtkm::Vec3f(bounds.MinCorner());
this->MaxPoint = vtkm::Vec3f(bounds.MaxCorner());
auto size = this->MaxPoint - this->Origin;
vtkm::Vec3f spacing(
size[0] / this->UniformDims[0], size[1] / this->UniformDims[1], size[2] / this->UniformDims[2]);
for (vtkm::IdComponent i = 0; i < 3; i++)
{
if (vtkm::Abs(spacing[i]) > 0)
this->InvSpacing[i] = 1.0f / spacing[i];
else
this->InvSpacing[i] = 0;
}
// The following example will be used in the explanation below
// Dataset with 3 cells: c0, c1, c2
// 2x2 uniform grid: b0, b1, b2, b3
// Assume that the bounding box for each cell overlaps as follows:
// c0: b0, b1, b2
// c1: b1
// c2: b2
//
// The acceleration structure is an array of cell ids that are grouped
// by the overlapping bin. This information can be represented using
// vtkm::cont::ArrayHandleGroupVecVariable
// In the example above:
// CellIds = {c0, c0,c1, c0,c2, }
// b0 b1 b2 b3
//
// The algorithm runs as follows:
// Given a point p, find the bin (b) that contains p.
// Do a point-in-cell test for each cell in bin b.
//
// Example: point p is in b=1
// vtkm::cont::ArrayHandleGroupVecVariable provides the offset and number
// cells in bin b=1. The offset is 1 and the count is 2.
// Do a point-in-cell test on the 2 cells that start at offset 1
// CellIds[ 1 + 0], which is c0
// CellIds[ 1 + 1], which is c1
//Computing this involves several steps which are described below.
//Step 1:
// For each cell in cellSet, count the number of bins that overlap with the cell bbox.
// For the example above
// binCountsPerCell = {3, 1, 1}
// cell0 overlaps with 3 bins
// cell1 overlaps with 1 bin
// cell2 overlaps with 1 bin
vtkm::cont::ArrayHandle<vtkm::Id> binCountsPerCell;
CountCellBins countCellBins(this->Origin, this->InvSpacing, this->MaxCellIds);
invoker(countCellBins, cellset, coords, binCountsPerCell);
//2: Compute number of unique cell/bin pairs and start indices.
//Step 2:
// Given the number of bins for each cell, we can compute the offset for each cell.
// For the example above, binOffset is:
// {0, 3, 4}
// and the total number, num is 5
// c0 begins at idx=0
// c1 begins at idx=3
// c2 begins at idx=4
vtkm::cont::ArrayHandle<vtkm::Id> binOffset;
auto num = vtkm::cont::Algorithm::ScanExclusive(binCountsPerCell, binOffset);
//Step 3:
// Now that we know the start index and numbers, we can fill an array of bin ids.
// binsPerCell is the list of binIds for each cell. In the example above,
// binsPerCell = {b0,b1,b2, b2, b2}
// \ cell0 / cell1 cell2
// We can also compute the cellIds and number of cells per bin, cellCount
// cids = {c0,c0,c0, c1, c2}
// cellCount = {3, 1, 1, 0}
// These are set using RecordBinsPerCell worklet, which does the following
// for each cell
// compute cell bbox and list of overlaping bins
// for each overlaping bin
// add the bin id to binsPerCell starting at binOffset
// add the cell id to the CellIds starting at binOffset
// increment CellCount for the bin (uses an atomic for thread safety).
vtkm::cont::ArrayHandle<vtkm::Id> binsPerCell, cids, cellCount;
binsPerCell.AllocateAndFill(num, 0);
cids.Allocate(num);
cellCount.AllocateAndFill(totalNumBins, 0);
RecordBinsPerCell recordBinsPerCell(
this->Origin, this->InvSpacing, this->UniformDims, this->MaxCellIds);
invoker(recordBinsPerCell, cellset, coords, binOffset, binsPerCell, cids, cellCount);
//Step 4:
// binsPerCell is the overlapping bins for each cell.
// We want to sort CellIds by the bin ID. SortByKey does this.
vtkm::cont::Algorithm::SortByKey(binsPerCell, cids);
// Convert the cell counts to offsets using the helper function
// vtkm::cont::ConvertNumComponentsToOffsets, and create the
// CellIds that are used as the acceleration structure.
this->CellIds = vtkm::cont::make_ArrayHandleGroupVecVariable(
cids, vtkm::cont::ConvertNumComponentsToOffsets(cellCount));
}
//----------------------------------------------------------------------------
struct CellLocatorUniformBins::MakeExecObject
{
template <typename CellSetType>
VTKM_CONT void operator()(const CellSetType& cellSet,
vtkm::cont::DeviceAdapterId device,
vtkm::cont::Token& token,
const CellLocatorUniformBins& self,
ExecObjType& execObject) const
{
using CellStructureType = CellSetContToExec<CellSetType>;
execObject = vtkm::exec::CellLocatorUniformBins<CellStructureType>(self.UniformDims,
self.Origin,
self.MaxPoint,
self.InvSpacing,
self.MaxCellIds,
self.CellIds,
cellSet,
self.GetCoordinates(),
device,
token);
}
};
CellLocatorUniformBins::ExecObjType CellLocatorUniformBins::PrepareForExecution(
vtkm::cont::DeviceAdapterId device,
vtkm::cont::Token& token) const
{
this->Update();
ExecObjType execObject;
vtkm::cont::CastAndCall(this->GetCellSet(), MakeExecObject{}, device, token, *this, execObject);
return execObject;
}
//----------------------------------------------------------------------------
void CellLocatorUniformBins::PrintSummary(std::ostream& out) const
{
out << std::endl;
out << "CellLocatorUniformBins" << std::endl;
out << " UniformDims: " << this->UniformDims << std::endl;
out << " Origin: " << this->Origin << std::endl;
out << " MaxPoint: " << this->MaxPoint << std::endl;
out << " InvSpacing: " << this->InvSpacing << std::endl;
out << " MaxCellIds: " << this->MaxCellIds << std::endl;
out << "Input CellSet: \n";
this->GetCellSet().PrintSummary(out);
out << "Input Coordinates: \n";
this->GetCoordinates().PrintSummary(out);
}
}
} // vtkm::cont

@ -0,0 +1,89 @@
//============================================================================
// 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.
//============================================================================
#ifndef vtk_m_cont_CellLocatorUniformBins_h
#define vtk_m_cont_CellLocatorUniformBins_h
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/CellSetList.h>
#include <vtkm/cont/ArrayHandleGroupVecVariable.h>
#include <vtkm/cont/internal/CellLocatorBase.h>
#include <vtkm/exec/CellLocatorMultiplexer.h>
#include <vtkm/exec/CellLocatorUniformBins.h>
namespace vtkm
{
namespace cont
{
/// \brief A locator that uses a uniform grid
///
/// `CellLocatorUniformBins` creates a cell search structure using a single uniform
/// grid. The size of the uniform grid is specified using the `SetDims` method.
/// In general, the `CellLocatorTwoLevel` has the better performance. However,
/// there are some cases where this is not the case. One example of this is
/// a uniformly dense triangle grid. In some cases the `CellLocatorUniformBins`
/// produces a more efficient search structure, especially for GPUs where memory
/// access patterns are critical to performance.
class VTKM_CONT_EXPORT CellLocatorUniformBins
: public vtkm::cont::internal::CellLocatorBase<CellLocatorUniformBins>
{
using Superclass = vtkm::cont::internal::CellLocatorBase<CellLocatorUniformBins>;
template <typename CellSetCont>
using CellSetContToExec =
typename CellSetCont::template ExecConnectivityType<vtkm::TopologyElementTagCell,
vtkm::TopologyElementTagPoint>;
public:
using SupportedCellSets = VTKM_DEFAULT_CELL_SET_LIST;
using CellExecObjectList = vtkm::ListTransform<SupportedCellSets, CellSetContToExec>;
using CellLocatorExecList =
vtkm::ListTransform<CellExecObjectList, vtkm::exec::CellLocatorUniformBins>;
using ExecObjType = vtkm::ListApply<CellLocatorExecList, vtkm::exec::CellLocatorMultiplexer>;
using LastCell = typename ExecObjType::LastCell;
CellLocatorUniformBins() {}
void SetDims(const vtkm::Id3& dims) { this->UniformDims = dims; }
vtkm::Id3 GetDims() const { return this->UniformDims; }
void PrintSummary(std::ostream& out) const;
public:
ExecObjType PrepareForExecution(vtkm::cont::DeviceAdapterId device,
vtkm::cont::Token& token) const;
private:
friend Superclass;
VTKM_CONT void Build();
vtkm::Vec3f InvSpacing;
vtkm::Vec3f MaxPoint;
vtkm::Vec3f Origin;
vtkm::Id3 UniformDims;
vtkm::Id3 MaxCellIds;
using CellIdArrayType = vtkm::cont::ArrayHandle<vtkm::Id>;
using CellIdOffsetArrayType = vtkm::cont::ArrayHandle<vtkm::Id>;
vtkm::cont::ArrayHandleGroupVecVariable<CellIdArrayType, CellIdOffsetArrayType> CellIds;
struct MakeExecObject;
};
}
} // vtkm::cont
#endif // vtk_m_cont_CellLocatorUniformBins_h

@ -97,8 +97,8 @@ set(unit_tests_device
UnitTestCellLocatorGeneral.cxx
UnitTestCellLocatorPartitioned.cxx
UnitTestCellLocatorRectilinearGrid.cxx
UnitTestCellLocatorTwoLevel.cxx
UnitTestCellLocatorUniformGrid.cxx
UnitTestCellLocatorUnstructured.cxx
UnitTestCellSet.cxx
UnitTestCellSetExplicit.cxx
UnitTestCellSetPermutation.cxx

@ -10,6 +10,7 @@
#include <vtkm/cont/ArrayCopy.h>
#include <vtkm/cont/CellLocatorTwoLevel.h>
#include <vtkm/cont/CellLocatorUniformBins.h>
#include <vtkm/cont/DataSetBuilderUniform.h>
#include <vtkm/cont/Invoker.h>
#include <vtkm/cont/testing/Testing.h>
@ -204,9 +205,10 @@ public:
}
};
void TestLastCell(vtkm::cont::CellLocatorTwoLevel& locator,
template <typename LocatorType>
void TestLastCell(LocatorType& locator,
vtkm::Id numPoints,
vtkm::cont::ArrayHandle<vtkm::cont::CellLocatorTwoLevel::LastCell>& lastCell,
vtkm::cont::ArrayHandle<typename LocatorType::LastCell>& lastCell,
const vtkm::cont::ArrayHandle<PointType>& points,
const vtkm::cont::ArrayHandle<vtkm::Id>& expCellIds,
const vtkm::cont::ArrayHandle<PointType>& expPCoords)
@ -231,16 +233,15 @@ void TestLastCell(vtkm::cont::CellLocatorTwoLevel& locator,
}
}
template <vtkm::IdComponent DIMENSIONS>
void TestCellLocator(const vtkm::Vec<vtkm::Id, DIMENSIONS>& dim, vtkm::Id numberOfPoints)
template <typename LocatorType, vtkm::IdComponent DIMENSIONS>
void TestCellLocator(LocatorType& locator,
const vtkm::Vec<vtkm::Id, DIMENSIONS>& dim,
vtkm::Id numberOfPoints)
{
auto ds = MakeTestDataSet(dim);
std::cout << "Testing " << DIMENSIONS << "D dataset with " << ds.GetNumberOfCells() << " cells\n";
vtkm::cont::CellLocatorTwoLevel locator;
locator.SetDensityL1(64.0f);
locator.SetDensityL2(1.0f);
locator.SetCellSet(ds.GetCellSet());
locator.SetCoordinates(ds.GetCoordinateSystem());
locator.Update();
@ -271,15 +272,15 @@ void TestCellLocator(const vtkm::Vec<vtkm::Id, DIMENSIONS>& dim, vtkm::Id number
//Test locator using lastCell
//Test it with initialized.
vtkm::cont::ArrayHandle<vtkm::cont::CellLocatorTwoLevel::LastCell> lastCell;
lastCell.AllocateAndFill(numberOfPoints, vtkm::cont::CellLocatorTwoLevel::LastCell{});
vtkm::cont::ArrayHandle<typename LocatorType::LastCell> lastCell;
lastCell.AllocateAndFill(numberOfPoints, typename LocatorType::LastCell{});
TestLastCell(locator, numberOfPoints, lastCell, points, expCellIds, pcoords);
//Call it again using the lastCell just computed to validate.
TestLastCell(locator, numberOfPoints, lastCell, points, expCellIds, pcoords);
//Test it with uninitialized array.
vtkm::cont::ArrayHandle<vtkm::cont::CellLocatorTwoLevel::LastCell> lastCell2;
vtkm::cont::ArrayHandle<typename LocatorType::LastCell> lastCell2;
lastCell2.Allocate(numberOfPoints);
TestLastCell(locator, numberOfPoints, lastCell2, points, expCellIds, pcoords);
@ -287,19 +288,35 @@ void TestCellLocator(const vtkm::Vec<vtkm::Id, DIMENSIONS>& dim, vtkm::Id number
TestLastCell(locator, numberOfPoints, lastCell2, points, expCellIds, pcoords);
}
void TestingCellLocatorTwoLevel()
void TestingCellLocatorUnstructured()
{
vtkm::UInt32 seed = static_cast<vtkm::UInt32>(std::time(nullptr));
std::cout << "Seed: " << seed << std::endl;
RandomGenerator.seed(seed);
TestCellLocator(vtkm::Id3(8), 512); // 3D dataset
TestCellLocator(vtkm::Id2(18), 512); // 2D dataset
//Test vtkm::cont::CellLocatorTwoLevel
vtkm::cont::CellLocatorTwoLevel locator2L;
locator2L.SetDensityL1(64.0f);
locator2L.SetDensityL2(1.0f);
TestCellLocator(locator2L, vtkm::Id3(8), 512); // 3D dataset
TestCellLocator(locator2L, vtkm::Id2(18), 512); // 2D dataset
//Test vtkm::cont::CellLocatorUniformBins
vtkm::cont::CellLocatorUniformBins locatorUB;
locatorUB.SetDims({ 32, 32, 32 });
TestCellLocator(locatorUB, vtkm::Id3(8), 512); // 3D dataset
TestCellLocator(locatorUB, vtkm::Id2(18), 512); // 2D dataset
//Test 2D dataset with 2D bins.
locatorUB.SetDims({ 32, 32, 1 });
TestCellLocator(locatorUB, vtkm::Id2(18), 512); // 2D dataset
}
} // anonymous
int UnitTestCellLocatorTwoLevel(int argc, char* argv[])
int UnitTestCellLocatorUnstructured(int argc, char* argv[])
{
return vtkm::cont::testing::Testing::Run(TestingCellLocatorTwoLevel, argc, argv);
return vtkm::cont::testing::Testing::Run(TestingCellLocatorUnstructured, argc, argv);
}

@ -21,6 +21,7 @@ set(headers
CellLocatorPartitioned.h
CellLocatorRectilinearGrid.h
CellLocatorTwoLevel.h
CellLocatorUniformBins.h
CellLocatorUniformGrid.h
CellMeasure.h
ColorTable.h

@ -0,0 +1,266 @@
//============================================================================
// 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.
//============================================================================
#ifndef vtk_m_exec_CellLocatorUniformBins_h
#define vtk_m_exec_CellLocatorUniformBins_h
#include <vtkm/exec/CellInside.h>
#include <vtkm/exec/ParametricCoordinates.h>
#include <vtkm/cont/ArrayHandle.h>
#include <vtkm/cont/CoordinateSystem.h>
#include <vtkm/TopologyElementTag.h>
#include <vtkm/VecFromPortalPermute.h>
#include <vtkm/VecTraits.h>
namespace vtkm
{
namespace exec
{
//--------------------------------------------------------------------
template <typename CellStructureType>
class VTKM_ALWAYS_EXPORT CellLocatorUniformBins
{
template <typename T>
using ReadPortal = typename vtkm::cont::ArrayHandle<T>::ReadPortalType;
using CoordsPortalType =
typename vtkm::cont::CoordinateSystem::MultiplexerArrayType::ReadPortalType;
using CellIdArrayType = vtkm::cont::ArrayHandle<vtkm::Id>;
using CellIdOffsetArrayType = vtkm::cont::ArrayHandle<vtkm::Id>;
using CellIdReadPortal =
typename vtkm::cont::ArrayHandleGroupVecVariable<CellIdArrayType,
CellIdOffsetArrayType>::ReadPortalType;
public:
template <typename CellSetType>
VTKM_CONT CellLocatorUniformBins(
const vtkm::Id3& cellDims,
const vtkm::Vec3f& origin,
const vtkm::Vec3f& maxPoint,
const vtkm::Vec3f& invSpacing,
const vtkm::Id3& maxCellIds,
const vtkm::cont::ArrayHandleGroupVecVariable<CellIdArrayType, CellIdOffsetArrayType>& cellIds,
const CellSetType& cellSet,
const vtkm::cont::CoordinateSystem& coords,
vtkm::cont::DeviceAdapterId device,
vtkm::cont::Token& token)
: CellDims(cellDims)
, Origin(origin)
, MaxPoint(maxPoint)
, InvSpacing(invSpacing)
, MaxCellIds(maxCellIds)
, CellIds(cellIds.PrepareForInput(device, token))
, CellSet(cellSet.PrepareForInput(device,
vtkm::TopologyElementTagCell{},
vtkm::TopologyElementTagPoint{},
token))
, Coords(coords.GetDataAsMultiplexer().PrepareForInput(device, token))
{
}
struct LastCell
{
vtkm::Id CellId = -1;
vtkm::Id BinIdx = -1;
};
VTKM_EXEC
vtkm::ErrorCode FindCell(const vtkm::Vec3f& point,
vtkm::Id& cellId,
vtkm::Vec3f& parametric) const
{
LastCell lastCell;
return this->FindCellImpl(point, cellId, parametric, lastCell);
}
VTKM_EXEC
vtkm::ErrorCode FindCell(const vtkm::Vec3f& point,
vtkm::Id& cellId,
vtkm::Vec3f& parametric,
LastCell& lastCell) const
{
//See if point is inside the last cell.
vtkm::Vec3f pc;
if ((lastCell.CellId >= 0) && (lastCell.CellId < this->CellSet.GetNumberOfElements()) &&
this->PointInCell(point, lastCell.CellId, pc) == vtkm::ErrorCode::Success)
{
parametric = pc;
cellId = lastCell.CellId;
return vtkm::ErrorCode::Success;
}
//See if it's in the last bin.
if ((lastCell.BinIdx >= 0) && (lastCell.BinIdx < this->CellIds.GetNumberOfValues()) &&
this->PointInBin(point, lastCell.BinIdx, cellId, pc) == vtkm::ErrorCode::Success)
{
parametric = pc;
lastCell.CellId = cellId;
return vtkm::ErrorCode::Success;
}
return this->FindCellImpl(point, cellId, parametric, lastCell);
}
VTKM_DEPRECATED(1.6, "Locators are no longer pointers. Use . operator.")
VTKM_EXEC CellLocatorUniformBins* operator->() { return this; }
VTKM_DEPRECATED(1.6, "Locators are no longer pointers. Use . operator.")
VTKM_EXEC const CellLocatorUniformBins* operator->() const { return this; }
private:
VTKM_EXEC bool IsInside(const vtkm::Vec3f& point) const
{
if (point[0] < this->Origin[0] || point[0] > this->MaxPoint[0])
return false;
if (point[1] < this->Origin[1] || point[1] > this->MaxPoint[1])
return false;
if (point[2] < this->Origin[2] || point[2] > this->MaxPoint[2])
return false;
return true;
}
VTKM_EXEC
vtkm::ErrorCode FindCellImpl(const vtkm::Vec3f& point,
vtkm::Id& cellId,
vtkm::Vec3f& parametric,
LastCell& lastCell) const
{
lastCell.CellId = -1;
lastCell.BinIdx = -1;
if (!this->IsInside(point))
{
cellId = -1;
return vtkm::ErrorCode::CellNotFound;
}
//Find the bin containing the point.
vtkm::Id3 logicalCell(0, 0, 0);
vtkm::Vec3f temp;
temp = point - this->Origin;
temp = temp * this->InvSpacing;
//make sure that if we border the upper edge, we sample the correct cell
logicalCell = vtkm::Min(vtkm::Id3(temp), this->MaxCellIds);
vtkm::Id binIdx =
(logicalCell[2] * this->CellDims[1] + logicalCell[1]) * this->CellDims[0] + logicalCell[0];
vtkm::Vec3f pc;
if (this->PointInBin(point, binIdx, cellId, pc) == vtkm::ErrorCode::Success)
{
parametric = pc;
lastCell.CellId = cellId;
lastCell.BinIdx = binIdx;
return vtkm::ErrorCode::Success;
}
return vtkm::ErrorCode::CellNotFound;
}
template <typename PointsVecType>
VTKM_EXEC vtkm::Bounds ComputeCellBounds(const PointsVecType& points) const
{
auto numPoints = vtkm::VecTraits<PointsVecType>::GetNumberOfComponents(points);
vtkm::Bounds bounds;
for (vtkm::IdComponent i = 0; i < numPoints; ++i)
bounds.Include(points[i]);
return bounds;
}
// TODO: This function may return false positives for non 3D cells as the
// tests are done on the projection of the point on the cell. Extra checks
// should be added to test if the point actually falls on the cell.
template <typename CellShapeTag, typename CoordsType>
VTKM_EXEC vtkm::ErrorCode PointInsideCell(vtkm::Vec3f point,
CellShapeTag cellShape,
CoordsType cellPoints,
vtkm::Vec3f& parametricCoordinates,
bool& inside) const
{
auto bounds = this->ComputeCellBounds(cellPoints);
if (bounds.Contains(point))
{
VTKM_RETURN_ON_ERROR(vtkm::exec::WorldCoordinatesToParametricCoordinates(
cellPoints, point, cellShape, parametricCoordinates));
inside = vtkm::exec::CellInside(parametricCoordinates, cellShape);
}
else
{
inside = false;
}
// Return success error code even point is not inside this cell
return vtkm::ErrorCode::Success;
}
VTKM_EXEC
vtkm::ErrorCode PointInBin(const vtkm::Vec3f& point,
const vtkm::Id& binIdx,
vtkm::Id& cellId,
vtkm::Vec3f& parametric) const
{
auto binIds = this->CellIds.Get(binIdx);
for (vtkm::IdComponent i = 0; i < binIds.GetNumberOfComponents(); i++)
{
vtkm::Id cid = binIds[i];
vtkm::Vec3f pc;
if (this->PointInCell(point, cid, pc) == vtkm::ErrorCode::Success)
{
cellId = cid;
parametric = pc;
return vtkm::ErrorCode::Success;
}
}
return vtkm::ErrorCode::CellNotFound;
}
VTKM_EXEC
vtkm::ErrorCode PointInCell(const vtkm::Vec3f& point,
const vtkm::Id& cid,
vtkm::Vec3f& parametric) const
{
auto indices = this->CellSet.GetIndices(cid);
auto pts = vtkm::make_VecFromPortalPermute(&indices, this->Coords);
vtkm::Vec3f pc;
bool inside;
auto status = this->PointInsideCell(point, this->CellSet.GetCellShape(cid), pts, pc, inside);
if (status == vtkm::ErrorCode::Success && inside)
{
parametric = pc;
return vtkm::ErrorCode::Success;
}
return vtkm::ErrorCode::CellNotFound;
}
vtkm::Id3 CellDims;
vtkm::Vec3f Origin;
vtkm::Vec3f MaxPoint;
vtkm::Vec3f InvSpacing;
vtkm::Id3 MaxCellIds;
CellIdReadPortal CellIds;
CellStructureType CellSet;
CoordsPortalType Coords;
};
}
} // vtkm::exec
#endif //vtk_m_exec_CellLocatorUniformBins_h