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