mirror of
https://gitlab.kitware.com/vtk/vtk-m
synced 2024-09-16 17:22:55 +00:00
Adding binary search for Rectilinear locator
This commit is contained in:
parent
06d05196b5
commit
5776945798
@ -59,26 +59,35 @@ public:
|
|||||||
// Linear search in the coordinates.
|
// Linear search in the coordinates.
|
||||||
vtkm::Id index;
|
vtkm::Id index;
|
||||||
/*Get floor X location*/
|
/*Get floor X location*/
|
||||||
for (index = 0; index < this->Dims[0] - 1; index++)
|
if (point[0] == xAxis.Get(this->Dims[0] - 1))
|
||||||
if (xAxis.Get(index) <= point[0] && point[0] < xAxis.Get(index + 1))
|
logical[0] = this->Dims[0] - 1;
|
||||||
{
|
else
|
||||||
logical[0] = index;
|
for (index = 0; index < this->Dims[0] - 1; index++)
|
||||||
break;
|
if (xAxis.Get(index) <= point[0] && point[0] < xAxis.Get(index + 1))
|
||||||
}
|
{
|
||||||
|
logical[0] = index;
|
||||||
|
break;
|
||||||
|
}
|
||||||
/*Get floor Y location*/
|
/*Get floor Y location*/
|
||||||
for (index = 0; index < this->Dims[1] - 1; index++)
|
if (point[1] == yAxis.Get(this->Dims[1] - 1))
|
||||||
if (yAxis.Get(index) <= point[1] && point[1] < yAxis.Get(index + 1))
|
logical[1] = this->Dims[1] - 1;
|
||||||
{
|
else
|
||||||
logical[1] = index;
|
for (index = 0; index < this->Dims[1] - 1; index++)
|
||||||
break;
|
if (yAxis.Get(index) <= point[1] && point[1] < yAxis.Get(index + 1))
|
||||||
}
|
{
|
||||||
|
logical[1] = index;
|
||||||
|
break;
|
||||||
|
}
|
||||||
/*Get floor Z location*/
|
/*Get floor Z location*/
|
||||||
for (index = 0; index < this->Dims[2] - 1; index++)
|
if (point[2] == zAxis.Get(this->Dims[2] - 1))
|
||||||
if (zAxis.Get(index) <= point[2] && point[2] < zAxis.Get(index + 1))
|
logical[2] = this->Dims[2] - 1;
|
||||||
{
|
else
|
||||||
logical[2] = index;
|
for (index = 0; index < this->Dims[2] - 1; index++)
|
||||||
break;
|
if (zAxis.Get(index) <= point[2] && point[2] < zAxis.Get(index + 1))
|
||||||
}
|
{
|
||||||
|
logical[2] = index;
|
||||||
|
break;
|
||||||
|
}
|
||||||
if (logical[0] == -1 || logical[1] == -1 || logical[2] == -1)
|
if (logical[0] == -1 || logical[1] == -1 || logical[2] == -1)
|
||||||
return -1;
|
return -1;
|
||||||
return logical[2] * (Dims[0] - 1) * (Dims[1] - 1) + logical[1] * (Dims[0] - 1) + logical[0];
|
return logical[2] * (Dims[0] - 1) * (Dims[1] - 1) + logical[1] * (Dims[0] - 1) + logical[0];
|
||||||
@ -136,12 +145,8 @@ public:
|
|||||||
vtkm::cont::CoordinateSystem coords = dataset.GetCoordinateSystem();
|
vtkm::cont::CoordinateSystem coords = dataset.GetCoordinateSystem();
|
||||||
vtkm::cont::DynamicCellSet cellSet = dataset.GetCellSet();
|
vtkm::cont::DynamicCellSet cellSet = dataset.GetCellSet();
|
||||||
vtkm::Bounds bounds = coords.GetBounds();
|
vtkm::Bounds bounds = coords.GetBounds();
|
||||||
std::cout << "X bounds : " << bounds.X.Min << " to " << bounds.X.Max << std::endl;
|
|
||||||
std::cout << "Y bounds : " << bounds.Y.Min << " to " << bounds.Y.Max << std::endl;
|
|
||||||
std::cout << "Z bounds : " << bounds.Z.Min << " to " << bounds.Z.Max << std::endl;
|
|
||||||
vtkm::Id3 dims =
|
vtkm::Id3 dims =
|
||||||
cellSet.Cast<StructuredType>().GetSchedulingRange(vtkm::TopologyElementTagPoint());
|
cellSet.Cast<StructuredType>().GetSchedulingRange(vtkm::TopologyElementTagPoint());
|
||||||
std::cout << "Dimensions of dataset : " << dims << std::endl;
|
|
||||||
|
|
||||||
// Generate some sample points.
|
// Generate some sample points.
|
||||||
using PointType = vtkm::Vec3f;
|
using PointType = vtkm::Vec3f;
|
||||||
@ -180,7 +185,6 @@ public:
|
|||||||
{
|
{
|
||||||
VTKM_TEST_ASSERT(matchPortal.Get(index), "Points do not match");
|
VTKM_TEST_ASSERT(matchPortal.Get(index), "Points do not match");
|
||||||
}
|
}
|
||||||
std::cout << "Test finished successfully." << std::endl;
|
|
||||||
}
|
}
|
||||||
|
|
||||||
void operator()() const
|
void operator()() const
|
||||||
|
@ -57,22 +57,23 @@ public:
|
|||||||
, Coords(coords.PrepareForInput(DeviceAdapter()))
|
, Coords(coords.PrepareForInput(DeviceAdapter()))
|
||||||
, PointDimensions(cellSet.GetPointDimensions())
|
, PointDimensions(cellSet.GetPointDimensions())
|
||||||
{
|
{
|
||||||
auto contPortal = coords.GetPortalConstControl();
|
this->AxisPortals[0] = Coords.GetFirstPortal();
|
||||||
for (vtkm::IdComponent dim = 0; dim < dimensions; dim++)
|
this->MinPoint[0] = coords.GetPortalConstControl().GetFirstPortal().Get(0);
|
||||||
|
this->MaxPoint[0] = coords.GetPortalConstControl().GetFirstPortal().Get(PointDimensions[0] - 1);
|
||||||
|
|
||||||
|
this->AxisPortals[1] = Coords.GetSecondPortal();
|
||||||
|
this->MinPoint[1] = coords.GetPortalConstControl().GetSecondPortal().Get(0);
|
||||||
|
this->MaxPoint[1] =
|
||||||
|
coords.GetPortalConstControl().GetSecondPortal().Get(PointDimensions[1] - 1);
|
||||||
|
if (dimensions == 3)
|
||||||
{
|
{
|
||||||
auto componentContPortal = (dim == 0) ? contPortal.GetFirstPortal() : (dim == 1)
|
this->AxisPortals[2] = Coords.GetThirdPortal();
|
||||||
? contPortal.GetSecondPortal()
|
this->MinPoint[2] = coords.GetPortalConstControl().GetThirdPortal().Get(0);
|
||||||
: contPortal.GetThirdPortal();
|
this->MaxPoint[2] =
|
||||||
auto componentExecPortal = (dim == 0) ? Coords.GetFirstPortal() : (dim == 1)
|
coords.GetPortalConstControl().GetThirdPortal().Get(PointDimensions[2] - 1);
|
||||||
? Coords.GetSecondPortal()
|
|
||||||
: Coords.GetThirdPortal();
|
|
||||||
this->AxisPortals[dim] = componentExecPortal;
|
|
||||||
this->MinPoint[dim] = componentContPortal.Get(0);
|
|
||||||
this->MaxPoint[dim] = componentContPortal.Get(PointDimensions[0] - 1);
|
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
VTKM_EXEC_CONT virtual ~CellLocatorRectilinearGrid() noexcept
|
VTKM_EXEC_CONT virtual ~CellLocatorRectilinearGrid() noexcept
|
||||||
{
|
{
|
||||||
// This must not be defaulted, since defaulted virtual destructors are
|
// This must not be defaulted, since defaulted virtual destructors are
|
||||||
@ -100,6 +101,7 @@ public:
|
|||||||
vtkm::Vec3f& parametric,
|
vtkm::Vec3f& parametric,
|
||||||
const vtkm::exec::FunctorBase& worklet) const override
|
const vtkm::exec::FunctorBase& worklet) const override
|
||||||
{
|
{
|
||||||
|
(void)worklet; //suppress unused warning
|
||||||
if (!IsInside(point))
|
if (!IsInside(point))
|
||||||
{
|
{
|
||||||
cellId = -1;
|
cellId = -1;
|
||||||
@ -121,47 +123,35 @@ public:
|
|||||||
continue;
|
continue;
|
||||||
}
|
}
|
||||||
|
|
||||||
bool found = false;
|
vtkm::Id minIndex = 0;
|
||||||
vtkm::FloatDefault minVal = this->AxisPortals[dim].Get(logicalCell[dim]);
|
vtkm::Id maxIndex = PointDimensions[dim] - 1;
|
||||||
const vtkm::Id searchDir = (point[dim] - minVal >= 0.f) ? 1 : -1;
|
vtkm::FloatDefault minVal;
|
||||||
vtkm::FloatDefault maxVal = this->AxisPortals[dim].Get(logicalCell[dim] + 1);
|
vtkm::FloatDefault maxVal;
|
||||||
|
minVal = this->AxisPortals[dim].Get(minIndex);
|
||||||
while (!found)
|
maxVal = this->AxisPortals[dim].Get(maxIndex);
|
||||||
|
while (maxIndex > minIndex + 1)
|
||||||
{
|
{
|
||||||
if (point[dim] >= minVal && point[dim] < maxVal)
|
vtkm::Id midIndex = (minIndex + maxIndex) / 2;
|
||||||
|
vtkm::FloatDefault midVal = this->AxisPortals[dim].Get(midIndex);
|
||||||
|
if (point[dim] <= midVal)
|
||||||
{
|
{
|
||||||
found = true;
|
maxIndex = midIndex;
|
||||||
continue;
|
maxVal = midVal;
|
||||||
}
|
|
||||||
|
|
||||||
logicalCell[dim] += searchDir;
|
|
||||||
vtkm::Id nextCellId = searchDir == 1 ? logicalCell[dim] + 1 : logicalCell[dim];
|
|
||||||
vtkm::FloatDefault next = this->AxisPortals[dim].Get(nextCellId);
|
|
||||||
if (searchDir == 1)
|
|
||||||
{
|
|
||||||
minVal = maxVal;
|
|
||||||
maxVal = next;
|
|
||||||
}
|
}
|
||||||
else
|
else
|
||||||
{
|
{
|
||||||
maxVal = minVal;
|
minIndex = midIndex;
|
||||||
minVal = next;
|
minVal = midVal;
|
||||||
}
|
}
|
||||||
}
|
}
|
||||||
|
logicalCell[dim] = minIndex;
|
||||||
|
//printf("Min Index [%d] : %lld\n", dim, minIndex);
|
||||||
|
//printf("Max Index [%d] : %lld\n", dim, maxIndex);
|
||||||
|
//printf("Logical [%d] : %lld\n", dim, logicalCell[dim]);
|
||||||
|
parametric[dim] = (point[dim] - minVal) / (maxVal - minVal);
|
||||||
}
|
}
|
||||||
|
|
||||||
// Get the actual cellId, from the logical cell index of the cell
|
// Get the actual cellId, from the logical cell index of the cell
|
||||||
cellId = logicalCell[2] * this->PlaneSize + logicalCell[1] * this->RowSize + logicalCell[0];
|
cellId = logicalCell[2] * this->PlaneSize + logicalCell[1] * this->RowSize + logicalCell[0];
|
||||||
|
|
||||||
bool success = false;
|
|
||||||
using IndicesType = typename CellSetPortal::IndicesType;
|
|
||||||
IndicesType cellPointIndices = this->CellSet.GetIndices(cellId);
|
|
||||||
vtkm::VecFromPortalPermute<IndicesType, RectilinearPortalType> cellPoints(&cellPointIndices,
|
|
||||||
Coords);
|
|
||||||
auto cellShape = this->CellSet.GetCellShape(cellId);
|
|
||||||
// Get Parametric Coordinates from the cell, for the point.
|
|
||||||
parametric = vtkm::exec::WorldCoordinatesToParametricCoordinates(
|
|
||||||
cellPoints, point, cellShape, success, worklet);
|
|
||||||
}
|
}
|
||||||
|
|
||||||
private:
|
private:
|
||||||
|
Loading…
Reference in New Issue
Block a user