VTK-m  2.0
exec/PointLocatorSparseGrid.h
Go to the documentation of this file.
1 //============================================================================
2 // Copyright (c) Kitware, Inc.
3 // All rights reserved.
4 // See LICENSE.txt for details.
5 //
6 // This software is distributed WITHOUT ANY WARRANTY; without even
7 // the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
8 // PURPOSE. See the above copyright notice for more information.
9 //============================================================================
10 #ifndef vtk_m_exec_PointLocatorSparseGrid_h
11 #define vtk_m_exec_PointLocatorSparseGrid_h
12 
14 
15 #include <vtkm/VectorAnalysis.h>
16 
17 namespace vtkm
18 {
19 namespace exec
20 {
21 
23 {
24 public:
25  using CoordPortalType =
26  typename vtkm::cont::CoordinateSystem::MultiplexerArrayType::ReadPortalType;
28 
29 
31  const vtkm::Vec3f& max,
32  const vtkm::Id3& dims,
33  const CoordPortalType& coords,
34  const IdPortalType& pointIds,
35  const IdPortalType& cellLower,
36  const IdPortalType& cellUpper)
37  : Min(min)
38  , Dims(dims)
39  , Dxdydz((max - Min) / Dims)
40  , Coords(coords)
41  , PointIds(pointIds)
42  , CellLower(cellLower)
43  , CellUpper(cellUpper)
44  {
45  }
46 
57  VTKM_EXEC void FindNearestNeighbor(const vtkm::Vec3f& queryPoint,
58  vtkm::Id& nearestNeighborId,
59  vtkm::FloatDefault& distance2) const
60  {
61  //std::cout << "FindNeareastNeighbor: " << queryPoint << std::endl;
62  vtkm::Id3 ijk = (queryPoint - this->Min) / this->Dxdydz;
63  ijk = vtkm::Max(ijk, vtkm::Id3(0));
64  ijk = vtkm::Min(ijk, this->Dims - vtkm::Id3(1));
65 
66  nearestNeighborId = -1;
67  distance2 = vtkm::Infinity<vtkm::FloatDefault>();
68 
69  this->FindInCell(queryPoint, ijk, nearestNeighborId, distance2);
70 
71  // TODO: This might stop looking before the absolute nearest neighbor is found.
72  vtkm::Id maxLevel = vtkm::Max(vtkm::Max(this->Dims[0], this->Dims[1]), this->Dims[2]);
73  vtkm::Id level;
74  for (level = 1; (nearestNeighborId < 0) && (level < maxLevel); ++level)
75  {
76  this->FindInBox(queryPoint, ijk, level, nearestNeighborId, distance2);
77  }
78 
79  // Search one more level out. This is still not guaranteed to find the closest point
80  // in all cases (past level 2), but it will catch most cases where the closest point
81  // is just on the other side of a cell boundary.
82  this->FindInBox(queryPoint, ijk, level, nearestNeighborId, distance2);
83  }
84 
85 private:
89 
91 
95 
96  VTKM_EXEC void FindInCell(const vtkm::Vec3f& queryPoint,
97  const vtkm::Id3& ijk,
98  vtkm::Id& nearestNeighborId,
99  vtkm::FloatDefault& nearestDistance2) const
100  {
101  vtkm::Id cellId = ijk[0] + (ijk[1] * this->Dims[0]) + (ijk[2] * this->Dims[0] * this->Dims[1]);
102  vtkm::Id lower = this->CellLower.Get(cellId);
103  vtkm::Id upper = this->CellUpper.Get(cellId);
104  for (vtkm::Id index = lower; index < upper; index++)
105  {
106  vtkm::Id pointid = this->PointIds.Get(index);
107  vtkm::Vec3f point = this->Coords.Get(pointid);
108  vtkm::FloatDefault distance2 = vtkm::MagnitudeSquared(point - queryPoint);
109  if (distance2 < nearestDistance2)
110  {
111  nearestNeighborId = pointid;
112  nearestDistance2 = distance2;
113  }
114  }
115  }
116 
117  VTKM_EXEC void FindInBox(const vtkm::Vec3f& queryPoint,
118  const vtkm::Id3& boxCenter,
119  vtkm::Id level,
120  vtkm::Id& nearestNeighborId,
121  vtkm::FloatDefault& nearestDistance2) const
122  {
123  if ((boxCenter[0] - level) >= 0)
124  {
125  this->FindInXPlane(
126  queryPoint, boxCenter - vtkm::Id3(level, 0, 0), level, nearestNeighborId, nearestDistance2);
127  }
128  if ((boxCenter[0] + level) < this->Dims[0])
129  {
130  this->FindInXPlane(
131  queryPoint, boxCenter + vtkm::Id3(level, 0, 0), level, nearestNeighborId, nearestDistance2);
132  }
133 
134  if ((boxCenter[1] - level) >= 0)
135  {
136  this->FindInYPlane(
137  queryPoint, boxCenter - vtkm::Id3(0, level, 0), level, nearestNeighborId, nearestDistance2);
138  }
139  if ((boxCenter[1] + level) < this->Dims[1])
140  {
141  this->FindInYPlane(
142  queryPoint, boxCenter + vtkm::Id3(0, level, 0), level, nearestNeighborId, nearestDistance2);
143  }
144 
145  if ((boxCenter[2] - level) >= 0)
146  {
147  this->FindInZPlane(
148  queryPoint, boxCenter - vtkm::Id3(0, 0, level), level, nearestNeighborId, nearestDistance2);
149  }
150  if ((boxCenter[2] + level) < this->Dims[2])
151  {
152  this->FindInZPlane(
153  queryPoint, boxCenter + vtkm::Id3(0, 0, level), level, nearestNeighborId, nearestDistance2);
154  }
155  }
156 
157  VTKM_EXEC void FindInPlane(const vtkm::Vec3f& queryPoint,
158  const vtkm::Id3& planeCenter,
159  const vtkm::Id3& div,
160  const vtkm::Id3& mod,
161  const vtkm::Id3& origin,
162  vtkm::Id numInPlane,
163  vtkm::Id& nearestNeighborId,
164  vtkm::FloatDefault& nearestDistance2) const
165  {
166  for (vtkm::Id index = 0; index < numInPlane; ++index)
167  {
168  vtkm::Id3 ijk = planeCenter + vtkm::Id3(index) / div +
169  vtkm::Id3(index % mod[0], index % mod[1], index % mod[2]) + origin;
170  if ((ijk[0] >= 0) && (ijk[0] < this->Dims[0]) && (ijk[1] >= 0) && (ijk[1] < this->Dims[1]) &&
171  (ijk[2] >= 0) && (ijk[2] < this->Dims[2]))
172  {
173  this->FindInCell(queryPoint, ijk, nearestNeighborId, nearestDistance2);
174  }
175  }
176  }
177 
178  VTKM_EXEC void FindInXPlane(const vtkm::Vec3f& queryPoint,
179  const vtkm::Id3& planeCenter,
180  vtkm::Id level,
181  vtkm::Id& nearestNeighborId,
182  vtkm::FloatDefault& nearestDistance2) const
183  {
184  vtkm::Id yWidth = (2 * level) + 1;
185  vtkm::Id zWidth = (2 * level) + 1;
186  vtkm::Id3 div = { yWidth * zWidth, yWidth * zWidth, yWidth };
187  vtkm::Id3 mod = { 1, yWidth, 1 };
188  vtkm::Id3 origin = { 0, -level, -level };
189  vtkm::Id numInPlane = yWidth * zWidth;
190  this->FindInPlane(
191  queryPoint, planeCenter, div, mod, origin, numInPlane, nearestNeighborId, nearestDistance2);
192  }
193 
194  VTKM_EXEC void FindInYPlane(const vtkm::Vec3f& queryPoint,
195  vtkm::Id3 planeCenter,
196  vtkm::Id level,
197  vtkm::Id& nearestNeighborId,
198  vtkm::FloatDefault& nearestDistance2) const
199  {
200  vtkm::Id xWidth = (2 * level) - 1;
201  vtkm::Id zWidth = (2 * level) + 1;
202  vtkm::Id3 div = { xWidth * zWidth, xWidth * zWidth, xWidth };
203  vtkm::Id3 mod = { xWidth, 1, 1 };
204  vtkm::Id3 origin = { -level + 1, 0, -level };
205  vtkm::Id numInPlane = xWidth * zWidth;
206  this->FindInPlane(
207  queryPoint, planeCenter, div, mod, origin, numInPlane, nearestNeighborId, nearestDistance2);
208  }
209 
210  VTKM_EXEC void FindInZPlane(const vtkm::Vec3f& queryPoint,
211  vtkm::Id3 planeCenter,
212  vtkm::Id level,
213  vtkm::Id& nearestNeighborId,
214  vtkm::FloatDefault& nearestDistance2) const
215  {
216  vtkm::Id xWidth = (2 * level) - 1;
217  vtkm::Id yWidth = (2 * level) - 1;
218  vtkm::Id3 div = { xWidth * yWidth, xWidth, xWidth * yWidth };
219  vtkm::Id3 mod = { xWidth, 1, 1 };
220  vtkm::Id3 origin = { -level + 1, -level + 1, 0 };
221  vtkm::Id numInPlane = xWidth * yWidth;
222  this->FindInPlane(
223  queryPoint, planeCenter, div, mod, origin, numInPlane, nearestNeighborId, nearestDistance2);
224  }
225 };
226 
227 } // vtkm::exec
228 } // vtkm
229 
230 #endif // vtk_m_exec_PointLocatorSparseGrid_h
vtkm::exec::PointLocatorSparseGrid
Definition: exec/PointLocatorSparseGrid.h:22
vtkm::exec::PointLocatorSparseGrid::FindInZPlane
VTKM_EXEC void FindInZPlane(const vtkm::Vec3f &queryPoint, vtkm::Id3 planeCenter, vtkm::Id level, vtkm::Id &nearestNeighborId, vtkm::FloatDefault &nearestDistance2) const
Definition: exec/PointLocatorSparseGrid.h:210
vtkm::MagnitudeSquared
VTKM_EXEC_CONT detail::FloatingPointReturnType< T >::Type MagnitudeSquared(const T &x)
Returns the square of the magnitude of a vector.
Definition: VectorAnalysis.h:64
vtkm::exec::PointLocatorSparseGrid::FindNearestNeighbor
VTKM_EXEC void FindNearestNeighbor(const vtkm::Vec3f &queryPoint, vtkm::Id &nearestNeighborId, vtkm::FloatDefault &distance2) const
Nearest neighbor search using a Uniform Grid.
Definition: exec/PointLocatorSparseGrid.h:57
VTKM_EXEC
#define VTKM_EXEC
Definition: ExportMacros.h:51
vtkm
Groups connected points that have the same field value.
Definition: Atomic.h:19
vtkm::exec::PointLocatorSparseGrid::IdPortalType
typename vtkm::cont::ArrayHandle< vtkm::Id >::ReadPortalType IdPortalType
Definition: exec/PointLocatorSparseGrid.h:27
vtkm::exec::PointLocatorSparseGrid::FindInYPlane
VTKM_EXEC void FindInYPlane(const vtkm::Vec3f &queryPoint, vtkm::Id3 planeCenter, vtkm::Id level, vtkm::Id &nearestNeighborId, vtkm::FloatDefault &nearestDistance2) const
Definition: exec/PointLocatorSparseGrid.h:194
vtkm::exec::PointLocatorSparseGrid::CellLower
IdPortalType CellLower
Definition: exec/PointLocatorSparseGrid.h:93
vtkm::exec::PointLocatorSparseGrid::CellUpper
IdPortalType CellUpper
Definition: exec/PointLocatorSparseGrid.h:94
vtkm::exec::PointLocatorSparseGrid::FindInXPlane
VTKM_EXEC void FindInXPlane(const vtkm::Vec3f &queryPoint, const vtkm::Id3 &planeCenter, vtkm::Id level, vtkm::Id &nearestNeighborId, vtkm::FloatDefault &nearestDistance2) const
Definition: exec/PointLocatorSparseGrid.h:178
CoordinateSystem.h
vtkm::cont::ArrayHandle::ReadPortalType
typename StorageType::ReadPortalType ReadPortalType
Definition: ArrayHandle.h:294
vtkm::Id
vtkm::Int32 Id
Represents an ID (index into arrays).
Definition: Types.h:191
vtkm::exec::PointLocatorSparseGrid::FindInBox
VTKM_EXEC void FindInBox(const vtkm::Vec3f &queryPoint, const vtkm::Id3 &boxCenter, vtkm::Id level, vtkm::Id &nearestNeighborId, vtkm::FloatDefault &nearestDistance2) const
Definition: exec/PointLocatorSparseGrid.h:117
VectorAnalysis.h
vtkm::exec::PointLocatorSparseGrid::Min
vtkm::Vec3f Min
Definition: exec/PointLocatorSparseGrid.h:86
vtkm::exec::PointLocatorSparseGrid::Dims
vtkm::Id3 Dims
Definition: exec/PointLocatorSparseGrid.h:87
vtkm::Id3
vtkm::Vec< vtkm::Id, 3 > Id3
Id3 corresponds to a 3-dimensional index for 3d arrays.
Definition: Types.h:1003
vtkm::exec::PointLocatorSparseGrid::Dxdydz
vtkm::Vec3f Dxdydz
Definition: exec/PointLocatorSparseGrid.h:88
vtkm::Vec< vtkm::FloatDefault, 3 >
vtkm::FloatDefault
vtkm::Float32 FloatDefault
The floating point type to use when no other precision is specified.
Definition: Types.h:198
vtkm::exec::PointLocatorSparseGrid::PointIds
IdPortalType PointIds
Definition: exec/PointLocatorSparseGrid.h:92
VTKM_ALWAYS_EXPORT
#define VTKM_ALWAYS_EXPORT
Definition: ExportMacros.h:92
vtkm::exec::PointLocatorSparseGrid::FindInPlane
VTKM_EXEC void FindInPlane(const vtkm::Vec3f &queryPoint, const vtkm::Id3 &planeCenter, const vtkm::Id3 &div, const vtkm::Id3 &mod, const vtkm::Id3 &origin, vtkm::Id numInPlane, vtkm::Id &nearestNeighborId, vtkm::FloatDefault &nearestDistance2) const
Definition: exec/PointLocatorSparseGrid.h:157
vtkm::exec::PointLocatorSparseGrid::FindInCell
VTKM_EXEC void FindInCell(const vtkm::Vec3f &queryPoint, const vtkm::Id3 &ijk, vtkm::Id &nearestNeighborId, vtkm::FloatDefault &nearestDistance2) const
Definition: exec/PointLocatorSparseGrid.h:96
vtkm::exec::PointLocatorSparseGrid::Coords
CoordPortalType Coords
Definition: exec/PointLocatorSparseGrid.h:90
vtkm::exec::PointLocatorSparseGrid::PointLocatorSparseGrid
PointLocatorSparseGrid(const vtkm::Vec3f &min, const vtkm::Vec3f &max, const vtkm::Id3 &dims, const CoordPortalType &coords, const IdPortalType &pointIds, const IdPortalType &cellLower, const IdPortalType &cellUpper)
Definition: exec/PointLocatorSparseGrid.h:30
vtkm::exec::PointLocatorSparseGrid::CoordPortalType
typename vtkm::cont::CoordinateSystem::MultiplexerArrayType::ReadPortalType CoordPortalType
Definition: exec/PointLocatorSparseGrid.h:26