1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkStaticPointLocator.cxx
5
6 Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
7 All rights reserved.
8 See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
9
10 This software is distributed WITHOUT ANY WARRANTY; without even
11 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
12 PURPOSE. See the above copyright notice for more information.
13
14 =========================================================================*/
15 #include "vtkStaticPointLocator.h"
16
17 #include "vtkBoundingBox.h"
18 #include "vtkBox.h"
19 #include "vtkCellArray.h"
20 #include "vtkIdList.h"
21 #include "vtkIntArray.h"
22 #include "vtkLine.h"
23 #include "vtkMath.h"
24 #include "vtkObjectFactory.h"
25 #include "vtkPoints.h"
26 #include "vtkPolyData.h"
27 #include "vtkSMPThreadLocalObject.h"
28 #include "vtkSMPTools.h"
29
30 #include <vector>
31
32 vtkStandardNewMacro(vtkStaticPointLocator);
33
34 // There are stack-allocated bucket neighbor lists. This is the initial
35 // value. Too small and heap allocation kicks in.
36 #define VTK_INITIAL_BUCKET_SIZE 10000
37
38 //------------------------------------------------------------------------------
39 // The following code supports threaded point locator construction. The locator
40 // is assumed to be constructed once (i.e., it does not allow incremental point
41 // insertion). The algorithm proceeds in three steps:
42 // 1) All points are assigned a bucket index (combined i-j-k bucket location).
43 // The index is computed in parallel. This requires a one time allocation of an
44 // index array (which is also associated with the originating point ids).
45 // 2) vtkSMPTools::Sort() is used to sort the index array. Note that the sort
46 // carries along the point ids as well. This creates contiguous runs of points
47 // all resident in the same bucket.
48 // 3) The bucket offsets are updated to refer to the right entry location into
49 // the sorted point ids array. This enables quick access, and an indirect count
50 // of the number of points in each bucket.
51
52 // Believe it or not I had to change the name because MS Visual Studio was
53 // mistakenly linking the hidden, scoped classes (vtkNeighborBuckets) found
54 // in vtkPointLocator and vtkStaticPointLocator and causing weird faults.
55 struct NeighborBuckets;
56
57 //------------------------------------------------------------------------------
58 // The bucketed points, including the sorted map. This is just a PIMPLd
59 // wrapper around the classes that do the real work.
60 struct vtkBucketList
61 {
62 vtkStaticPointLocator* Locator; // locater
63 vtkIdType NumPts; // the number of points to bucket
64 vtkIdType NumBuckets;
65 int BatchSize;
66
67 // These are internal data members used for performance reasons
68 vtkDataSet* DataSet;
69 int Divisions[3];
70 double Bounds[6];
71 double H[3];
72 double hX, hY, hZ;
73 double fX, fY, fZ, bX, bY, bZ;
74 vtkIdType xD, yD, zD, xyD;
75
76 // Construction
vtkBucketListvtkBucketList77 vtkBucketList(vtkStaticPointLocator* loc, vtkIdType numPts, int numBuckets)
78 {
79 this->Locator = loc;
80 this->NumPts = numPts;
81 this->NumBuckets = numBuckets;
82 this->BatchSize = 10000; // building the offset array
83 this->DataSet = loc->GetDataSet();
84 loc->GetDivisions(this->Divisions);
85
86 // Setup internal data members for more efficient processing.
87 double spacing[3], bounds[6];
88 loc->GetDivisions(this->Divisions);
89 loc->GetSpacing(spacing);
90 loc->GetBounds(bounds);
91 this->hX = this->H[0] = spacing[0];
92 this->hY = this->H[1] = spacing[1];
93 this->hZ = this->H[2] = spacing[2];
94 this->fX = 1.0 / spacing[0];
95 this->fY = 1.0 / spacing[1];
96 this->fZ = 1.0 / spacing[2];
97 this->bX = this->Bounds[0] = bounds[0];
98 this->Bounds[1] = bounds[1];
99 this->bY = this->Bounds[2] = bounds[2];
100 this->Bounds[3] = bounds[3];
101 this->bZ = this->Bounds[4] = bounds[4];
102 this->Bounds[5] = bounds[5];
103 this->xD = this->Divisions[0];
104 this->yD = this->Divisions[1];
105 this->zD = this->Divisions[2];
106 this->xyD = this->Divisions[0] * this->Divisions[1];
107 }
108
109 // Virtuals for templated subclasses
110 virtual ~vtkBucketList() = default;
111 virtual void BuildLocator() = 0;
112
113 // place points in appropriate buckets
114 void GetBucketNeighbors(
115 NeighborBuckets* buckets, const int ijk[3], const int ndivs[3], int level);
116 void GenerateFace(int face, int i, int j, int k, vtkPoints* pts, vtkCellArray* polys);
117 double Distance2ToBucket(const double x[3], const int nei[3]);
118 double Distance2ToBounds(const double x[3], const double bounds[6]);
119
120 //-----------------------------------------------------------------------------
121 // Inlined for performance. These function invocations must be called after
122 // BuildLocator() is invoked, otherwise the output is indeterminate.
GetBucketIndicesvtkBucketList123 void GetBucketIndices(const double* x, int ijk[3]) const
124 {
125 // Compute point index. Make sure it lies within range of locator.
126 vtkIdType tmp0 = static_cast<vtkIdType>(((x[0] - bX) * fX));
127 vtkIdType tmp1 = static_cast<vtkIdType>(((x[1] - bY) * fY));
128 vtkIdType tmp2 = static_cast<vtkIdType>(((x[2] - bZ) * fZ));
129
130 ijk[0] = tmp0 < 0 ? 0 : (tmp0 >= xD ? xD - 1 : tmp0);
131 ijk[1] = tmp1 < 0 ? 0 : (tmp1 >= yD ? yD - 1 : tmp1);
132 ijk[2] = tmp2 < 0 ? 0 : (tmp2 >= zD ? zD - 1 : tmp2);
133 }
134
135 //-----------------------------------------------------------------------------
GetBucketIndexvtkBucketList136 vtkIdType GetBucketIndex(const double* x) const
137 {
138 int ijk[3];
139 this->GetBucketIndices(x, ijk);
140 return ijk[0] + ijk[1] * xD + ijk[2] * xyD;
141 }
142 };
143
144 //------------------------------------------------------------------------------
145 // Utility class to store an array of ijk values
146 struct NeighborBuckets
147 {
NeighborBucketsNeighborBuckets148 NeighborBuckets()
149 {
150 this->Count = 0;
151 this->P = this->InitialBuffer;
152 this->MaxSize = VTK_INITIAL_BUCKET_SIZE;
153 }
~NeighborBucketsNeighborBuckets154 ~NeighborBuckets()
155 {
156 this->Count = 0;
157 if (this->P != this->InitialBuffer)
158 {
159 delete[] this->P;
160 }
161 }
GetNumberOfNeighborsNeighborBuckets162 int GetNumberOfNeighbors() { return this->Count; }
ResetNeighborBuckets163 void Reset() { this->Count = 0; }
164
GetPointNeighborBuckets165 int* GetPoint(vtkIdType i) { return this->P + 3 * i; }
166
InsertNextBucketNeighborBuckets167 vtkIdType InsertNextBucket(const int x[3])
168 {
169 // Re-allocate if beyond the current max size.
170 // (Increase by VTK_INITIAL_BUCKET_SIZE)
171 int* tmp;
172 vtkIdType offset = this->Count * 3;
173
174 if (this->Count >= this->MaxSize)
175 {
176 tmp = this->P;
177 this->MaxSize *= 2;
178 this->P = new int[this->MaxSize * 3];
179
180 memcpy(this->P, tmp, offset * sizeof(int));
181
182 if (tmp != this->InitialBuffer)
183 {
184 delete[] tmp;
185 }
186 }
187
188 tmp = this->P + offset;
189 *tmp++ = *x++;
190 *tmp++ = *x++;
191 *tmp = *x;
192 this->Count++;
193 return this->Count - 1;
194 }
195
196 protected:
197 // Start with an array to avoid memory allocation overhead
198 int InitialBuffer[VTK_INITIAL_BUCKET_SIZE * 3];
199 int* P;
200 vtkIdType Count;
201 vtkIdType MaxSize;
202 };
203
204 //------------------------------------------------------------------------------
205 // Internal function to get bucket neighbors at specified level
206 //
GetBucketNeighbors(NeighborBuckets * buckets,const int ijk[3],const int ndivs[3],int level)207 void vtkBucketList::GetBucketNeighbors(
208 NeighborBuckets* buckets, const int ijk[3], const int ndivs[3], int level)
209 {
210 int i, j, k, min, max, minLevel[3], maxLevel[3];
211 int nei[3];
212
213 // Initialize
214 //
215 buckets->Reset();
216
217 // If at this bucket, just place into list
218 //
219 if (level == 0)
220 {
221 buckets->InsertNextBucket(ijk);
222 return;
223 }
224
225 // Create permutations of the ijk indices that are at the level
226 // required. If these are legal buckets, add to list for searching.
227 //
228 for (i = 0; i < 3; i++)
229 {
230 min = ijk[i] - level;
231 max = ijk[i] + level;
232 minLevel[i] = (min > 0 ? min : 0);
233 maxLevel[i] = (max < (ndivs[i] - 1) ? max : (ndivs[i] - 1));
234 }
235
236 for (i = minLevel[0]; i <= maxLevel[0]; i++)
237 {
238 for (j = minLevel[1]; j <= maxLevel[1]; j++)
239 {
240 for (k = minLevel[2]; k <= maxLevel[2]; k++)
241 {
242 if (i == (ijk[0] + level) || i == (ijk[0] - level) || j == (ijk[1] + level) ||
243 j == (ijk[1] - level) || k == (ijk[2] + level) || k == (ijk[2] - level))
244 {
245 nei[0] = i;
246 nei[1] = j;
247 nei[2] = k;
248 buckets->InsertNextBucket(nei);
249 }
250 }
251 }
252 }
253 }
254
255 //------------------------------------------------------------------------------
GenerateFace(int face,int i,int j,int k,vtkPoints * pts,vtkCellArray * polys)256 void vtkBucketList::GenerateFace(int face, int i, int j, int k, vtkPoints* pts, vtkCellArray* polys)
257 {
258 vtkIdType ids[4];
259 double origin[3], x[3];
260
261 // define first corner
262 origin[0] = this->bX + i * this->hX;
263 origin[1] = this->bY + j * this->hY;
264 origin[2] = this->bZ + k * this->hZ;
265 ids[0] = pts->InsertNextPoint(origin);
266
267 if (face == 0) // x face
268 {
269 x[0] = origin[0];
270 x[1] = origin[1] + this->hY;
271 x[2] = origin[2];
272 ids[1] = pts->InsertNextPoint(x);
273
274 x[0] = origin[0];
275 x[1] = origin[1] + this->hY;
276 x[2] = origin[2] + this->hZ;
277 ids[2] = pts->InsertNextPoint(x);
278
279 x[0] = origin[0];
280 x[1] = origin[1];
281 x[2] = origin[2] + this->hZ;
282 ids[3] = pts->InsertNextPoint(x);
283 }
284
285 else if (face == 1) // y face
286 {
287 x[0] = origin[0] + this->hX;
288 x[1] = origin[1];
289 x[2] = origin[2];
290 ids[1] = pts->InsertNextPoint(x);
291
292 x[0] = origin[0] + this->hX;
293 x[1] = origin[1];
294 x[2] = origin[2] + this->hZ;
295 ids[2] = pts->InsertNextPoint(x);
296
297 x[0] = origin[0];
298 x[1] = origin[1];
299 x[2] = origin[2] + this->hZ;
300 ids[3] = pts->InsertNextPoint(x);
301 }
302
303 else // z face
304 {
305 x[0] = origin[0] + this->hX;
306 x[1] = origin[1];
307 x[2] = origin[2];
308 ids[1] = pts->InsertNextPoint(x);
309
310 x[0] = origin[0] + this->hX;
311 x[1] = origin[1] + this->hY;
312 x[2] = origin[2];
313 ids[2] = pts->InsertNextPoint(x);
314
315 x[0] = origin[0];
316 x[1] = origin[1] + this->hY;
317 x[2] = origin[2];
318 ids[3] = pts->InsertNextPoint(x);
319 }
320
321 polys->InsertNextCell(4, ids);
322 }
323
324 //------------------------------------------------------------------------------
325 // Calculate the distance between the point x to the bucket "nei".
326 //
327 // WARNING!!!!! Be very careful altering this routine. Simple changes to this
328 // routine can make is 25% slower!!!!
329 //
Distance2ToBucket(const double x[3],const int nei[3])330 double vtkBucketList::Distance2ToBucket(const double x[3], const int nei[3])
331 {
332 double bounds[6];
333
334 bounds[0] = nei[0] * this->hX + this->bX;
335 bounds[1] = (nei[0] + 1) * this->hX + this->bX;
336 bounds[2] = nei[1] * this->hY + this->bY;
337 bounds[3] = (nei[1] + 1) * this->hY + this->bY;
338 bounds[4] = nei[2] * this->hZ + this->bZ;
339 bounds[5] = (nei[2] + 1) * this->hZ + this->bZ;
340
341 return this->Distance2ToBounds(x, bounds);
342 }
343
344 //------------------------------------------------------------------------------
345 // Calculate the distance between the point x and the specified bounds
346 //
347 // WARNING!!!!! Be very careful altering this routine. Simple changes to this
348 // routine can make is 25% slower!!!!
Distance2ToBounds(const double x[3],const double bounds[6])349 double vtkBucketList::Distance2ToBounds(const double x[3], const double bounds[6])
350 {
351 double distance;
352 double deltas[3];
353
354 // Are we within the bounds?
355 if (x[0] >= bounds[0] && x[0] <= bounds[1] && x[1] >= bounds[2] && x[1] <= bounds[3] &&
356 x[2] >= bounds[4] && x[2] <= bounds[5])
357 {
358 return 0.0;
359 }
360
361 deltas[0] = deltas[1] = deltas[2] = 0.0;
362
363 // dx
364 //
365 if (x[0] < bounds[0])
366 {
367 deltas[0] = bounds[0] - x[0];
368 }
369 else if (x[0] > bounds[1])
370 {
371 deltas[0] = x[0] - bounds[1];
372 }
373
374 // dy
375 //
376 if (x[1] < bounds[2])
377 {
378 deltas[1] = bounds[2] - x[1];
379 }
380 else if (x[1] > bounds[3])
381 {
382 deltas[1] = x[1] - bounds[3];
383 }
384
385 // dz
386 //
387 if (x[2] < bounds[4])
388 {
389 deltas[2] = bounds[4] - x[2];
390 }
391 else if (x[2] > bounds[5])
392 {
393 deltas[2] = x[2] - bounds[5];
394 }
395
396 distance = vtkMath::Dot(deltas, deltas);
397 return distance;
398 }
399
400 //------------------------------------------------------------------------------
401 // The following tuple is what is sorted in the map. Note that it is templated
402 // because depending on the number of points / buckets to process we may want
403 // to use vtkIdType. Otherwise for performance reasons it's best to use an int
404 // (or other integral type). Typically sort() is 25-30% faster on smaller
405 // integral types, plus it takes a heck less memory (when vtkIdType is 64-bit
406 // and int is 32-bit).
407 template <typename TTuple>
408 struct LocatorTuple
409 {
410 TTuple PtId; // originating point id
411 TTuple Bucket; // i-j-k index into bucket space
412
413 // Operator< used to support the subsequent sort operation. There are two
414 // implementations, one gives a stable sort (points ordered by id within
415 // each bucket) and the other a little faster but less stable (in parallel
416 // sorting the order of sorted points in a bucket may vary).
417 // bool operator< (const LocatorTuple& tuple) const
418 // {return Bucket < tuple.Bucket;}
operator <LocatorTuple419 bool operator<(const LocatorTuple& tuple) const
420 {
421 if (Bucket < tuple.Bucket)
422 return true;
423 if (tuple.Bucket < Bucket)
424 return false;
425 if (PtId < tuple.PtId)
426 return true;
427 return false;
428 }
429 };
430
431 //------------------------------------------------------------------------------
432 // This templates class manages the creation of the static locator
433 // structures. It also implements the operator() functors which are supplied
434 // to vtkSMPTools for threaded processesing.
435 template <typename TIds>
436 struct BucketList : public vtkBucketList
437 {
438 // Okay the various ivars
439 LocatorTuple<TIds>* Map; // the map to be sorted
440 TIds* Offsets; // offsets for each bucket into the map
441
442 // Construction
BucketListBucketList443 BucketList(vtkStaticPointLocator* loc, vtkIdType numPts, int numBuckets)
444 : vtkBucketList(loc, numPts, numBuckets)
445 {
446 // one extra to simplify traversal
447 this->Map = new LocatorTuple<TIds>[numPts + 1];
448 this->Map[numPts].Bucket = numBuckets;
449 this->Offsets = new TIds[numBuckets + 1];
450 this->Offsets[numBuckets] = numPts;
451 }
452
453 // Release allocated memory
~BucketListBucketList454 ~BucketList() override
455 {
456 delete[] this->Map;
457 delete[] this->Offsets;
458 }
459
460 // The number of point ids in a bucket is determined by computing the
461 // difference between the offsets into the sorted points array.
GetNumberOfIdsBucketList462 vtkIdType GetNumberOfIds(vtkIdType bucketNum)
463 {
464 return (this->Offsets[bucketNum + 1] - this->Offsets[bucketNum]);
465 }
466
467 // Given a bucket number, return the point ids in that bucket.
GetIdsBucketList468 const LocatorTuple<TIds>* GetIds(vtkIdType bucketNum)
469 {
470 return this->Map + this->Offsets[bucketNum];
471 }
472
473 // Given a bucket number, return the point ids in that bucket.
GetIdsBucketList474 void GetIds(vtkIdType bucketNum, vtkIdList* bList)
475 {
476 const LocatorTuple<TIds>* ids = this->GetIds(bucketNum);
477 vtkIdType numIds = this->GetNumberOfIds(bucketNum);
478 bList->SetNumberOfIds(numIds);
479 for (int i = 0; i < numIds; i++)
480 {
481 bList->SetId(i, ids[i].PtId);
482 }
483 }
484
485 // Templated implementations of the locator
486 vtkIdType FindClosestPoint(const double x[3]);
487 vtkIdType FindClosestPointWithinRadius(
488 double radius, const double x[3], double inputDataLength, double& dist2);
489 void FindClosestNPoints(int N, const double x[3], vtkIdList* result);
490 void FindPointsWithinRadius(double R, const double x[3], vtkIdList* result);
491 int IntersectWithLine(double a0[3], double a1[3], double tol, double& t, double lineX[3],
492 double ptX[3], vtkIdType& ptId);
493 void MergePoints(double tol, vtkIdType* pointMap);
494 void GenerateRepresentation(int vtkNotUsed(level), vtkPolyData* pd);
495
496 // Internal methods
497 void GetOverlappingBuckets(
498 NeighborBuckets* buckets, const double x[3], const int ijk[3], double dist, int level);
499 void GetOverlappingBuckets(NeighborBuckets* buckets, const double x[3], double dist,
500 int prevMinLevel[3], int prevMaxLevel[3]);
501
502 // Implicit point representation, slower path
503 template <typename T>
504 struct MapDataSet
505 {
506 BucketList<T>* BList;
507 vtkDataSet* DataSet;
508
MapDataSetBucketList::MapDataSet509 MapDataSet(BucketList<T>* blist, vtkDataSet* ds)
510 : BList(blist)
511 , DataSet(ds)
512 {
513 }
514
operator ()BucketList::MapDataSet515 void operator()(vtkIdType ptId, vtkIdType end)
516 {
517 double p[3];
518 LocatorTuple<T>* t = this->BList->Map + ptId;
519 for (; ptId < end; ++ptId, ++t)
520 {
521 this->DataSet->GetPoint(ptId, p);
522 t->PtId = ptId;
523 t->Bucket = this->BList->GetBucketIndex(p);
524 } // for all points in this batch
525 }
526 };
527
528 // Explicit point representation (e.g., vtkPointSet), faster path
529 template <typename T, typename TPts>
530 struct MapPointsArray
531 {
532 BucketList<T>* BList;
533 const TPts* Points;
534
MapPointsArrayBucketList::MapPointsArray535 MapPointsArray(BucketList<T>* blist, const TPts* pts)
536 : BList(blist)
537 , Points(pts)
538 {
539 }
540
operator ()BucketList::MapPointsArray541 void operator()(vtkIdType ptId, vtkIdType end)
542 {
543 double p[3];
544 const TPts* x = this->Points + 3 * ptId;
545 LocatorTuple<T>* t = this->BList->Map + ptId;
546 for (; ptId < end; ++ptId, x += 3, ++t)
547 {
548 p[0] = static_cast<double>(x[0]);
549 p[1] = static_cast<double>(x[1]);
550 p[2] = static_cast<double>(x[2]);
551 t->PtId = ptId;
552 t->Bucket = this->BList->GetBucketIndex(p);
553 } // for all points in this batch
554 }
555 };
556
557 // A clever way to build offsets in parallel. Basically each thread builds
558 // offsets across a range of the sorted map. Recall that offsets are an
559 // integral value referring to the locations of the sorted points that
560 // reside in each bucket.
561 template <typename T>
562 struct MapOffsets
563 {
564 BucketList<T>* BList;
565 vtkIdType NumPts;
566 int NumBuckets;
567
MapOffsetsBucketList::MapOffsets568 MapOffsets(BucketList<T>* blist)
569 : BList(blist)
570 {
571 this->NumPts = this->BList->NumPts;
572 this->NumBuckets = this->BList->NumBuckets;
573 }
574
575 // Traverse sorted points (i.e., tuples) and update bucket offsets.
operator ()BucketList::MapOffsets576 void operator()(vtkIdType batch, vtkIdType batchEnd)
577 {
578 T* offsets = this->BList->Offsets;
579 const LocatorTuple<T>* curPt = this->BList->Map + batch * this->BList->BatchSize;
580 const LocatorTuple<T>* endBatchPt = this->BList->Map + batchEnd * this->BList->BatchSize;
581 const LocatorTuple<T>* endPt = this->BList->Map + this->NumPts;
582 const LocatorTuple<T>* prevPt;
583 endBatchPt = (endBatchPt > endPt ? endPt : endBatchPt);
584
585 // Special case at the very beginning of the mapped points array. If
586 // the first point is in bucket# N, then all buckets up and including
587 // N must refer to the first point.
588 if (curPt == this->BList->Map)
589 {
590 prevPt = this->BList->Map;
591 std::fill_n(offsets, curPt->Bucket + 1, 0); // point to the first points
592 } // at the very beginning of the map (sorted points array)
593
594 // We are entering this functor somewhere in the interior of the
595 // mapped points array. All we need to do is point to the entry
596 // position because we are interested only in prevPt->Bucket.
597 else
598 {
599 prevPt = curPt;
600 } // else in the middle of a batch
601
602 // Okay we have a starting point for a bucket run. Now we can begin
603 // filling in the offsets in this batch. A previous thread should
604 // have/will have completed the previous and subsequent runs outside
605 // of the [batch,batchEnd) range
606 for (curPt = prevPt; curPt < endBatchPt;)
607 {
608 for (; curPt->Bucket == prevPt->Bucket && curPt <= endBatchPt; ++curPt)
609 {
610 ; // advance
611 }
612 // Fill in any gaps in the offset array
613 std::fill_n(
614 offsets + prevPt->Bucket + 1, curPt->Bucket - prevPt->Bucket, curPt - this->BList->Map);
615 prevPt = curPt;
616 } // for all batches in this range
617 } // operator()
618 };
619
620 // Merge points that are pecisely coincident. Operates in parallel on
621 // locator buckets. Does not need to check neighbor buckets.
622 template <typename T>
623 struct MergePrecise
624 {
625 BucketList<T>* BList;
626 vtkDataSet* DataSet;
627 vtkIdType* MergeMap;
628
MergePreciseBucketList::MergePrecise629 MergePrecise(BucketList<T>* blist, vtkIdType* mergeMap)
630 : BList(blist)
631 , MergeMap(mergeMap)
632 {
633 this->DataSet = blist->DataSet;
634 }
635
operator ()BucketList::MergePrecise636 void operator()(vtkIdType bucket, vtkIdType endBucket)
637 {
638 BucketList<T>* bList = this->BList;
639 vtkIdType* mergeMap = this->MergeMap;
640 int i, j;
641 const LocatorTuple<TIds>* ids;
642 double p[3], p2[3];
643 vtkIdType ptId, ptId2, numIds;
644
645 for (; bucket < endBucket; ++bucket)
646 {
647 if ((numIds = bList->GetNumberOfIds(bucket)) > 0)
648 {
649 ids = bList->GetIds(bucket);
650 for (i = 0; i < numIds; i++)
651 {
652 ptId = ids[i].PtId;
653 if (mergeMap[ptId] < 0)
654 {
655 mergeMap[ptId] = ptId;
656 this->DataSet->GetPoint(ptId, p);
657 for (j = i + 1; j < numIds; j++)
658 {
659 ptId2 = ids[j].PtId;
660 if (mergeMap[ptId2] < 0)
661 {
662 this->DataSet->GetPoint(ptId2, p2);
663 if (p[0] == p2[0] && p[1] == p2[1] && p[2] == p2[2])
664 {
665 mergeMap[ptId2] = ptId;
666 }
667 }
668 }
669 } // if point not yet visited
670 }
671 }
672 }
673 }
674 };
675
676 // Merge points that are coincident within a tolerance. Operates in
677 // parallel on points. Needs to check neighbor buckets which slows it down
678 // considerably. Note that merging is one direction: larger ids are merged
679 // to lower.
680 template <typename T>
681 struct MergeClose
682 {
683 BucketList<T>* BList;
684 vtkDataSet* DataSet;
685 vtkIdType* MergeMap;
686 double Tol;
687
688 vtkSMPThreadLocalObject<vtkIdList> PIds;
689
MergeCloseBucketList::MergeClose690 MergeClose(BucketList<T>* blist, double tol, vtkIdType* mergeMap)
691 : BList(blist)
692 , MergeMap(mergeMap)
693 , Tol(tol)
694 {
695 this->DataSet = blist->DataSet;
696 }
697
698 // Just allocate a little bit of memory to get started.
InitializeBucketList::MergeClose699 void Initialize()
700 {
701 vtkIdList*& pIds = this->PIds.Local();
702 pIds->Allocate(128); // allocate some memory
703 }
704
operator ()BucketList::MergeClose705 void operator()(vtkIdType ptId, vtkIdType endPtId)
706 {
707 BucketList<T>* bList = this->BList;
708 vtkIdType* mergeMap = this->MergeMap;
709 int i;
710 double p[3];
711 vtkIdType nearId, numIds;
712 vtkIdList*& nearby = this->PIds.Local();
713
714 for (; ptId < endPtId; ++ptId)
715 {
716 if (mergeMap[ptId] < 0)
717 {
718 mergeMap[ptId] = ptId;
719 this->DataSet->GetPoint(ptId, p);
720 bList->FindPointsWithinRadius(this->Tol, p, nearby);
721 if ((numIds = nearby->GetNumberOfIds()) > 0)
722 {
723 for (i = 0; i < numIds; i++)
724 {
725 nearId = nearby->GetId(i);
726 if (ptId < nearId && (mergeMap[nearId] < 0 || ptId < mergeMap[nearId]))
727 {
728 mergeMap[nearId] = ptId;
729 }
730 }
731 }
732 } // if point not yet processed
733 } // for all points in this batch
734 }
735
ReduceBucketList::MergeClose736 void Reduce() {}
737 };
738
739 // Build the map and other structures to support locator operations
BuildLocatorBucketList740 void BuildLocator() override
741 {
742 // Place each point in a bucket
743 //
744 vtkPointSet* ps = vtkPointSet::SafeDownCast(this->DataSet);
745 int mapped = 0;
746 if (ps)
747 { // map points array: explicit points representation of float or double
748 int dataType = ps->GetPoints()->GetDataType();
749 void* pts = ps->GetPoints()->GetVoidPointer(0);
750 if (dataType == VTK_FLOAT)
751 {
752 MapPointsArray<TIds, float> mapper(this, static_cast<float*>(pts));
753 vtkSMPTools::For(0, this->NumPts, mapper);
754 mapped = 1;
755 }
756 else if (dataType == VTK_DOUBLE)
757 {
758 MapPointsArray<TIds, double> mapper(this, static_cast<double*>(pts));
759 vtkSMPTools::For(0, this->NumPts, mapper);
760 mapped = 1;
761 }
762 }
763
764 if (!mapped)
765 { // map dataset points: non-float points or implicit points representation
766 MapDataSet<TIds> mapper(this, this->DataSet);
767 vtkSMPTools::For(0, this->NumPts, mapper);
768 }
769
770 // Now gather the points into contiguous runs in buckets
771 //
772 vtkSMPTools::Sort(this->Map, this->Map + this->NumPts);
773
774 // Build the offsets into the Map. The offsets are the positions of
775 // each bucket into the sorted list. They mark the beginning of the
776 // list of points in each bucket. Amazingly, this can be done in
777 // parallel.
778 //
779 int numBatches = static_cast<int>(ceil(static_cast<double>(this->NumPts) / this->BatchSize));
780 MapOffsets<TIds> offMapper(this);
781 vtkSMPTools::For(0, numBatches, offMapper);
782 }
783 };
784
785 //------------------------------------------------------------------------------
786 // Given a position x, return the id of the point closest to it.
787 template <typename TIds>
FindClosestPoint(const double x[3])788 vtkIdType BucketList<TIds>::FindClosestPoint(const double x[3])
789 {
790 int i, j;
791 double minDist2;
792 double dist2 = VTK_DOUBLE_MAX;
793 double pt[3];
794 int closest, level;
795 vtkIdType ptId, cno, numIds;
796 int ijk[3], *nei;
797 NeighborBuckets buckets;
798 const LocatorTuple<TIds>* ids;
799
800 // Find bucket point is in.
801 //
802 this->GetBucketIndices(x, ijk);
803
804 // Need to search this bucket for the closest point. If there are no
805 // points in this bucket, search 1st level neighbors, and so on, until
806 // closest point found.
807 //
808 for (closest = (-1), minDist2 = VTK_DOUBLE_MAX, level = 0; (closest == -1) &&
809 (level < this->Divisions[0] || level < this->Divisions[1] || level < this->Divisions[2]);
810 level++)
811 {
812 this->GetBucketNeighbors(&buckets, ijk, this->Divisions, level);
813
814 for (i = 0; i < buckets.GetNumberOfNeighbors(); i++)
815 {
816 nei = buckets.GetPoint(i);
817 cno = nei[0] + nei[1] * this->xD + nei[2] * this->xyD;
818
819 if ((numIds = this->GetNumberOfIds(cno)) > 0)
820 {
821 ids = this->GetIds(cno);
822 for (j = 0; j < numIds; j++)
823 {
824 ptId = ids[j].PtId;
825 this->DataSet->GetPoint(ptId, pt);
826 if ((dist2 = vtkMath::Distance2BetweenPoints(x, pt)) < minDist2)
827 {
828 closest = ptId;
829 minDist2 = dist2;
830 }
831 }
832 }
833 }
834 }
835
836 //
837 // Because of the relative location of the points in the buckets, the
838 // point found previously may not be the closest point. We have to
839 // search those bucket neighbors that might also contain the point.
840 //
841 if (minDist2 > 0.0)
842 {
843 this->GetOverlappingBuckets(&buckets, x, ijk, sqrt(minDist2), 0);
844 for (i = 0; i < buckets.GetNumberOfNeighbors(); i++)
845 {
846 nei = buckets.GetPoint(i);
847 cno = nei[0] + nei[1] * this->xD + nei[2] * this->xyD;
848
849 if ((numIds = this->GetNumberOfIds(cno)) > 0)
850 {
851 ids = this->GetIds(cno);
852 for (j = 0; j < numIds; j++)
853 {
854 ptId = ids[j].PtId;
855 this->DataSet->GetPoint(ptId, pt);
856 if ((dist2 = vtkMath::Distance2BetweenPoints(x, pt)) < minDist2)
857 {
858 closest = ptId;
859 minDist2 = dist2;
860 }
861 } // for each point
862 } // if points in bucket
863 } // for each overlapping bucket
864 } // if not identical point
865
866 return closest;
867 }
868
869 //------------------------------------------------------------------------------
870 template <typename TIds>
FindClosestPointWithinRadius(double radius,const double x[3],double inputDataLength,double & dist2)871 vtkIdType BucketList<TIds>::FindClosestPointWithinRadius(
872 double radius, const double x[3], double inputDataLength, double& dist2)
873 {
874 int i, j;
875 double pt[3];
876 vtkIdType ptId, closest = -1;
877 int ijk[3], *nei;
878 double minDist2;
879
880 double refinedRadius, radius2, refinedRadius2;
881 double currentRadius;
882 double distance2ToDataBounds, maxDistance;
883 int ii, radiusLevels[3], radiusLevel, prevMinLevel[3], prevMaxLevel[3];
884 NeighborBuckets buckets;
885 const LocatorTuple<TIds>* ids;
886
887 // Initialize
888 dist2 = -1.0;
889 radius2 = radius * radius;
890 minDist2 = 1.01 * radius2; // something slightly bigger....
891
892 vtkDataArray* pointData = static_cast<vtkPointSet*>(this->DataSet)->GetPoints()->GetData();
893
894 // Find the bucket the point is in.
895 //
896 this->GetBucketIndices(x, ijk);
897
898 // Start by searching the bucket that the point is in.
899 //
900 vtkIdType numIds;
901 vtkIdType cno = ijk[0] + ijk[1] * this->xD + ijk[2] * this->xyD;
902 if ((numIds = this->GetNumberOfIds(cno)) > 0)
903 {
904 ids = this->GetIds(cno);
905 for (j = 0; j < numIds; j++)
906 {
907 ptId = ids[j].PtId;
908 pointData->GetTuple(ptId, pt);
909 if ((dist2 = vtkMath::Distance2BetweenPoints(x, pt)) < minDist2)
910 {
911 closest = ptId;
912 minDist2 = dist2;
913 }
914 }
915 }
916
917 // Now, search only those buckets that are within a radius. The radius used
918 // is the smaller of sqrt(minDist2) and the radius that is passed in. To avoid
919 // checking a large number of buckets unnecessarily, if the radius is
920 // larger than the dimensions of a bucket, we search outward using a
921 // simple heuristic of rings. This heuristic ends up collecting inner
922 // buckets multiple times, but this only happens in the case where these
923 // buckets are empty, so they are discarded quickly.
924 //
925 if (minDist2 < radius2)
926 {
927 refinedRadius = sqrt(minDist2);
928 refinedRadius2 = dist2;
929 }
930 else
931 {
932 refinedRadius = radius;
933 refinedRadius2 = radius2;
934 }
935
936 if (inputDataLength != 0.0)
937 {
938 distance2ToDataBounds = this->Distance2ToBounds(x, this->Bounds);
939 maxDistance = sqrt(distance2ToDataBounds) + inputDataLength;
940 if (refinedRadius > maxDistance)
941 {
942 refinedRadius = maxDistance;
943 refinedRadius2 = maxDistance * maxDistance;
944 }
945 }
946
947 for (i = 0; i < 3; i++)
948 {
949 radiusLevels[i] = static_cast<int>(refinedRadius / this->H[i]);
950 if (radiusLevels[i] > this->Divisions[i] / 2)
951 {
952 radiusLevels[i] = this->Divisions[i] / 2;
953 }
954 }
955
956 radiusLevel = radiusLevels[0];
957 radiusLevel = radiusLevels[1] > radiusLevel ? radiusLevels[1] : radiusLevel;
958 radiusLevel = radiusLevels[2] > radiusLevel ? radiusLevels[2] : radiusLevel;
959 if (radiusLevel == 0)
960 {
961 radiusLevel = 1;
962 }
963
964 // radius schedule increases the radius each iteration, this is currently
965 // implemented by decreasing ii by 1 each iteration. another alternative
966 // is to double the radius each iteration, i.e. ii = ii >> 1
967 // In practice, reducing ii by one has been found to be more efficient.
968 prevMinLevel[0] = prevMaxLevel[0] = ijk[0];
969 prevMinLevel[1] = prevMaxLevel[1] = ijk[1];
970 prevMinLevel[2] = prevMaxLevel[2] = ijk[2];
971 for (ii = radiusLevel; ii >= 1; ii--)
972 {
973 currentRadius = refinedRadius; // used in if at bottom of this for loop
974
975 // Build up a list of buckets that are arranged in rings
976 this->GetOverlappingBuckets(&buckets, x, refinedRadius / ii, prevMinLevel, prevMaxLevel);
977
978 for (i = 0; i < buckets.GetNumberOfNeighbors(); i++)
979 {
980 nei = buckets.GetPoint(i);
981 // do we still need to test this bucket?
982 if (this->Distance2ToBucket(x, nei) < refinedRadius2)
983 {
984 cno = nei[0] + nei[1] * this->xD + nei[2] * this->xyD;
985 if ((numIds = this->GetNumberOfIds(cno)) > 0)
986 {
987 ids = this->GetIds(cno);
988 for (j = 0; j < numIds; j++)
989 {
990 ptId = ids[j].PtId;
991 pointData->GetTuple(ptId, pt);
992 if ((dist2 = vtkMath::Distance2BetweenPoints(x, pt)) < minDist2)
993 {
994 closest = ptId;
995 minDist2 = dist2;
996 refinedRadius = sqrt(minDist2);
997 refinedRadius2 = minDist2;
998 }
999 } // for each pt in bucket
1000 } // if ids
1001 } // if bucket is within the current best distance
1002 } // for each overlapping bucket
1003
1004 // Don't want to check a smaller radius than we just checked so update
1005 // it appropriately
1006 if (refinedRadius < currentRadius && ii > 2) // always check ii==1
1007 {
1008 ii = static_cast<int>(static_cast<double>(ii) * (refinedRadius / currentRadius)) + 1;
1009 if (ii < 2)
1010 {
1011 ii = 2;
1012 }
1013 }
1014 } // for each radius in the radius schedule
1015
1016 if ((closest != -1) && (minDist2 <= radius2))
1017 {
1018 dist2 = minDist2;
1019 }
1020 else
1021 {
1022 closest = -1;
1023 }
1024
1025 return closest;
1026 }
1027
1028 namespace
1029 {
1030 //------------------------------------------------------------------------------
1031 // Obtaining closest points requires sorting nearby points
1032 struct IdTuple
1033 {
1034 vtkIdType PtId;
1035 double Dist2;
1036
operator <__anon7e6903b50111::IdTuple1037 bool operator<(const IdTuple& tuple) const { return Dist2 < tuple.Dist2; }
1038 };
1039 }
1040
1041 //------------------------------------------------------------------------------
1042 template <typename TIds>
FindClosestNPoints(int N,const double x[3],vtkIdList * result)1043 void BucketList<TIds>::FindClosestNPoints(int N, const double x[3], vtkIdList* result)
1044 {
1045 int i, j;
1046 double dist2;
1047 double pt[3];
1048 int level;
1049 vtkIdType ptId, cno, numIds;
1050 int ijk[3], *nei;
1051 NeighborBuckets buckets;
1052 const LocatorTuple<TIds>* ids;
1053
1054 // Clear out any previous results
1055 result->Reset();
1056
1057 // Find the bucket the point is in.
1058 //
1059 this->GetBucketIndices(x, ijk);
1060
1061 // There are two steps, first a simple expanding wave of buckets until
1062 // we have enough points. Then a refinement to make sure we have the
1063 // N closest points.
1064 level = 0;
1065 double maxDistance = 0.0;
1066 int currentCount = 0;
1067 std::vector<IdTuple> res(N);
1068
1069 this->GetBucketNeighbors(&buckets, ijk, this->Divisions, level);
1070 while (buckets.GetNumberOfNeighbors() && currentCount < N)
1071 {
1072 for (i = 0; i < buckets.GetNumberOfNeighbors(); i++)
1073 {
1074 nei = buckets.GetPoint(i);
1075 cno = nei[0] + nei[1] * this->xD + nei[2] * this->xyD;
1076
1077 if ((numIds = this->GetNumberOfIds(cno)) > 0)
1078 {
1079 ids = this->GetIds(cno);
1080 for (j = 0; j < numIds; j++)
1081 {
1082 ptId = ids[j].PtId;
1083 this->DataSet->GetPoint(ptId, pt);
1084 dist2 = vtkMath::Distance2BetweenPoints(x, pt);
1085 if (currentCount < N)
1086 {
1087 res[currentCount].Dist2 = dist2;
1088 res[currentCount].PtId = ptId;
1089 if (dist2 > maxDistance)
1090 {
1091 maxDistance = dist2;
1092 }
1093 currentCount++;
1094 if (currentCount == N)
1095 {
1096 std::sort(res.begin(), res.begin() + currentCount);
1097 }
1098 }
1099 else if (dist2 < maxDistance)
1100 {
1101 res[N - 1].Dist2 = dist2;
1102 res[N - 1].PtId = ptId;
1103 std::sort(res.begin(), res.begin() + N);
1104 maxDistance = res[N - 1].Dist2;
1105 }
1106 }
1107 }
1108 }
1109 level++;
1110 this->GetBucketNeighbors(&buckets, ijk, this->Divisions, level);
1111 }
1112
1113 // do a sort
1114 std::sort(res.begin(), res.begin() + currentCount);
1115
1116 // Now do the refinement
1117 this->GetOverlappingBuckets(&buckets, x, ijk, sqrt(maxDistance), level - 1);
1118
1119 for (i = 0; i < buckets.GetNumberOfNeighbors(); i++)
1120 {
1121 nei = buckets.GetPoint(i);
1122 cno = nei[0] + nei[1] * this->xD + nei[2] * this->xyD;
1123
1124 if ((numIds = this->GetNumberOfIds(cno)) > 0)
1125 {
1126 ids = this->GetIds(cno);
1127 for (j = 0; j < numIds; j++)
1128 {
1129 ptId = ids[j].PtId;
1130 this->DataSet->GetPoint(ptId, pt);
1131 dist2 = vtkMath::Distance2BetweenPoints(x, pt);
1132 if (dist2 < maxDistance)
1133 {
1134 res[N - 1].Dist2 = dist2;
1135 res[N - 1].PtId = ptId;
1136 std::sort(res.begin(), res.begin() + N);
1137 maxDistance = res[N - 1].Dist2;
1138 }
1139 }
1140 }
1141 }
1142
1143 // Fill in the IdList
1144 result->SetNumberOfIds(currentCount);
1145 for (i = 0; i < currentCount; i++)
1146 {
1147 result->SetId(i, res[i].PtId);
1148 }
1149 }
1150
1151 //------------------------------------------------------------------------------
1152 // The Radius defines a block of buckets which the sphere of radis R may
1153 // touch.
1154 template <typename TIds>
FindPointsWithinRadius(double R,const double x[3],vtkIdList * result)1155 void BucketList<TIds>::FindPointsWithinRadius(double R, const double x[3], vtkIdList* result)
1156 {
1157 double dist2;
1158 double pt[3];
1159 vtkIdType ptId, cno, numIds;
1160 double R2 = R * R;
1161 const LocatorTuple<TIds>* ids;
1162 double xMin[3], xMax[3];
1163 int i, j, k, ii, jOffset, kOffset, ijkMin[3], ijkMax[3];
1164
1165 // Determine the range of indices in each direction based on radius R
1166 xMin[0] = x[0] - R;
1167 xMin[1] = x[1] - R;
1168 xMin[2] = x[2] - R;
1169 xMax[0] = x[0] + R;
1170 xMax[1] = x[1] + R;
1171 xMax[2] = x[2] + R;
1172
1173 // Find the footprint in the locator
1174 this->GetBucketIndices(xMin, ijkMin);
1175 this->GetBucketIndices(xMax, ijkMax);
1176
1177 // Clear out previous results
1178 result->Reset();
1179
1180 // Add points within footprint and radius
1181 for (k = ijkMin[2]; k <= ijkMax[2]; ++k)
1182 {
1183 kOffset = k * this->xyD;
1184 for (j = ijkMin[1]; j <= ijkMax[1]; ++j)
1185 {
1186 jOffset = j * this->xD;
1187 for (i = ijkMin[0]; i <= ijkMax[0]; ++i)
1188 {
1189 cno = i + jOffset + kOffset;
1190
1191 if ((numIds = this->GetNumberOfIds(cno)) > 0)
1192 {
1193 ids = this->GetIds(cno);
1194 for (ii = 0; ii < numIds; ii++)
1195 {
1196 ptId = ids[ii].PtId;
1197 this->DataSet->GetPoint(ptId, pt);
1198 dist2 = vtkMath::Distance2BetweenPoints(x, pt);
1199 if (dist2 <= R2)
1200 {
1201 result->InsertNextId(ptId);
1202 }
1203 } // for all points in bucket
1204 } // if points in bucket
1205 } // i-footprint
1206 } // j-footprint
1207 } // k-footprint
1208 }
1209
1210 //------------------------------------------------------------------------------
1211 // Find the point within tol of the finite line, and closest to the starting
1212 // point of the line (i.e., min parametric coordinate t).
1213 //
1214 // Note that we have to traverse more than just the buckets (aka bins)
1215 // containing the line since the closest point could be in a neighboring
1216 // bin. To keep the code simple here's the straightforward approach used in
1217 // the code below. Imagine tracing a sphere of radius tol along the finite
1218 // line, and processing all bins (and of course the points in the bins) which
1219 // intersect the sphere. We use a typical ray tracing approach (see
1220 // vtkStaticCellLocator for references) and update the current voxels/bins at
1221 // boundaries, including intersecting the sphere with neighboring bins. Since
1222 // this simple approach may visit bins multiple times, we keep an array that
1223 // marks whether the bin has been visited previously and skip it if we have.
1224 template <typename TIds>
IntersectWithLine(double a0[3],double a1[3],double tol,double & t,double lineX[3],double ptX[3],vtkIdType & ptId)1225 int BucketList<TIds>::IntersectWithLine(double a0[3], double a1[3], double tol, double& t,
1226 double lineX[3], double ptX[3], vtkIdType& ptId)
1227 {
1228 double* bounds = this->Bounds;
1229 int* ndivs = this->Divisions;
1230 vtkIdType prod = ndivs[0] * ndivs[1];
1231 double* h = this->H;
1232 TIds ii, numPtsInBin;
1233 double x[3], xl[3], rayDir[3], xmin[3], xmax[3];
1234 vtkMath::Subtract(a1, a0, rayDir);
1235 double curPos[3], curT, tHit, tMin = VTK_FLOAT_MAX;
1236 int i, j, k, enterExitCount;
1237 int ijk[3], ijkMin[3], ijkMax[3];
1238 vtkIdType idx, pId, bestPtId = (-1);
1239 double step[3], next[3], tMax[3], tDelta[3];
1240 double tol2 = tol * tol;
1241
1242 // Make sure the bounding box of the locator is hit
1243 if (vtkBox::IntersectBox(bounds, a0, rayDir, curPos, curT))
1244 {
1245 // Initialize intersection query array if necessary. This is done
1246 // locally to ensure thread safety.
1247 std::vector<unsigned char> bucketHasBeenVisited(this->NumBuckets, 0);
1248
1249 // Get the i-j-k point of intersection and bin index. This is
1250 // clamped to the boundary of the locator.
1251 this->GetBucketIndices(curPos, ijk);
1252
1253 // Set up some parameters for traversing through bins
1254 step[0] = (rayDir[0] >= 0.0) ? 1.0 : -1.0;
1255 step[1] = (rayDir[1] >= 0.0) ? 1.0 : -1.0;
1256 step[2] = (rayDir[2] >= 0.0) ? 1.0 : -1.0;
1257
1258 // If the ray is going in the negative direction, then the next voxel boundary
1259 // is on the "-" direction so we stay in the current voxel.
1260 next[0] = bounds[0] + h[0] * (rayDir[0] >= 0.0 ? (ijk[0] + step[0]) : ijk[0]);
1261 next[1] = bounds[2] + h[1] * (rayDir[1] >= 0.0 ? (ijk[1] + step[1]) : ijk[1]);
1262 next[2] = bounds[4] + h[2] * (rayDir[2] >= 0.0 ? (ijk[2] + step[2]) : ijk[2]);
1263
1264 tMax[0] = (rayDir[0] != 0.0) ? (next[0] - curPos[0]) / rayDir[0] : VTK_FLOAT_MAX;
1265 tMax[1] = (rayDir[1] != 0.0) ? (next[1] - curPos[1]) / rayDir[1] : VTK_FLOAT_MAX;
1266 tMax[2] = (rayDir[2] != 0.0) ? (next[2] - curPos[2]) / rayDir[2] : VTK_FLOAT_MAX;
1267
1268 tDelta[0] = (rayDir[0] != 0.0) ? (h[0] / rayDir[0]) * step[0] : VTK_FLOAT_MAX;
1269 tDelta[1] = (rayDir[1] != 0.0) ? (h[1] / rayDir[1]) * step[1] : VTK_FLOAT_MAX;
1270 tDelta[2] = (rayDir[2] != 0.0) ? (h[2] / rayDir[2]) * step[2] : VTK_FLOAT_MAX;
1271
1272 // Process current position including the bins in the sphere
1273 // footprint. Note there is a rare pathological case where the footprint
1274 // on voxel exit must also be considered.
1275 for (bestPtId = (-1), enterExitCount = 0; bestPtId < 0 || enterExitCount < 2;)
1276 {
1277 // Get the "footprint" of bins containing the sphere defined by the
1278 // current position and a radius of tol.
1279 xmin[0] = curPos[0] - tol;
1280 xmin[1] = curPos[1] - tol;
1281 xmin[2] = curPos[2] - tol;
1282 xmax[0] = curPos[0] + tol;
1283 xmax[1] = curPos[1] + tol;
1284 xmax[2] = curPos[2] + tol;
1285 this->GetBucketIndices(xmin, ijkMin);
1286 this->GetBucketIndices(xmax, ijkMax);
1287
1288 // Start walking through the bins, find the best point of
1289 // intersection. Note that the ray may not penetrate all of the way
1290 // through the locator so may terminate when (t > 1.0).
1291 for (k = ijkMin[2]; k <= ijkMax[2]; ++k)
1292 {
1293 for (j = ijkMin[1]; j <= ijkMax[1]; ++j)
1294 {
1295 for (i = ijkMin[0]; i <= ijkMax[0]; ++i)
1296 {
1297 // Current bin index
1298 idx = i + j * ndivs[0] + k * prod;
1299
1300 if (!bucketHasBeenVisited[idx])
1301 {
1302 bucketHasBeenVisited[idx] = 1;
1303 if ((numPtsInBin = this->GetNumberOfIds(idx)) > 0) // there are some points here
1304 {
1305 const LocatorTuple<TIds>* ptIds = this->GetIds(idx);
1306 for (ii = 0; ii < numPtsInBin; ii++)
1307 {
1308 pId = ptIds[ii].PtId;
1309 this->DataSet->GetPoint(pId, x);
1310 if (vtkLine::DistanceToLine(x, a0, a1, tHit, xl) <= tol2 && t < tMin)
1311 {
1312 tMin = t;
1313 bestPtId = pId;
1314 } // point is within tolerance and closer
1315 } // over all points in bin
1316 } // if points in bin
1317 } // bucket not visited
1318 } // i bins
1319 } // j bins
1320 } // k bins
1321
1322 // Make sure to evaluate exit footprint as well. Must evaluate entrance
1323 // and exit of current voxel.
1324 if (bestPtId >= 0)
1325 {
1326 enterExitCount++;
1327 }
1328
1329 // Advance to next voxel / bin
1330 if (tMax[0] < tMax[1])
1331 {
1332 if (tMax[0] < tMax[2])
1333 {
1334 ijk[0] += static_cast<int>(step[0]);
1335 tMax[0] += tDelta[0];
1336 curT = tMax[0];
1337 }
1338 else
1339 {
1340 ijk[2] += static_cast<int>(step[2]);
1341 tMax[2] += tDelta[2];
1342 curT = tMax[2];
1343 }
1344 }
1345 else
1346 {
1347 if (tMax[1] < tMax[2])
1348 {
1349 ijk[1] += static_cast<int>(step[1]);
1350 tMax[1] += tDelta[1];
1351 curT = tMax[1];
1352 }
1353 else
1354 {
1355 ijk[2] += static_cast<int>(step[2]);
1356 tMax[2] += tDelta[2];
1357 curT = tMax[2];
1358 }
1359 }
1360
1361 // Check exit conditions
1362 if (curT > 1.0 || ijk[0] < 0 || ijk[0] >= ndivs[0] || ijk[1] < 0 || ijk[1] >= ndivs[1] ||
1363 ijk[2] < 0 || ijk[2] >= ndivs[2])
1364 {
1365 break;
1366 }
1367 else
1368 {
1369 curPos[0] = a0[0] + curT * rayDir[0];
1370 curPos[1] = a0[1] + curT * rayDir[1];
1371 curPos[2] = a0[2] + curT * rayDir[2];
1372 }
1373
1374 } // for looking for valid intersected point
1375 } // if (vtkBox::IntersectBox(...))
1376
1377 // If a point has been intersected, recover the information and return.
1378 // This information could be cached....
1379 if (bestPtId >= 0)
1380 {
1381 // update the return information
1382 ptId = bestPtId;
1383 this->DataSet->GetPoint(ptId, ptX);
1384 vtkLine::DistanceToLine(ptX, a0, a1, t, lineX);
1385
1386 return 1;
1387 }
1388
1389 return 0;
1390 }
1391
1392 //------------------------------------------------------------------------------
1393 // Merge points based on tolerance. Return a point map. There are two
1394 // separate paths: when the tolerance is precisely 0.0, and when tol >
1395 // 0.0. Both are executed in parallel, although the second uses a
1396 // checkerboard approach to avoid write collisions.
1397 template <typename TIds>
MergePoints(double tol,vtkIdType * mergeMap)1398 void BucketList<TIds>::MergePoints(double tol, vtkIdType* mergeMap)
1399 {
1400 // First mark all points as uninitialized
1401 std::fill_n(mergeMap, this->NumPts, (-1));
1402
1403 // If tol=0, then just process points bucket by bucket. Don't have to worry
1404 // about points in other buckets.
1405 if (tol <= 0.0)
1406 {
1407 MergePrecise<TIds> merge(this, mergeMap);
1408 vtkSMPTools::For(0, this->NumBuckets, merge);
1409 }
1410
1411 // Merge within a tolerance. This is a greedy algorithm that can give
1412 // weird results since exactly which points to merge with is not an
1413 // obvious answer (without doing fancy clustering etc).
1414 else
1415 {
1416 MergeClose<TIds> merge(this, tol, mergeMap);
1417 vtkSMPTools::For(0, this->NumPts, merge);
1418 }
1419 }
1420
1421 //------------------------------------------------------------------------------
1422 // Internal method to find those buckets that are within distance specified
1423 // only those buckets outside of level radiuses of ijk are returned
1424 template <typename TIds>
GetOverlappingBuckets(NeighborBuckets * buckets,const double x[3],const int ijk[3],double dist,int level)1425 void BucketList<TIds>::GetOverlappingBuckets(
1426 NeighborBuckets* buckets, const double x[3], const int ijk[3], double dist, int level)
1427 {
1428 int i, j, k, nei[3], minLevel[3], maxLevel[3];
1429 double xMin[3], xMax[3];
1430
1431 // Initialize
1432 buckets->Reset();
1433
1434 // Determine the range of indices in each direction
1435 xMin[0] = x[0] - dist;
1436 xMin[1] = x[1] - dist;
1437 xMin[2] = x[2] - dist;
1438 xMax[0] = x[0] + dist;
1439 xMax[1] = x[1] + dist;
1440 xMax[2] = x[2] + dist;
1441
1442 this->GetBucketIndices(xMin, minLevel);
1443 this->GetBucketIndices(xMax, maxLevel);
1444
1445 for (i = minLevel[0]; i <= maxLevel[0]; i++)
1446 {
1447 for (j = minLevel[1]; j <= maxLevel[1]; j++)
1448 {
1449 for (k = minLevel[2]; k <= maxLevel[2]; k++)
1450 {
1451 if (i < (ijk[0] - level) || i > (ijk[0] + level) || j < (ijk[1] - level) ||
1452 j > (ijk[1] + level) || k < (ijk[2] - level) || k > (ijk[2] + level))
1453 {
1454 nei[0] = i;
1455 nei[1] = j;
1456 nei[2] = k;
1457 buckets->InsertNextBucket(nei);
1458 }
1459 }
1460 }
1461 }
1462 }
1463
1464 //------------------------------------------------------------------------------
1465 // Internal method to find those buckets that are within distance specified
1466 // only those buckets outside of level radiuses of ijk are returned
1467 template <typename TIds>
GetOverlappingBuckets(NeighborBuckets * buckets,const double x[3],double dist,int prevMinLevel[3],int prevMaxLevel[3])1468 void BucketList<TIds>::GetOverlappingBuckets(NeighborBuckets* buckets, const double x[3],
1469 double dist, int prevMinLevel[3], int prevMaxLevel[3])
1470 {
1471 int i, j, k, nei[3], minLevel[3], maxLevel[3];
1472 int kFactor, jFactor;
1473 int jkSkipFlag, kSkipFlag;
1474 double xMin[3], xMax[3];
1475
1476 // Initialize
1477 buckets->Reset();
1478
1479 // Determine the range of indices in each direction
1480 xMin[0] = x[0] - dist;
1481 xMin[1] = x[1] - dist;
1482 xMin[2] = x[2] - dist;
1483 xMax[0] = x[0] + dist;
1484 xMax[1] = x[1] + dist;
1485 xMax[2] = x[2] + dist;
1486
1487 this->GetBucketIndices(xMin, minLevel);
1488 this->GetBucketIndices(xMax, maxLevel);
1489
1490 if (minLevel[0] == prevMinLevel[0] && maxLevel[0] == prevMaxLevel[0] &&
1491 minLevel[1] == prevMinLevel[1] && maxLevel[1] == prevMaxLevel[1] &&
1492 minLevel[2] == prevMinLevel[2] && maxLevel[2] == prevMaxLevel[2])
1493 {
1494 return;
1495 }
1496
1497 for (k = minLevel[2]; k <= maxLevel[2]; k++)
1498 {
1499 kFactor = k * this->xyD;
1500 if (k >= prevMinLevel[2] && k <= prevMaxLevel[2])
1501 {
1502 kSkipFlag = 1;
1503 }
1504 else
1505 {
1506 kSkipFlag = 0;
1507 }
1508 for (j = minLevel[1]; j <= maxLevel[1]; j++)
1509 {
1510 if (kSkipFlag && j >= prevMinLevel[1] && j <= prevMaxLevel[1])
1511 {
1512 jkSkipFlag = 1;
1513 }
1514 else
1515 {
1516 jkSkipFlag = 0;
1517 }
1518 jFactor = j * this->xD;
1519 for (i = minLevel[0]; i <= maxLevel[0]; i++)
1520 {
1521 if (jkSkipFlag && i == prevMinLevel[0])
1522 {
1523 i = prevMaxLevel[0];
1524 continue;
1525 }
1526 // if this bucket has any cells, add it to the list
1527 if (this->GetNumberOfIds(i + jFactor + kFactor) > 0)
1528 {
1529 nei[0] = i;
1530 nei[1] = j;
1531 nei[2] = k;
1532 buckets->InsertNextBucket(nei);
1533 }
1534 }
1535 }
1536 }
1537
1538 prevMinLevel[0] = minLevel[0];
1539 prevMinLevel[1] = minLevel[1];
1540 prevMinLevel[2] = minLevel[2];
1541 prevMaxLevel[0] = maxLevel[0];
1542 prevMaxLevel[1] = maxLevel[1];
1543 prevMaxLevel[2] = maxLevel[2];
1544 }
1545
1546 //------------------------------------------------------------------------------
1547 // Build polygonal representation of locator. Create faces that separate
1548 // inside/outside buckets, or separate inside/boundary of locator.
1549 template <typename TIds>
GenerateRepresentation(int vtkNotUsed (level),vtkPolyData * pd)1550 void BucketList<TIds>::GenerateRepresentation(int vtkNotUsed(level), vtkPolyData* pd)
1551 {
1552 vtkPoints* pts;
1553 vtkCellArray* polys;
1554 int ii, i, j, k, idx, offset[3], minusOffset[3], inside, sliceSize;
1555
1556 pts = vtkPoints::New();
1557 pts->Allocate(5000);
1558 polys = vtkCellArray::New();
1559 polys->AllocateEstimate(2048, 3);
1560
1561 // loop over all buckets, creating appropriate faces
1562 sliceSize = this->Divisions[0] * this->Divisions[1];
1563 for (k = 0; k < this->Divisions[2]; k++)
1564 {
1565 offset[2] = k * sliceSize;
1566 minusOffset[2] = (k - 1) * sliceSize;
1567 for (j = 0; j < this->Divisions[1]; j++)
1568 {
1569 offset[1] = j * this->Divisions[0];
1570 minusOffset[1] = (j - 1) * this->Divisions[0];
1571 for (i = 0; i < this->Divisions[0]; i++)
1572 {
1573 offset[0] = i;
1574 minusOffset[0] = i - 1;
1575 idx = offset[0] + offset[1] + offset[2];
1576 if (this->GetNumberOfIds(idx) > 0)
1577 {
1578 inside = 0;
1579 }
1580 else
1581 {
1582 inside = 1;
1583 }
1584
1585 // check "negative" neighbors
1586 for (ii = 0; ii < 3; ii++)
1587 {
1588 if (minusOffset[ii] < 0)
1589 {
1590 if (inside)
1591 {
1592 this->GenerateFace(ii, i, j, k, pts, polys);
1593 }
1594 }
1595 else
1596 {
1597 if (ii == 0)
1598 {
1599 idx = minusOffset[0] + offset[1] + offset[2];
1600 }
1601 else if (ii == 1)
1602 {
1603 idx = offset[0] + minusOffset[1] + offset[2];
1604 }
1605 else
1606 {
1607 idx = offset[0] + offset[1] + minusOffset[2];
1608 }
1609
1610 if ((this->GetNumberOfIds(idx) > 0 && inside) ||
1611 (this->GetNumberOfIds(idx) > 0 && !inside))
1612 {
1613 this->GenerateFace(ii, i, j, k, pts, polys);
1614 }
1615 }
1616 // those buckets on "positive" boundaries can generate faces specially
1617 if ((i + 1) >= this->Divisions[0] && inside)
1618 {
1619 this->GenerateFace(0, i + 1, j, k, pts, polys);
1620 }
1621 if ((j + 1) >= this->Divisions[1] && inside)
1622 {
1623 this->GenerateFace(1, i, j + 1, k, pts, polys);
1624 }
1625 if ((k + 1) >= this->Divisions[2] && inside)
1626 {
1627 this->GenerateFace(2, i, j, k + 1, pts, polys);
1628 }
1629
1630 } // over negative faces
1631 } // over i divisions
1632 } // over j divisions
1633 } // over k divisions
1634
1635 pd->SetPoints(pts);
1636 pts->Delete();
1637 pd->SetPolys(polys);
1638 polys->Delete();
1639 pd->Squeeze();
1640 }
1641
1642 //------------------------------------------------------------------------------
1643 // Here is the VTK class proper. It's implemented with the templated
1644 // BucketList class.
1645
1646 //------------------------------------------------------------------------------
1647 // Construct with automatic computation of divisions, averaging
1648 // 5 points per bucket.
vtkStaticPointLocator()1649 vtkStaticPointLocator::vtkStaticPointLocator()
1650 {
1651 this->NumberOfPointsPerBucket = 1;
1652 this->Divisions[0] = this->Divisions[1] = this->Divisions[2] = 50;
1653 this->H[0] = this->H[1] = this->H[2] = 0.0;
1654 this->Buckets = nullptr;
1655 this->MaxNumberOfBuckets = VTK_INT_MAX;
1656 this->LargeIds = false;
1657 }
1658
1659 //------------------------------------------------------------------------------
~vtkStaticPointLocator()1660 vtkStaticPointLocator::~vtkStaticPointLocator()
1661 {
1662 this->FreeSearchStructure();
1663 }
1664
1665 //------------------------------------------------------------------------------
Initialize()1666 void vtkStaticPointLocator::Initialize()
1667 {
1668 this->FreeSearchStructure();
1669 }
1670
1671 //------------------------------------------------------------------------------
FreeSearchStructure()1672 void vtkStaticPointLocator::FreeSearchStructure()
1673 {
1674 if (this->Buckets)
1675 {
1676 delete this->Buckets;
1677 this->Buckets = nullptr;
1678 }
1679 }
1680
1681 //------------------------------------------------------------------------------
1682 // Method to form subdivision of space based on the points provided and
1683 // subject to the constraints of levels and NumberOfPointsPerBucket.
1684 // The result is directly addressable and of uniform subdivision.
1685 //
BuildLocator()1686 void vtkStaticPointLocator::BuildLocator()
1687 {
1688 int ndivs[3];
1689 int i;
1690 vtkIdType numPts;
1691
1692 if ((this->Buckets != nullptr) && (this->BuildTime > this->MTime) &&
1693 (this->BuildTime > this->DataSet->GetMTime()))
1694 {
1695 return;
1696 }
1697
1698 vtkDebugMacro(<< "Hashing points...");
1699 this->Level = 1; // only single lowest level - from superclass
1700
1701 if (!this->DataSet || (numPts = this->DataSet->GetNumberOfPoints()) < 1)
1702 {
1703 vtkErrorMacro(<< "No points to locate");
1704 return;
1705 }
1706
1707 // Make sure the appropriate data is available
1708 //
1709 if (this->Buckets)
1710 {
1711 this->FreeSearchStructure();
1712 }
1713
1714 // Size the root bucket. Initialize bucket data structure, compute
1715 // level and divisions. The GetBounds() method below can be very slow;
1716 // hopefully it is cached or otherwise accelerated.
1717 //
1718 const double* bounds = this->DataSet->GetBounds();
1719 vtkIdType numBuckets = static_cast<vtkIdType>(
1720 static_cast<double>(numPts) / static_cast<double>(this->NumberOfPointsPerBucket));
1721 numBuckets = (numBuckets > this->MaxNumberOfBuckets ? this->MaxNumberOfBuckets : numBuckets);
1722
1723 vtkBoundingBox bbox(bounds);
1724 if (this->Automatic)
1725 {
1726 bbox.ComputeDivisions(numBuckets, this->Bounds, ndivs);
1727 }
1728 else
1729 {
1730 bbox.Inflate(); // make sure non-zero volume
1731 bbox.GetBounds(this->Bounds);
1732 for (i = 0; i < 3; i++)
1733 {
1734 ndivs[i] = (this->Divisions[i] < 1 ? 1 : this->Divisions[i]);
1735 }
1736 }
1737
1738 this->Divisions[0] = ndivs[0];
1739 this->Divisions[1] = ndivs[1];
1740 this->Divisions[2] = ndivs[2];
1741 this->NumberOfBuckets = numBuckets = static_cast<vtkIdType>(ndivs[0]) *
1742 static_cast<vtkIdType>(ndivs[1]) * static_cast<vtkIdType>(ndivs[2]);
1743
1744 // Compute width of bucket in three directions
1745 //
1746 for (i = 0; i < 3; i++)
1747 {
1748 this->H[i] = (this->Bounds[2 * i + 1] - this->Bounds[2 * i]) / static_cast<double>(ndivs[i]);
1749 }
1750
1751 // Instantiate the locator. The type is related to the maximum point id.
1752 // This is done for performance (e.g., the sort is faster) and significant
1753 // memory savings.
1754 //
1755 if (numPts >= VTK_INT_MAX || numBuckets >= VTK_INT_MAX)
1756 {
1757 this->LargeIds = true;
1758 this->Buckets = new BucketList<vtkIdType>(this, numPts, numBuckets);
1759 }
1760 else
1761 {
1762 this->LargeIds = false;
1763 this->Buckets = new BucketList<int>(this, numPts, numBuckets);
1764 }
1765
1766 // Actually construct the locator
1767 this->Buckets->BuildLocator();
1768
1769 this->BuildTime.Modified();
1770 }
1771
1772 //------------------------------------------------------------------------------
1773 // Method to form subdivision of space based on the points provided and
1774 // subject to the constraints of levels and NumberOfPointsPerBucket.
1775 // The result is directly addressable and of uniform subdivision.
1776 //
BuildLocator(const double * inBounds)1777 void vtkStaticPointLocator::BuildLocator(const double* inBounds)
1778 {
1779 int ndivs[3];
1780 int i;
1781 vtkIdType numPts;
1782
1783 if ((this->Buckets != nullptr) && (this->BuildTime > this->MTime) &&
1784 (this->BuildTime > this->DataSet->GetMTime()))
1785 {
1786 return;
1787 }
1788
1789 vtkDebugMacro(<< "Hashing points...");
1790 this->Level = 1; // only single lowest level - from superclass
1791
1792 if (!this->DataSet || (numPts = this->DataSet->GetNumberOfPoints()) < 1)
1793 {
1794 vtkErrorMacro(<< "No points to locate");
1795 return;
1796 }
1797
1798 // Make sure the appropriate data is available
1799 //
1800 if (this->Buckets)
1801 {
1802 this->FreeSearchStructure();
1803 }
1804
1805 // Size the root bucket. Initialize bucket data structure, compute
1806 // level and divisions. The GetBounds() method below can be very slow;
1807 // hopefully it is cached or otherwise accelerated.
1808 //
1809 const double* bounds = (inBounds == nullptr ? this->DataSet->GetBounds() : inBounds);
1810 vtkIdType numBuckets = static_cast<vtkIdType>(
1811 static_cast<double>(numPts) / static_cast<double>(this->NumberOfPointsPerBucket));
1812 numBuckets = (numBuckets > this->MaxNumberOfBuckets ? this->MaxNumberOfBuckets : numBuckets);
1813
1814 vtkBoundingBox bbox(bounds);
1815 if (this->Automatic)
1816 {
1817 bbox.ComputeDivisions(numBuckets, this->Bounds, ndivs);
1818 }
1819 else
1820 {
1821 bbox.Inflate(); // make sure non-zero volume
1822 bbox.GetBounds(this->Bounds);
1823 for (i = 0; i < 3; i++)
1824 {
1825 ndivs[i] = (this->Divisions[i] < 1 ? 1 : this->Divisions[i]);
1826 }
1827 }
1828
1829 this->Divisions[0] = ndivs[0];
1830 this->Divisions[1] = ndivs[1];
1831 this->Divisions[2] = ndivs[2];
1832 this->NumberOfBuckets = numBuckets = static_cast<vtkIdType>(ndivs[0]) *
1833 static_cast<vtkIdType>(ndivs[1]) * static_cast<vtkIdType>(ndivs[2]);
1834
1835 // Compute width of bucket in three directions
1836 //
1837 for (i = 0; i < 3; i++)
1838 {
1839 this->H[i] = (this->Bounds[2 * i + 1] - this->Bounds[2 * i]) / static_cast<double>(ndivs[i]);
1840 }
1841
1842 // Instantiate the locator. The type is related to the maximum point id.
1843 // This is done for performance (e.g., the sort is faster) and significant
1844 // memory savings.
1845 //
1846 if (numPts >= VTK_INT_MAX || numBuckets >= VTK_INT_MAX)
1847 {
1848 this->LargeIds = true;
1849 this->Buckets = new BucketList<vtkIdType>(this, numPts, numBuckets);
1850 }
1851 else
1852 {
1853 this->LargeIds = false;
1854 this->Buckets = new BucketList<int>(this, numPts, numBuckets);
1855 }
1856
1857 // Actually construct the locator
1858 this->Buckets->BuildLocator();
1859
1860 this->BuildTime.Modified();
1861 }
1862
1863 //------------------------------------------------------------------------------
1864 // These methods satisfy the vtkStaticPointLocator API. The implementation is
1865 // with the templated BucketList class. Note that a lot of the complexity here
1866 // is due to the desire to use different id types (int versus vtkIdType) for the
1867 // purposes of increasing speed and reducing memory.
1868 //
1869 // You're probably wondering why an if check (on LargeIds) is used to
1870 // static_cast on BukcetList<T> type, when virtual inheritance could be
1871 // used. Benchmarking shows a small speed difference due to inlining, which
1872 // the use of virtual methods short circuits.
1873
1874 //------------------------------------------------------------------------------
1875 // Given a position x, return the id of the point closest to it.
FindClosestPoint(const double x[3])1876 vtkIdType vtkStaticPointLocator::FindClosestPoint(const double x[3])
1877 {
1878 this->BuildLocator(); // will subdivide if modified; otherwise returns
1879 if (!this->Buckets)
1880 {
1881 return -1;
1882 }
1883
1884 if (this->LargeIds)
1885 {
1886 return static_cast<BucketList<vtkIdType>*>(this->Buckets)->FindClosestPoint(x);
1887 }
1888 else
1889 {
1890 return static_cast<BucketList<int>*>(this->Buckets)->FindClosestPoint(x);
1891 }
1892 }
1893
1894 //------------------------------------------------------------------------------
FindClosestPointWithinRadius(double radius,const double x[3],double inputDataLength,double & dist2)1895 vtkIdType vtkStaticPointLocator::FindClosestPointWithinRadius(
1896 double radius, const double x[3], double inputDataLength, double& dist2)
1897 {
1898 this->BuildLocator(); // will subdivide if modified; otherwise returns
1899 if (!this->Buckets)
1900 {
1901 return -1;
1902 }
1903
1904 if (this->LargeIds)
1905 {
1906 return static_cast<BucketList<vtkIdType>*>(this->Buckets)
1907 ->FindClosestPointWithinRadius(radius, x, inputDataLength, dist2);
1908 }
1909 else
1910 {
1911 return static_cast<BucketList<int>*>(this->Buckets)
1912 ->FindClosestPointWithinRadius(radius, x, inputDataLength, dist2);
1913 }
1914 }
1915
1916 //------------------------------------------------------------------------------
FindClosestPointWithinRadius(double radius,const double x[3],double & dist2)1917 vtkIdType vtkStaticPointLocator::FindClosestPointWithinRadius(
1918 double radius, const double x[3], double& dist2)
1919 {
1920 return this->FindClosestPointWithinRadius(radius, x, this->DataSet->GetLength(), dist2);
1921 }
1922
1923 //------------------------------------------------------------------------------
FindClosestNPoints(int N,const double x[3],vtkIdList * result)1924 void vtkStaticPointLocator::FindClosestNPoints(int N, const double x[3], vtkIdList* result)
1925 {
1926 this->BuildLocator(); // will subdivide if modified; otherwise returns
1927 if (!this->Buckets)
1928 {
1929 return;
1930 }
1931
1932 if (this->LargeIds)
1933 {
1934 return static_cast<BucketList<vtkIdType>*>(this->Buckets)->FindClosestNPoints(N, x, result);
1935 }
1936 else
1937 {
1938 return static_cast<BucketList<int>*>(this->Buckets)->FindClosestNPoints(N, x, result);
1939 }
1940 }
1941
1942 //------------------------------------------------------------------------------
FindPointsWithinRadius(double R,const double x[3],vtkIdList * result)1943 void vtkStaticPointLocator::FindPointsWithinRadius(double R, const double x[3], vtkIdList* result)
1944 {
1945 this->BuildLocator(); // will subdivide if modified; otherwise returns
1946 if (!this->Buckets)
1947 {
1948 return;
1949 }
1950
1951 if (this->LargeIds)
1952 {
1953 return static_cast<BucketList<vtkIdType>*>(this->Buckets)->FindPointsWithinRadius(R, x, result);
1954 }
1955 else
1956 {
1957 return static_cast<BucketList<int>*>(this->Buckets)->FindPointsWithinRadius(R, x, result);
1958 }
1959 }
1960
1961 //------------------------------------------------------------------------------
1962 // This method traverses the locator along the defined ray, finding the
1963 // closest point to a0 when projected onto the line (a0,a1) (i.e., min
1964 // parametric coordinate t) and within the tolerance tol (measured in the
1965 // world coordinate system).
IntersectWithLine(double a0[3],double a1[3],double tol,double & t,double lineX[3],double ptX[3],vtkIdType & ptId)1966 int vtkStaticPointLocator::IntersectWithLine(double a0[3], double a1[3], double tol, double& t,
1967 double lineX[3], double ptX[3], vtkIdType& ptId)
1968 {
1969 this->BuildLocator(); // will subdivide if modified; otherwise returns
1970 if (!this->Buckets)
1971 {
1972 return 0;
1973 }
1974
1975 if (this->LargeIds)
1976 {
1977 return static_cast<BucketList<vtkIdType>*>(this->Buckets)
1978 ->IntersectWithLine(a0, a1, tol, t, lineX, ptX, ptId);
1979 }
1980 else
1981 {
1982 return static_cast<BucketList<int>*>(this->Buckets)
1983 ->IntersectWithLine(a0, a1, tol, t, lineX, ptX, ptId);
1984 }
1985 }
1986
1987 //------------------------------------------------------------------------------
1988 // Build a representation for the locator.
GenerateRepresentation(int level,vtkPolyData * pd)1989 void vtkStaticPointLocator::GenerateRepresentation(int level, vtkPolyData* pd)
1990 {
1991 this->BuildLocator(); // will subdivide if modified; otherwise returns
1992 if (!this->Buckets)
1993 {
1994 return;
1995 }
1996
1997 if (this->LargeIds)
1998 {
1999 return static_cast<BucketList<vtkIdType>*>(this->Buckets)->GenerateRepresentation(level, pd);
2000 }
2001 else
2002 {
2003 return static_cast<BucketList<int>*>(this->Buckets)->GenerateRepresentation(level, pd);
2004 }
2005 }
2006
2007 //------------------------------------------------------------------------------
2008 // Given a bucket, return the number of points inside of it.
GetNumberOfPointsInBucket(vtkIdType bNum)2009 vtkIdType vtkStaticPointLocator::GetNumberOfPointsInBucket(vtkIdType bNum)
2010 {
2011 this->BuildLocator(); // will subdivide if modified; otherwise returns
2012 if (!this->Buckets)
2013 {
2014 return 0;
2015 }
2016
2017 if (this->LargeIds)
2018 {
2019 return static_cast<BucketList<vtkIdType>*>(this->Buckets)->GetNumberOfIds(bNum);
2020 }
2021 else
2022 {
2023 return static_cast<BucketList<int>*>(this->Buckets)->GetNumberOfIds(bNum);
2024 }
2025 }
2026
2027 //------------------------------------------------------------------------------
2028 // Given a bucket, return the ids in the bucket.
GetBucketIds(vtkIdType bNum,vtkIdList * bList)2029 void vtkStaticPointLocator::GetBucketIds(vtkIdType bNum, vtkIdList* bList)
2030 {
2031 this->BuildLocator(); // will subdivide if modified; otherwise returns
2032 if (!this->Buckets)
2033 {
2034 bList->Reset();
2035 return;
2036 }
2037
2038 if (this->LargeIds)
2039 {
2040 return static_cast<BucketList<vtkIdType>*>(this->Buckets)->GetIds(bNum, bList);
2041 }
2042 else
2043 {
2044 return static_cast<BucketList<int>*>(this->Buckets)->GetIds(bNum, bList);
2045 }
2046 }
2047
2048 //------------------------------------------------------------------------------
2049 // Given a bucket, return the ids in the bucket.
MergePoints(double tol,vtkIdType * pointMap)2050 void vtkStaticPointLocator::MergePoints(double tol, vtkIdType* pointMap)
2051 {
2052 this->BuildLocator(); // will subdivide if modified; otherwise returns
2053 if (!this->Buckets)
2054 {
2055 return;
2056 }
2057
2058 if (this->LargeIds)
2059 {
2060 return static_cast<BucketList<vtkIdType>*>(this->Buckets)->MergePoints(tol, pointMap);
2061 }
2062 else
2063 {
2064 return static_cast<BucketList<int>*>(this->Buckets)->MergePoints(tol, pointMap);
2065 }
2066 }
2067
2068 //------------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)2069 void vtkStaticPointLocator::PrintSelf(ostream& os, vtkIndent indent)
2070 {
2071 this->Superclass::PrintSelf(os, indent);
2072
2073 os << indent << "Number of Points Per Bucket: " << this->NumberOfPointsPerBucket << "\n";
2074
2075 os << indent << "Divisions: (" << this->Divisions[0] << ", " << this->Divisions[1] << ", "
2076 << this->Divisions[2] << ")\n";
2077
2078 os << indent << "Max Number Of Buckets: " << this->MaxNumberOfBuckets << "\n";
2079
2080 os << indent << "Large IDs: " << this->LargeIds << "\n";
2081 }
2082