1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkStaticCellLocator.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 "vtkStaticCellLocator.h"
16
17 #include "vtkBoundingBox.h"
18 #include "vtkBox.h"
19 #include "vtkCellArray.h"
20 #include "vtkDoubleArray.h"
21 #include "vtkGenericCell.h"
22 #include "vtkIdList.h"
23 #include "vtkIntArray.h"
24 #include "vtkMath.h"
25 #include "vtkMergePoints.h"
26 #include "vtkObjectFactory.h"
27 #include "vtkPlane.h"
28 #include "vtkPoints.h"
29 #include "vtkPolyData.h"
30 #include "vtkSMPThreadLocal.h"
31 #include "vtkSMPThreadLocalObject.h"
32 #include "vtkSMPTools.h"
33
34 #include <functional>
35 #include <queue>
36 #include <vector>
37
38 vtkStandardNewMacro(vtkStaticCellLocator);
39
40 //------------------------------------------------------------------------------
41 // Helper classes to support efficient computing, and threaded execution.
42 //
43 // Note that are two key classes: the vtkCellBinner and the vtkCellProcessor.
44 // the Binner is used to perform binning operations are cells are placed into
45 // the uniformally subdivided bin space. The Processor is a templated class
46 // (templated over ID type to reduce memory and speed sorting when the cell
47 // ids are small).
48 //
49 // The algorithm is multipass. First, the overall bounds of the data is
50 // determined, and then the space is subdivided into uniform bins. Next each
51 // cell is visited, it's bounds are obtained. and the ijk footprint into the
52 // binning is obtained. The footprint implicitly indicates the number of bins
53 // the cell touches (i.e., the number of cell fragment tuples
54 // (cellId,binId)), and this number is stored in an array. Once all cells
55 // have been visited (in parallel), a prefix sum is executed on the counts to
56 // determine the total number of fragments. Next another (parallel) pass is
57 // made over each cell and the fragments are placed into a tuple array (using
58 // the count offsets), which is then (parallel) sorted on binIds. This
59 // produces contiguous runs of cellIds for each bin. Finally, integral
60 // offsets are created that for each cell point at the beginning of each
61 // run.
62 //
63 // This algorithm is implemented into two parts as mentioned previously. The
64 // Binner is non templated and simply deals with cell bounds and eventually
65 // computes the number of cell fragments. Depending on the size of the
66 // fragment count, a templated class of either int or vtkIdType is
67 // created. Different types are used because of 1) significant reduction in
68 // memory and 2) significant speed up in the parallel sort.
69
70 //========================= CELL LOCATOR MACHINERY ====================================
71
72 // PIMPLd class which wraps binning functionality.
73 struct vtkCellBinner
74 {
75 vtkStaticCellLocator* Locator; // locater
76 vtkIdType NumCells; // the number of cells to bin
77 vtkIdType NumBins;
78 vtkIdType NumFragments; // total number of (cellId,binId) tuples
79
80 // These are internal data members used for performance reasons
81 vtkDataSet* DataSet;
82 int Divisions[3];
83 double Bounds[6];
84 double* CellBounds;
85 vtkIdType* Counts;
86 double H[3];
87 double hX, hY, hZ;
88 double fX, fY, fZ, bX, bY, bZ;
89 vtkIdType xD, yD, zD, xyD;
90 double binTol;
91
92 // Construction
vtkCellBinnervtkCellBinner93 vtkCellBinner(vtkStaticCellLocator* loc, vtkIdType numCells, vtkIdType numBins)
94 {
95 this->Locator = loc;
96 this->NumCells = numCells;
97 this->NumBins = numBins;
98 this->NumFragments = 0;
99 this->DataSet = loc->GetDataSet();
100 loc->GetDivisions(this->Divisions);
101
102 // Allocate data. Note that these arrays are deleted elsewhere
103 this->CellBounds = new double[numCells * 6];
104 this->Counts = new vtkIdType[numCells + 1]; // one extra holds total count
105
106 // This is done to cause non-thread safe initialization to occur due to
107 // side effects from GetCellBounds().
108 this->DataSet->GetCellBounds(0, this->CellBounds);
109
110 // Setup internal data members for more efficient processing.
111 this->hX = this->H[0] = loc->H[0];
112 this->hY = this->H[1] = loc->H[1];
113 this->hZ = this->H[2] = loc->H[2];
114 this->fX = 1.0 / loc->H[0];
115 this->fY = 1.0 / loc->H[1];
116 this->fZ = 1.0 / loc->H[2];
117 this->bX = this->Bounds[0] = loc->Bounds[0];
118 this->Bounds[1] = loc->Bounds[1];
119 this->bY = this->Bounds[2] = loc->Bounds[2];
120 this->Bounds[3] = loc->Bounds[3];
121 this->bZ = this->Bounds[4] = loc->Bounds[4];
122 this->Bounds[5] = loc->Bounds[5];
123 this->xD = this->Divisions[0];
124 this->yD = this->Divisions[1];
125 this->zD = this->Divisions[2];
126 this->xyD = this->Divisions[0] * this->Divisions[1];
127
128 this->binTol = loc->GetUseDiagonalLengthTolerance()
129 ? loc->GetTolerance() * sqrt(this->hX * this->hX + this->hY * this->hY + this->hZ * this->hZ)
130 : loc->GetTolerance();
131 }
132
~vtkCellBinnervtkCellBinner133 ~vtkCellBinner()
134 {
135 delete[] this->CellBounds;
136 delete[] this->Counts;
137 }
138
GetBinIndicesvtkCellBinner139 void GetBinIndices(const double* x, int ijk[3]) const
140 {
141 // Compute point index. Make sure it lies within range of locator.
142 ijk[0] = static_cast<int>(((x[0] - bX) * fX));
143 ijk[1] = static_cast<int>(((x[1] - bY) * fY));
144 ijk[2] = static_cast<int>(((x[2] - bZ) * fZ));
145
146 ijk[0] = (ijk[0] < 0 ? 0 : (ijk[0] >= xD ? xD - 1 : ijk[0]));
147 ijk[1] = (ijk[1] < 0 ? 0 : (ijk[1] >= yD ? yD - 1 : ijk[1]));
148 ijk[2] = (ijk[2] < 0 ? 0 : (ijk[2] >= zD ? zD - 1 : ijk[2]));
149 }
150
GetBinIndicesvtkCellBinner151 void GetBinIndices(vtkIdType binId, int ijk[3]) const
152 {
153 ijk[0] = binId % xD;
154 vtkIdType tmp = binId / xD;
155 ijk[1] = tmp % yD;
156 ijk[2] = tmp / yD;
157 }
158
159 // Given a point x, determine which bin it is in. Note that points
160 // are clamped to lie inside of the locator.
GetBinIndexvtkCellBinner161 vtkIdType GetBinIndex(const double* x) const
162 {
163 int ijk[3];
164 this->GetBinIndices(x, ijk);
165 return ijk[0] + ijk[1] * xD + ijk[2] * xyD;
166 }
167
GetBinIndexvtkCellBinner168 vtkIdType GetBinIndex(int ijk[3]) const { return ijk[0] + ijk[1] * xD + ijk[2] * xyD; }
169
170 // These are helper functions
CountBinsvtkCellBinner171 vtkIdType CountBins(const int ijkMin[3], const int ijkMax[3])
172 {
173 // Ensure all temporary values are vtkIdType:
174 vtkIdType result = ijkMax[0] - ijkMin[0] + 1;
175 result *= ijkMax[1] - ijkMin[1] + 1;
176 result *= ijkMax[2] - ijkMin[2] + 1;
177
178 return result;
179 }
180
InitializevtkCellBinner181 void Initialize() {}
182
operator ()vtkCellBinner183 void operator()(vtkIdType cellId, vtkIdType endCellId)
184 {
185 double* bds = this->CellBounds + cellId * 6;
186 vtkIdType* counts = this->Counts + cellId;
187 double xmin[3], xmax[3];
188 int ijkMin[3], ijkMax[3];
189
190 for (; cellId < endCellId; ++cellId, bds += 6)
191 {
192 this->DataSet->GetCellBounds(cellId, bds);
193 double tol = this->binTol;
194 xmin[0] = bds[0] - tol;
195 xmin[1] = bds[2] - tol;
196 xmin[2] = bds[4] - tol;
197 xmax[0] = bds[1] + tol;
198 xmax[1] = bds[3] + tol;
199 xmax[2] = bds[5] + tol;
200
201 this->GetBinIndices(xmin, ijkMin);
202 this->GetBinIndices(xmax, ijkMax);
203
204 *counts++ = this->CountBins(ijkMin, ijkMax);
205 }
206 }
207
ReducevtkCellBinner208 void Reduce()
209 {
210 // Perform prefix sum
211 vtkIdType* counts = this->Counts;
212 vtkIdType numBins, total = 0, numCells = this->NumCells;
213 for (vtkIdType i = 0; i < numCells; ++i)
214 {
215 numBins = *counts;
216 *counts++ = total;
217 total += numBins;
218 }
219 this->NumFragments = total;
220 }
221
222 }; // vtkCellBinner
223
224 //------------------------------------------------------------------------------
225 // The following tuple is what is sorted in the map. Note that it is templated
226 // because depending on the number of points / buckets to process we may want
227 // to use vtkIdType. Otherwise for performance reasons it's best to use an int
228 // (or other integral type). Typically sort() is 25-30% faster on smaller
229 // integral types, plus it takes a heck less memory (when vtkIdType is 64-bit
230 // and int is 32-bit).
231 template <typename TId>
232 struct CellFragments
233 {
234 TId CellId; // originating cell id
235 TId BinId; // i-j-k index into bin space
236
237 // Operator< used to support the subsequent sort operation.
operator <CellFragments238 bool operator<(const CellFragments& tuple) const { return BinId < tuple.BinId; }
239 };
240
241 // Perform locator operations like FindCell. Uses templated subclasses
242 // to reduce memory and enhance speed.
243 struct vtkCellProcessor
244 {
245 vtkCellBinner* Binner;
246 vtkDataSet* DataSet;
247 double* CellBounds;
248 vtkIdType* Counts;
249 vtkIdType NumFragments;
250 vtkIdType NumCells;
251 vtkIdType NumBins;
252 int BatchSize;
253 int NumBatches;
254 vtkIdType xD, xyD;
255
vtkCellProcessorvtkCellProcessor256 vtkCellProcessor(vtkCellBinner* cb)
257 : Binner(cb)
258 {
259 this->DataSet = cb->DataSet;
260 this->CellBounds = cb->CellBounds;
261 this->Counts = cb->Counts;
262 this->NumCells = cb->NumCells;
263 this->NumFragments = cb->NumFragments;
264 this->NumBins = cb->NumBins;
265 this->BatchSize = 10000; // building the offset array
266 this->NumBatches =
267 static_cast<int>(ceil(static_cast<double>(this->NumFragments) / this->BatchSize));
268
269 xD = cb->xD; // for speeding up computation
270 xyD = cb->xyD;
271 }
272
273 virtual ~vtkCellProcessor() = default;
274
275 // Satisfy cell locator API
276 virtual vtkIdType FindCell(
277 const double pos[3], vtkGenericCell* cell, double pcoords[3], double* weights) = 0;
278 virtual void FindCellsWithinBounds(double* bbox, vtkIdList* cells) = 0;
279 virtual void FindCellsAlongLine(
280 const double p1[3], const double p2[3], double tol, vtkIdList* cells) = 0;
281 virtual void FindCellsAlongPlane(
282 const double o[3], const double n[3], double tolerance, vtkIdList* cells) = 0;
283 virtual int IntersectWithLine(const double a0[3], const double a1[3], double tol, double& t,
284 double x[3], double pcoords[3], int& subId, vtkIdType& cellId, vtkGenericCell* cell) = 0;
285 virtual vtkIdType FindClosestPointWithinRadius(const double x[3], double radius,
286 double closestPoint[3], vtkGenericCell* cell, vtkIdType& cellId, int& subId, double& dist2,
287 int& inside) = 0;
288
289 // Convenience for computing
290 virtual int IsEmpty(vtkIdType binId) = 0;
291 };
292
293 namespace
294 { // anonymous to wrap non-public stuff
295
296 // Typed subclass
297 template <typename T>
298 struct CellProcessor : public vtkCellProcessor
299 {
300 // Type dependent members
301 CellFragments<T>* Map; // the map to be sorted
302 T* Offsets; // offsets for each bin into the map
303
CellProcessor__anond231c8eb0111::CellProcessor304 CellProcessor(vtkCellBinner* cb)
305 : vtkCellProcessor(cb)
306 {
307 // Prepare to sort
308 // one extra to simplify traversal
309 this->Map = new CellFragments<T>[this->NumFragments + 1];
310 this->Map[this->NumFragments].BinId = this->NumBins;
311 this->Offsets = new T[this->NumBins + 1];
312 this->Offsets[this->NumBins] = this->NumFragments;
313 }
314
~CellProcessor__anond231c8eb0111::CellProcessor315 ~CellProcessor() override
316 {
317 delete[] this->Map;
318 delete[] this->Offsets;
319 }
320
321 // The number of cell ids in a bin is determined by computing the
322 // difference between the offsets into the sorted cell fragments array.
GetNumberOfIds__anond231c8eb0111::CellProcessor323 T GetNumberOfIds(vtkIdType binNum) { return (this->Offsets[binNum + 1] - this->Offsets[binNum]); }
324
325 // Given a bin number, return the cells ids in that bin.
GetIds__anond231c8eb0111::CellProcessor326 const CellFragments<T>* GetIds(vtkIdType binNum) { return this->Map + this->Offsets[binNum]; }
327
ComputeBinBounds__anond231c8eb0111::CellProcessor328 void ComputeBinBounds(int i, int j, int k, double binBounds[6])
329 {
330 double* bds = this->Binner->Bounds;
331 double* h = this->Binner->H;
332 binBounds[0] = bds[0] + i * h[0];
333 binBounds[1] = binBounds[0] + h[0];
334 binBounds[2] = bds[2] + j * h[1];
335 binBounds[3] = binBounds[2] + h[1];
336 binBounds[4] = bds[4] + k * h[2];
337 binBounds[5] = binBounds[4] + h[2];
338 }
339
IsInBinBounds__anond231c8eb0111::CellProcessor340 int IsInBinBounds(double binBounds[6], double x[3], double binTol = 0.0)
341 {
342 if ((binBounds[0] - binTol) <= x[0] && x[0] <= (binBounds[1] + binTol) &&
343 (binBounds[2] - binTol) <= x[1] && x[1] <= (binBounds[3] + binTol) &&
344 (binBounds[4] - binTol) <= x[2] && x[2] <= (binBounds[5] + binTol))
345 {
346 return 1;
347 }
348 else
349 {
350 return 0;
351 }
352 }
353
354 // Methods to satisfy vtkCellProcessor virtual API
355 vtkIdType FindCell(
356 const double pos[3], vtkGenericCell* cell, double pcoords[3], double* weights) override;
357 void FindCellsWithinBounds(double* bbox, vtkIdList* cells) override;
358 void FindCellsAlongLine(
359 const double p1[3], const double p2[3], double tol, vtkIdList* cells) override;
360 void FindCellsAlongPlane(
361 const double o[3], const double n[3], double tolerance, vtkIdList* cells) override;
362 int IntersectWithLine(const double a0[3], const double a1[3], double tol, double& t, double x[3],
363 double pcoords[3], int& subId, vtkIdType& cellId, vtkGenericCell* cell) override;
364 vtkIdType FindClosestPointWithinRadius(const double x[3], double radius, double closestPoint[3],
365 vtkGenericCell* cell, vtkIdType& cellId, int& subId, double& dist2, int& inside) override;
IsEmpty__anond231c8eb0111::CellProcessor366 int IsEmpty(vtkIdType binId) override
367 {
368 return (this->GetNumberOfIds(static_cast<T>(binId)) > 0 ? 0 : 1);
369 }
370
371 // This functor is used to perform the final cell binning
Initialize__anond231c8eb0111::CellProcessor372 void Initialize() {}
373
operator ()__anond231c8eb0111::CellProcessor374 void operator()(vtkIdType cellId, vtkIdType endCellId)
375 {
376 const double* bds = this->CellBounds + cellId * 6;
377 CellFragments<T>* t = this->Map + *(this->Counts + cellId);
378 double xmin[3], xmax[3];
379 int ijkMin[3], ijkMax[3];
380 int i, j, k;
381 vtkIdType binId;
382
383 for (; cellId < endCellId; ++cellId, bds += 6)
384 {
385 double tol = this->Binner->binTol;
386 xmin[0] = bds[0] - tol;
387 xmin[1] = bds[2] - tol;
388 xmin[2] = bds[4] - tol;
389 xmax[0] = bds[1] + tol;
390 xmax[1] = bds[3] + tol;
391 xmax[2] = bds[5] + tol;
392
393 this->Binner->GetBinIndices(xmin, ijkMin);
394 this->Binner->GetBinIndices(xmax, ijkMax);
395
396 for (k = ijkMin[2]; k <= ijkMax[2]; ++k)
397 {
398 for (j = ijkMin[1]; j <= ijkMax[1]; ++j)
399 {
400 for (i = ijkMin[0]; i <= ijkMax[0]; ++i)
401 {
402 binId = i + j * xD + k * xyD;
403 t->CellId = cellId;
404 t->BinId = binId;
405 t++;
406 }
407 }
408 }
409 }
410 }
411
Reduce__anond231c8eb0111::CellProcessor412 void Reduce() {}
413 }; // CellProcessor
414
415 // This functor class creates offsets for each cell into the sorted tuple
416 // array. The offsets enable random access to cells.
417 template <typename TId>
418 struct MapOffsets
419 {
420 CellProcessor<TId>* Processor;
421 CellFragments<TId>* Map;
422 TId* Offsets;
423 vtkIdType NumCells;
424 int NumBins;
425 vtkIdType NumFragments;
426 int BatchSize;
427
MapOffsets__anond231c8eb0111::MapOffsets428 MapOffsets(CellProcessor<TId>* p)
429 : Processor(p)
430 {
431 this->Map = p->Map;
432 this->Offsets = p->Offsets;
433 this->NumCells = p->NumCells;
434 this->NumBins = p->NumBins;
435 this->NumFragments = p->NumFragments;
436 this->BatchSize = p->BatchSize;
437 }
438
439 // Traverse sorted points (i.e., tuples) and update bin offsets.
operator ()__anond231c8eb0111::MapOffsets440 void operator()(vtkIdType batch, vtkIdType batchEnd)
441 {
442 TId* offsets = this->Offsets;
443 const CellFragments<TId>* curPt = this->Map + batch * this->BatchSize;
444 const CellFragments<TId>* endBatchPt = this->Map + batchEnd * this->BatchSize;
445 const CellFragments<TId>* endPt = this->Map + this->NumFragments;
446 const CellFragments<TId>* prevPt;
447 endBatchPt = (endBatchPt > endPt ? endPt : endBatchPt);
448
449 // Special case at the very beginning of the mapped points array. If
450 // the first point is in bin# N, then all bins up and including
451 // N must refer to the first point.
452 if (curPt == this->Map)
453 {
454 prevPt = this->Map;
455 std::fill_n(offsets, curPt->BinId + 1, 0); // point to the first points
456 } // at the very beginning of the map (sorted points array)
457
458 // We are entering this functor somewhere in the interior of the
459 // mapped points array. All we need to do is point to the entry
460 // position because we are interested only in prevPt->BinId.
461 else
462 {
463 prevPt = curPt;
464 } // else in the middle of a batch
465
466 // Okay we have a starting point for a bin run. Now we can begin
467 // filling in the offsets in this batch. A previous thread should
468 // have/will have completed the previous and subsequent runs outside
469 // of the [batch,batchEnd) range
470 for (curPt = prevPt; curPt < endBatchPt;)
471 {
472 for (; curPt->BinId == prevPt->BinId && curPt <= endBatchPt; ++curPt)
473 {
474 ; // advance
475 }
476 // Fill in any gaps in the offset array
477 std::fill_n(offsets + prevPt->BinId + 1, curPt->BinId - prevPt->BinId, curPt - this->Map);
478 prevPt = curPt;
479 } // for all batches in this range
480 } // operator()
481
482 }; // MapOffsets
483
484 //------------------------------------------------------------------------------
485 template <typename T>
FindCell(const double pos[3],vtkGenericCell * cell,double pcoords[3],double * weights)486 vtkIdType CellProcessor<T>::FindCell(
487 const double pos[3], vtkGenericCell* cell, double pcoords[3], double* weights)
488 {
489 vtkIdType binId = this->Binner->GetBinIndex(pos);
490 T numIds = this->GetNumberOfIds(binId);
491
492 // Only thread the evaluation if enough cells need to be processed
493 if (numIds < 1)
494 {
495 return -1;
496 }
497 // Run through serially. A parallel implementation is possible but does
498 // not seem to be much faster.
499 else
500 {
501 const CellFragments<T>* cellIds = this->GetIds(binId);
502 double tol = this->Binner->binTol;
503 double dist2, *bounds, delta[3] = { tol, tol, tol };
504 int subId;
505 vtkIdType cellId;
506
507 for (int j = 0; j < numIds; j++)
508 {
509 cellId = cellIds[j].CellId;
510 bounds = this->CellBounds + 6 * cellId;
511
512 if (vtkMath::PointIsWithinBounds(pos, bounds, delta))
513 {
514 this->DataSet->GetCell(cellId, cell);
515 if (cell->EvaluatePosition(pos, nullptr, subId, pcoords, dist2, weights) == 1)
516 {
517 return cellId;
518 }
519 } // in bounding box
520 } // for cells in this bin
521
522 return -1; // nothing found
523 } // serial
524 }
525
526 //------------------------------------------------------------------------------
527 template <typename T>
FindCellsWithinBounds(double * bbox,vtkIdList * cells)528 void CellProcessor<T>::FindCellsWithinBounds(double* bbox, vtkIdList* cells)
529 {
530 vtkIdType binNum, numIds, jOffset, kOffset;
531 int i, j, k, ii, ijkMin[3], ijkMax[3];
532 double pMin[3], pMax[3];
533 const CellFragments<T>* ids;
534
535 cells->Reset();
536
537 // Get the locator locations for the two extreme corners of the bounding box
538 pMin[0] = bbox[0];
539 pMin[1] = bbox[2];
540 pMin[2] = bbox[4];
541 pMax[0] = bbox[1];
542 pMax[1] = bbox[3];
543 pMax[2] = bbox[5];
544
545 this->Binner->GetBinIndices(pMin, ijkMin);
546 this->Binner->GetBinIndices(pMax, ijkMax);
547
548 // Loop over the block of bins and add cells that have not yet been visited.
549 for (k = ijkMin[2]; k <= ijkMax[2]; ++k)
550 {
551 kOffset = k * this->xyD;
552 for (j = ijkMin[1]; j <= ijkMax[1]; ++j)
553 {
554 jOffset = j * this->xD;
555 for (i = ijkMin[0]; i <= ijkMax[0]; ++i)
556 {
557 binNum = i + jOffset + kOffset;
558
559 if ((numIds = this->GetNumberOfIds(binNum)) > 0)
560 {
561 ids = this->GetIds(binNum);
562 for (ii = 0; ii < numIds; ii++)
563 {
564 // Could use query mechanism to speed up at some point
565 cells->InsertUniqueId(ids[ii].CellId);
566 } // for all points in bucket
567 } // if points in bucket
568 } // i-footprint
569 } // j-footprint
570 } // k-footprint
571 }
572
573 //------------------------------------------------------------------------------
574 // This code traverses the cell locator by following the intersection ray. All
575 // cells in intersected bins are placed into the output cellId vtkIdList. See
576 // the IntersectWithLine method for more information on voxel traversal.
577 template <typename T>
FindCellsAlongLine(const double a0[3],const double a1[3],double vtkNotUsed (tol),vtkIdList * cells)578 void CellProcessor<T>::FindCellsAlongLine(
579 const double a0[3], const double a1[3], double vtkNotUsed(tol), vtkIdList* cells)
580 {
581 // Initialize the list of cells
582 cells->Reset();
583
584 // Make sure the bounding box of the locator is hit. Also, determine the
585 // entry and exit points into and out of the locator. This is used to
586 // determine the bins where the ray starts and ends.
587 double* bounds = this->Binner->Bounds;
588 double t0, t1, x0[3], x1[3], rayDir[3];
589 int plane0, plane1;
590 vtkMath::Subtract(a1, a0, rayDir);
591
592 if (vtkBox::IntersectWithLine(bounds, a0, a1, t0, t1, x0, x1, plane0, plane1) == 0)
593 {
594 return; // No intersections possible, line is outside the locator
595 }
596
597 // Okay process line
598 int* ndivs = this->Binner->Divisions;
599 vtkIdType prod = ndivs[0] * ndivs[1];
600 double* h = this->Binner->H;
601 T i, numCellsInBin;
602 int ijk[3], ijkEnd[3];
603 vtkIdType idx, cId;
604 double hitCellBoundsPosition[3], tHitCell;
605 double step[3], next[3], tMax[3], tDelta[3];
606 double binBounds[6];
607
608 // Initialize intersection query array if necessary. This is done
609 // locally to ensure thread safety.
610 std::vector<unsigned char> cellHasBeenVisited(this->NumCells, 0);
611
612 // Get the i-j-k point of intersection and bin index. This is
613 // clamped to the boundary of the locator. Also get the exit bin
614 // of the line.
615 this->Binner->GetBinIndices(x0, ijk);
616 this->Binner->GetBinIndices(x1, ijkEnd);
617 idx = ijk[0] + ijk[1] * ndivs[0] + ijk[2] * prod;
618
619 // Set up some traversal parameters for traversing through bins
620 step[0] = (rayDir[0] >= 0.0) ? 1.0 : -1.0;
621 step[1] = (rayDir[1] >= 0.0) ? 1.0 : -1.0;
622 step[2] = (rayDir[2] >= 0.0) ? 1.0 : -1.0;
623
624 // If the ray is going in the negative direction, then the next voxel boundary
625 // is on the "-" direction so we stay in the current voxel.
626 next[0] = bounds[0] + h[0] * (rayDir[0] >= 0.0 ? (ijk[0] + step[0]) : ijk[0]);
627 next[1] = bounds[2] + h[1] * (rayDir[1] >= 0.0 ? (ijk[1] + step[1]) : ijk[1]);
628 next[2] = bounds[4] + h[2] * (rayDir[2] >= 0.0 ? (ijk[2] + step[2]) : ijk[2]);
629
630 tMax[0] = (rayDir[0] != 0.0) ? (next[0] - x0[0]) / rayDir[0] : VTK_FLOAT_MAX;
631 tMax[1] = (rayDir[1] != 0.0) ? (next[1] - x0[1]) / rayDir[1] : VTK_FLOAT_MAX;
632 tMax[2] = (rayDir[2] != 0.0) ? (next[2] - x0[2]) / rayDir[2] : VTK_FLOAT_MAX;
633
634 tDelta[0] = (rayDir[0] != 0.0) ? (h[0] / rayDir[0]) * step[0] : VTK_FLOAT_MAX;
635 tDelta[1] = (rayDir[1] != 0.0) ? (h[1] / rayDir[1]) * step[1] : VTK_FLOAT_MAX;
636 tDelta[2] = (rayDir[2] != 0.0) ? (h[2] / rayDir[2]) * step[2] : VTK_FLOAT_MAX;
637
638 // Start walking through the bins, continue until traversed the entire
639 // locator. Note that termination of the while(1) loop occurs when the ray
640 // passes out of the locator via the break command.
641 while (true)
642 {
643 if ((numCellsInBin = this->GetNumberOfIds(idx)) > 0) // there are some cell here
644 {
645 const CellFragments<T>* cellIds = this->GetIds(idx);
646 this->ComputeBinBounds(ijk[0], ijk[1], ijk[2], binBounds);
647 for (i = 0; i < numCellsInBin; i++)
648 {
649 cId = cellIds[i].CellId;
650 if (cellHasBeenVisited[cId] == 0)
651 {
652 cellHasBeenVisited[cId] = 1;
653
654 // check whether we intersect the cell bounds
655 int hitCellBounds = vtkBox::IntersectBox(this->CellBounds + (6 * cId), a0, rayDir,
656 hitCellBoundsPosition, tHitCell, this->Binner->binTol);
657
658 if (hitCellBounds)
659 {
660 // Note because of cellHasBeenVisited[], we know this cId is unique
661 cells->InsertNextId(cId);
662 } // if (hitCellBounds)
663 } // if (!cellHasBeenVisited[cId])
664 } // over all cells in bin
665 } // if cells in bin
666
667 // See if the traversal is complete (reached the end of the line).
668 if (ijk[0] == ijkEnd[0] && ijk[1] == ijkEnd[1] && ijk[2] == ijkEnd[2])
669 {
670 break;
671 }
672
673 // Advance to the next voxel
674 if (tMax[0] < tMax[1])
675 {
676 if (tMax[0] < tMax[2])
677 {
678 ijk[0] += static_cast<int>(step[0]);
679 tMax[0] += tDelta[0];
680 }
681 else
682 {
683 ijk[2] += static_cast<int>(step[2]);
684 tMax[2] += tDelta[2];
685 }
686 }
687 else
688 {
689 if (tMax[1] < tMax[2])
690 {
691 ijk[1] += static_cast<int>(step[1]);
692 tMax[1] += tDelta[1];
693 }
694 else
695 {
696 ijk[2] += static_cast<int>(step[2]);
697 tMax[2] += tDelta[2];
698 }
699 }
700
701 // If exiting the locator, we are done
702 if (ijk[0] < 0 || ijk[0] >= ndivs[0] || ijk[1] < 0 || ijk[1] >= ndivs[1] || ijk[2] < 0 ||
703 ijk[2] >= ndivs[2])
704 {
705 break;
706 }
707 else
708 {
709 idx = ijk[0] + ijk[1] * ndivs[0] + ijk[2] * prod;
710 }
711 } // while still in locator
712 }
713
714 // This functor identifies candidate cells as to whether they may intersect
715 // a specified plane. Locator bins are culled first, and if they intersect
716 // the plane, then the cell bounding boxes are used.
717 template <typename TId>
718 struct CellPlaneCandidates
719 {
720 CellProcessor<TId>* Processor;
721 const vtkCellBinner* Binner;
722 double Origin[3];
723 double Normal[3];
724 unsigned char* CellVisited;
725 double BinOffsetX, BinOffsetY, BinOffsetZ;
726 double BinRadius;
727
CellPlaneCandidates__anond231c8eb0111::CellPlaneCandidates728 CellPlaneCandidates(CellProcessor<TId>* p, const vtkCellBinner* b, const double* o,
729 const double* n, unsigned char* visited)
730 : Processor(p)
731 , Binner(b)
732 , CellVisited(visited)
733 {
734 this->Origin[0] = o[0];
735 this->Origin[1] = o[1];
736 this->Origin[2] = o[2];
737
738 this->Normal[0] = n[0];
739 this->Normal[1] = n[1];
740 this->Normal[2] = n[2];
741 vtkMath::Normalize(this->Normal);
742
743 // Offset from the bin origin to the bin center
744 this->BinOffsetX = this->Binner->hX / 2.0;
745 this->BinOffsetY = this->Binner->hY / 2.0;
746 this->BinOffsetZ = this->Binner->hZ / 2.0;
747
748 // The BinRadius is used to cull bins quickly. It's a variant of a sphere
749 // tree test (with the center of a sphere corrsponding to the center of a
750 // bin). Note that the plane orientation affects the radius: the end
751 // result is that a smaller sphere radius can typically be used (as
752 // compared to using the 0.5*(diagonal length) of a bin). This radius
753 // needs only to be computed once since the relative orientation of each
754 // bin to the plane is unchanged during processing. The bin radius is
755 // simply the maximum distance that one of the eight bin corner points is
756 // away from a plane passing through the center of the bin.
757 double x[3], d, dMax = 0.0;
758 x[0] = -this->BinOffsetX;
759 x[1] = -this->BinOffsetY;
760 x[2] = -this->BinOffsetZ;
761 d = vtkMath::Dot(x, this->Normal); // simplified because plane passes through origin
762 dMax = (d > dMax ? d : dMax);
763
764 x[0] = this->BinOffsetX;
765 x[1] = -this->BinOffsetY;
766 x[2] = -this->BinOffsetZ;
767 d = vtkMath::Dot(x, this->Normal);
768 dMax = (d > dMax ? d : dMax);
769
770 x[0] = -this->BinOffsetX;
771 x[1] = this->BinOffsetY;
772 x[2] = -this->BinOffsetZ;
773 d = vtkMath::Dot(x, this->Normal);
774 dMax = (d > dMax ? d : dMax);
775
776 x[0] = this->BinOffsetX;
777 x[1] = this->BinOffsetY;
778 x[2] = -this->BinOffsetZ;
779 d = vtkMath::Dot(x, this->Normal);
780 dMax = (d > dMax ? d : dMax);
781
782 x[0] = -this->BinOffsetX;
783 x[1] = -this->BinOffsetY;
784 x[2] = this->BinOffsetZ;
785 d = vtkMath::Dot(x, this->Normal);
786 dMax = (d > dMax ? d : dMax);
787
788 x[0] = this->BinOffsetX;
789 x[1] = -this->BinOffsetY;
790 x[2] = this->BinOffsetZ;
791 d = vtkMath::Dot(x, this->Normal);
792 dMax = (d > dMax ? d : dMax);
793
794 x[0] = -this->BinOffsetX;
795 x[1] = this->BinOffsetY;
796 x[2] = this->BinOffsetZ;
797 d = vtkMath::Dot(x, this->Normal);
798 dMax = (d > dMax ? d : dMax);
799
800 x[0] = this->BinOffsetX;
801 x[1] = this->BinOffsetY;
802 x[2] = this->BinOffsetZ;
803 d = vtkMath::Dot(x, this->Normal);
804 dMax = (d > dMax ? d : dMax);
805
806 this->BinRadius = dMax;
807 }
808
809 // Operate on slabs (e.g., z-slab) of bins. The algorithm works by checking
810 // whether the current bin is intersected by the plane; if it is, then the
811 // cell bounding box is evaluated as well. Note a potential data race
812 // situation exists since a cell may be marked simultaneously (using the
813 // same value). Since the same value will always be written this has never
814 // been found to be an issue, although it is a potential problem.
operator ()__anond231c8eb0111::CellPlaneCandidates815 void operator()(vtkIdType k, vtkIdType kEnd)
816 {
817 double *o = this->Origin, *n = this->Normal;
818 double center[3], d, *bounds;
819 vtkIdType i, j, iEnd = this->Binner->Divisions[0], jEnd = this->Binner->Divisions[1];
820
821 // Over all z slabs
822 for (; k < kEnd; ++k)
823 {
824 center[2] = this->Binner->Bounds[4] + k * this->Binner->hZ + this->BinOffsetZ;
825 for (j = 0; j < jEnd; ++j)
826 {
827 center[1] = this->Binner->Bounds[2] + j * this->Binner->hY + this->BinOffsetY;
828 for (i = 0; i < iEnd; ++i)
829 {
830 center[0] = this->Binner->Bounds[0] + i * this->Binner->hX + this->BinOffsetX;
831
832 // Now see if the bin could be intersected by the plane. It's a pseudo
833 // sphere tree operation.
834 d = vtkPlane::DistanceToPlane(center, n, o);
835 if (d <= this->BinRadius)
836 {
837 vtkIdType cId, ii, numCellsInBin;
838 vtkIdType bin = i + j * this->Binner->xD + k * this->Binner->xyD;
839 if ((numCellsInBin = this->Processor->GetNumberOfIds(bin)) >
840 0) // there are some cell here
841 {
842 const CellFragments<TId>* cellIds = this->Processor->GetIds(bin);
843 for (ii = 0; ii < numCellsInBin; ii++)
844 {
845 cId = cellIds[ii].CellId;
846 if (!this->CellVisited[cId])
847 {
848 bounds = this->Processor->CellBounds + (6 * cId);
849 this->CellVisited[cId] = (vtkBox::IntersectWithPlane(bounds, o, n) ? 2 : 1);
850 }
851 }
852 }
853 } // if bin intersects
854 } // i
855 } // j
856 } // k
857 } // operator()
858 };
859
860 //------------------------------------------------------------------------------
861 // This code evaluates cells in intersecting bins and places them in the
862 // output list.
863 template <typename T>
FindCellsAlongPlane(const double o[3],const double n[3],double vtkNotUsed (tol),vtkIdList * cells)864 void CellProcessor<T>::FindCellsAlongPlane(
865 const double o[3], const double n[3], double vtkNotUsed(tol), vtkIdList* cells)
866 {
867 // Initialize the list of cells
868 cells->Reset();
869
870 // Make sure that the bounding box of the locator is intersected.
871 double* bounds = this->Binner->Bounds;
872 double origin[3];
873 origin[0] = o[0];
874 origin[1] = o[1];
875 origin[2] = o[2];
876 double normal[3];
877 normal[0] = n[0];
878 normal[1] = n[1];
879 normal[2] = n[2];
880 if (vtkBox::IntersectWithPlane(bounds, origin, normal) == 0)
881 {
882 return;
883 }
884
885 // Okay now evaluate which bins intersect the plane, and then the cells in
886 // the bins. We do this in parallel and mark the cells. Later the marked
887 // cells will be added (in serial) to the output list. cellHasBeenVisited
888 // has three states: 0(not visited), 1(visited but not intersecting),
889 // 2(visited and potential intersection candidate).
890 std::vector<unsigned char> cellHasBeenVisited(this->NumCells, 0);
891
892 // For now we will parallelize over z-slabs of bins.
893 CellPlaneCandidates<T> cellCandidates(this, this->Binner, o, n, cellHasBeenVisited.data());
894 vtkSMPTools::For(0, this->Binner->Divisions[2], cellCandidates);
895
896 // Populate the output list.
897 for (vtkIdType cellId = 0; cellId < this->NumCells; ++cellId)
898 {
899 if (cellHasBeenVisited[cellId] >= 2) // candidate
900 {
901 cells->InsertNextId(cellId);
902 }
903 }
904 }
905
906 //
907 //------------------------------------------------------------------------------
908 // Calculate the distance between the point x and the specified bounds
909 //
910 // WARNING!!!!! Be very careful altering this routine. Simple changes to this
911 // routine can make it 25% slower!!!!
912 // Return closest point (if any) AND the cell on which this closest point lies
Distance2ToBounds(const double x[3],double bounds[6])913 double Distance2ToBounds(const double x[3], double bounds[6])
914 {
915 double distance;
916 double deltas[3];
917
918 // Are we within the bounds?
919 if (x[0] >= bounds[0] && x[0] <= bounds[1] && x[1] >= bounds[2] && x[1] <= bounds[3] &&
920 x[2] >= bounds[4] && x[2] <= bounds[5])
921 {
922 return 0.0;
923 }
924
925 deltas[0] = deltas[1] = deltas[2] = 0.0;
926
927 // dx
928 //
929 if (x[0] < bounds[0])
930 {
931 deltas[0] = bounds[0] - x[0];
932 }
933 else if (x[0] > bounds[1])
934 {
935 deltas[0] = x[0] - bounds[1];
936 }
937
938 // dy
939 //
940 if (x[1] < bounds[2])
941 {
942 deltas[1] = bounds[2] - x[1];
943 }
944 else if (x[1] > bounds[3])
945 {
946 deltas[1] = x[1] - bounds[3];
947 }
948
949 // dz
950 //
951 if (x[2] < bounds[4])
952 {
953 deltas[2] = bounds[4] - x[2];
954 }
955 else if (x[2] > bounds[5])
956 {
957 deltas[2] = x[2] - bounds[5];
958 }
959
960 distance = vtkMath::Dot(deltas, deltas);
961 return distance;
962 }
963
964 //------------------------------------------------------------------------------
965 // Return closest point (if any) AND the cell on which this closest point lies
966 template <typename T>
FindClosestPointWithinRadius(const double x[3],double radius,double closestPoint[3],vtkGenericCell * cell,vtkIdType & closestCellId,int & closestSubId,double & minDist2,int & inside)967 vtkIdType CellProcessor<T>::FindClosestPointWithinRadius(const double x[3], double radius,
968 double closestPoint[3], vtkGenericCell* cell, vtkIdType& closestCellId, int& closestSubId,
969 double& minDist2, int& inside)
970 {
971 std::vector<bool> binHasBeenQueued(this->NumBins, false);
972 std::vector<bool> cellHasBeenVisited(this->NumCells, false);
973 std::vector<double> weights(6);
974 double pcoords[3], point[3], bds[6];
975 double distance2ToCellBounds, dist2;
976 int subId;
977 int ijk[3];
978 vtkIdType retVal = 0;
979
980 using node = std::pair<double, vtkIdType>;
981 std::priority_queue<node, std::vector<node>, std::greater<node>> queue;
982
983 // first get ijk containing point
984 vtkIdType binId = this->Binner->GetBinIndex(x);
985 queue.push(std::make_pair(0.0, binId));
986 binHasBeenQueued[binId] = true;
987
988 // distance to closest point
989 minDist2 = radius * radius;
990
991 // Process the queue of candidate bins until the candidate bins are
992 // further away than the current closest point.
993 while (!queue.empty())
994 {
995 binId = queue.top().second;
996 double binDist2 = queue.top().first;
997 queue.pop();
998
999 // stop if bounding box is further away than current closest point
1000 if (binDist2 > minDist2)
1001 {
1002 break;
1003 }
1004
1005 // compute distance to cells in bin, if any
1006 T numIds = this->GetNumberOfIds(binId);
1007 if (numIds >= 1)
1008 {
1009 const CellFragments<T>* cellIds = this->GetIds(binId);
1010 double* bounds;
1011 vtkIdType cellId;
1012
1013 for (int j = 0; j < numIds; j++)
1014 {
1015 cellId = cellIds[j].CellId;
1016
1017 // skip if cell was already visited
1018 if (cellHasBeenVisited[cellId])
1019 {
1020 continue;
1021 }
1022 cellHasBeenVisited[cellId] = true;
1023
1024 // compute distance to cell bounding box
1025 bounds = this->CellBounds + 6 * cellId;
1026 distance2ToCellBounds = Distance2ToBounds(x, bounds);
1027
1028 // compute distance to cell only if distance to bounding box smaller than minDist2
1029 if (distance2ToCellBounds < minDist2)
1030 {
1031 this->DataSet->GetCell(cellId, cell);
1032
1033 // make sure we have enough storage space for the weights
1034 unsigned nPoints = static_cast<unsigned>(cell->GetPointIds()->GetNumberOfIds());
1035 if (nPoints > weights.size())
1036 {
1037 weights.resize(2 * nPoints);
1038 }
1039
1040 // evaluate the position to find the closest point
1041 // stat==(-1) is numerical error; stat==0 means outside;
1042 // stat=1 means inside. However, for real world performance,
1043 // we sometime select stat==0 cells if the distance is close
1044 // enough
1045 int stat = cell->EvaluatePosition(x, point, subId, pcoords, dist2, weights.data());
1046
1047 if (stat != -1 && dist2 < minDist2)
1048 {
1049 retVal = 1;
1050 inside = stat;
1051 minDist2 = dist2;
1052 closestCellId = cellId;
1053 closestSubId = subId;
1054 closestPoint[0] = point[0];
1055 closestPoint[1] = point[1];
1056 closestPoint[2] = point[2];
1057 }
1058 }
1059 }
1060 }
1061
1062 // Expand the search: queue neighbors if they are not already processed. (TODO: this
1063 // is pretty inefficient with sparse locators, since bins are revisited repeatedly.)
1064 this->Binner->GetBinIndices(binId, ijk);
1065 int ijkLo[3] = { std::max(ijk[0] - 1, 0), std::max(ijk[1] - 1, 0), std::max(ijk[2] - 1, 0) };
1066 int ijkHi[3] = { std::min(ijk[0] + 1, this->Binner->Divisions[0] - 1),
1067 std::min(ijk[1] + 1, this->Binner->Divisions[1] - 1),
1068 std::min(ijk[2] + 1, this->Binner->Divisions[2] - 1) };
1069
1070 for (ijk[0] = ijkLo[0]; ijk[0] <= ijkHi[0]; ++ijk[0])
1071 {
1072 for (ijk[1] = ijkLo[1]; ijk[1] <= ijkHi[1]; ++ijk[1])
1073 {
1074 for (ijk[2] = ijkLo[2]; ijk[2] <= ijkHi[2]; ++ijk[2])
1075 {
1076 binId = this->Binner->GetBinIndex(ijk);
1077 if (!binHasBeenQueued[binId])
1078 {
1079 binHasBeenQueued[binId] = true;
1080
1081 // get bin bounding box
1082 bds[0] = this->Binner->Bounds[0] + ijk[0] * this->Binner->hX;
1083 bds[2] = this->Binner->Bounds[2] + ijk[1] * this->Binner->hY;
1084 bds[4] = this->Binner->Bounds[4] + ijk[2] * this->Binner->hZ;
1085 bds[1] = bds[0] + this->Binner->hX;
1086 bds[3] = bds[2] + this->Binner->hY;
1087 bds[5] = bds[4] + this->Binner->hZ;
1088
1089 // compute distance to box
1090 distance2ToCellBounds = Distance2ToBounds(x, bds);
1091
1092 // add to queue
1093 queue.push(std::make_pair(distance2ToCellBounds, binId));
1094 }
1095 }
1096 }
1097 }
1098 }
1099 return retVal;
1100 }
1101
1102 //------------------------------------------------------------------------------
1103 // This code traverses the cell locator by following the intersection ray. As
1104 // each bin is intersected, the cells contained in the bin are
1105 // intersected. The cell with the smallest parametric coordinate t is
1106 // returned (assuming 0<=t<=1). Otherwise no intersection is returned. See
1107 // for reference: A Fast Voxel Traversal Algorithm for Ray Tracing by John
1108 // Amanatides & Andrew Woo. Also see the code repository which inspired some
1109 // of this code:
1110 // https://github.com/francisengelmann/fast_voxel_traversal/blob/master/main.cpp.
1111 template <typename T>
IntersectWithLine(const double a0[3],const double a1[3],double tol,double & t,double x[3],double pcoords[3],int & subId,vtkIdType & cellId,vtkGenericCell * cell)1112 int CellProcessor<T>::IntersectWithLine(const double a0[3], const double a1[3], double tol,
1113 double& t, double x[3], double pcoords[3], int& subId, vtkIdType& cellId, vtkGenericCell* cell)
1114 {
1115 cellId = (-1);
1116 subId = 0;
1117
1118 double* bounds = this->Binner->Bounds;
1119 double t0, t1, x0[3], x1[3], rayDir[3];
1120 int plane0, plane1;
1121 vtkMath::Subtract(a1, a0, rayDir);
1122
1123 // Make sure the bounding box of the locator is hit.
1124 if (vtkBox::IntersectWithLine(bounds, a0, a1, t0, t1, x0, x1, plane0, plane1) == 0)
1125 {
1126 return 0;
1127 }
1128
1129 int* ndivs = this->Binner->Divisions;
1130 vtkIdType prod = ndivs[0] * ndivs[1];
1131 double* h = this->Binner->H;
1132 T i, numCellsInBin;
1133 double tMin = VTK_FLOAT_MAX;
1134 int ijk[3], ijkEnd[3];
1135 vtkIdType idx, cId, bestCellId = (-1);
1136 double hitCellBoundsPosition[3], tHitCell;
1137 double step[3], next[3], tMax[3], tDelta[3];
1138 double binBounds[6], binTol = this->Binner->binTol;
1139
1140 // Initialize intersection query array if necessary. This is done
1141 // locally to ensure thread safety.
1142 std::vector<unsigned char> cellHasBeenVisited(this->NumCells, 0);
1143
1144 // Get the i-j-k point of intersection and bin index for the entry and exit
1145 // from the locator. This is clamped to the boundary of the locator.
1146 this->Binner->GetBinIndices(x0, ijk);
1147 this->Binner->GetBinIndices(x1, ijkEnd);
1148 idx = ijk[0] + ijk[1] * ndivs[0] + ijk[2] * prod;
1149
1150 // Set up some traversal parameters for traversing through bins
1151 step[0] = (rayDir[0] >= 0.0) ? 1.0 : -1.0;
1152 step[1] = (rayDir[1] >= 0.0) ? 1.0 : -1.0;
1153 step[2] = (rayDir[2] >= 0.0) ? 1.0 : -1.0;
1154
1155 // If the ray is going in the negative direction, then the next voxel boundary
1156 // is on the "-" direction so we stay in the current voxel.
1157 next[0] = bounds[0] + h[0] * (rayDir[0] >= 0.0 ? (ijk[0] + step[0]) : ijk[0]);
1158 next[1] = bounds[2] + h[1] * (rayDir[1] >= 0.0 ? (ijk[1] + step[1]) : ijk[1]);
1159 next[2] = bounds[4] + h[2] * (rayDir[2] >= 0.0 ? (ijk[2] + step[2]) : ijk[2]);
1160
1161 tMax[0] = (rayDir[0] != 0.0) ? (next[0] - x0[0]) / rayDir[0] : VTK_FLOAT_MAX;
1162 tMax[1] = (rayDir[1] != 0.0) ? (next[1] - x0[1]) / rayDir[1] : VTK_FLOAT_MAX;
1163 tMax[2] = (rayDir[2] != 0.0) ? (next[2] - x0[2]) / rayDir[2] : VTK_FLOAT_MAX;
1164
1165 tDelta[0] = (rayDir[0] != 0.0) ? (h[0] / rayDir[0]) * step[0] : VTK_FLOAT_MAX;
1166 tDelta[1] = (rayDir[1] != 0.0) ? (h[1] / rayDir[1]) * step[1] : VTK_FLOAT_MAX;
1167 tDelta[2] = (rayDir[2] != 0.0) ? (h[2] / rayDir[2]) * step[2] : VTK_FLOAT_MAX;
1168
1169 // Start walking through the bins, find the best cell of
1170 // intersection. Note that the ray may not penetrate all of the way
1171 // through the locator so may terminate when (t > 1.0).
1172 for (bestCellId = (-1); bestCellId < 0;)
1173 {
1174 if ((numCellsInBin = this->GetNumberOfIds(idx)) > 0) // there are some cell here
1175 {
1176 const CellFragments<T>* cellIds = this->GetIds(idx);
1177 this->ComputeBinBounds(ijk[0], ijk[1], ijk[2], binBounds);
1178 for (i = 0; i < numCellsInBin; i++)
1179 {
1180 cId = cellIds[i].CellId;
1181 if (cellHasBeenVisited[cId] == 0)
1182 {
1183 cellHasBeenVisited[cId] = 1;
1184
1185 // check whether we intersect the cell bounds
1186 int hitCellBounds = vtkBox::IntersectBox(
1187 this->CellBounds + (6 * cId), a0, rayDir, hitCellBoundsPosition, tHitCell);
1188
1189 if (hitCellBounds)
1190 {
1191 // now, do the expensive GetCell call and the expensive
1192 // intersect with line call
1193 this->DataSet->GetCell(cId, cell);
1194 if (cell->IntersectWithLine(a0, a1, tol, t, x, pcoords, subId) && t < tMin)
1195 {
1196 // Make sure that intersection occurs within this bin or else spurious cell
1197 // intersections can occur behind this bin which are not the correct answer.
1198 if (!this->IsInBinBounds(binBounds, x, binTol))
1199 {
1200 cellHasBeenVisited[cId] = 0; // mark the cell non-visited
1201 }
1202 else
1203 {
1204 tMin = t;
1205 bestCellId = cId;
1206 }
1207 } // if intersection
1208 } // if (hitCellBounds)
1209 } // if (!cellHasBeenVisited[cId])
1210 } // over all cells in bin
1211 } // if cells in bin
1212
1213 // Exit before end of ray, saves a few cycles
1214 if (bestCellId >= 0 || (ijk[0] == ijkEnd[0] && ijk[1] == ijkEnd[1] && ijk[2] == ijkEnd[2]))
1215 {
1216 break;
1217 }
1218
1219 // Advance to next voxel
1220 if (tMax[0] < tMax[1])
1221 {
1222 if (tMax[0] < tMax[2])
1223 {
1224 ijk[0] += static_cast<int>(step[0]);
1225 tMax[0] += tDelta[0];
1226 }
1227 else
1228 {
1229 ijk[2] += static_cast<int>(step[2]);
1230 tMax[2] += tDelta[2];
1231 }
1232 }
1233 else
1234 {
1235 if (tMax[1] < tMax[2])
1236 {
1237 ijk[1] += static_cast<int>(step[1]);
1238 tMax[1] += tDelta[1];
1239 }
1240 else
1241 {
1242 ijk[2] += static_cast<int>(step[2]);
1243 tMax[2] += tDelta[2];
1244 }
1245 }
1246
1247 // If exiting the locator, we are done
1248 if (ijk[0] < 0 || ijk[0] >= ndivs[0] || ijk[1] < 0 || ijk[1] >= ndivs[1] || ijk[2] < 0 ||
1249 ijk[2] >= ndivs[2])
1250 {
1251 break;
1252 }
1253 else
1254 {
1255 idx = ijk[0] + ijk[1] * ndivs[0] + ijk[2] * prod;
1256 }
1257
1258 } // for looking for valid intersected cell
1259
1260 // If a cell has been intersected, recover the information and return.
1261 // This information could be cached....
1262 if (bestCellId >= 0)
1263 {
1264 this->DataSet->GetCell(bestCellId, cell);
1265 cell->IntersectWithLine(a0, a1, tol, t, x, pcoords, subId);
1266
1267 // store the best cell id in the return "parameter"
1268 cellId = bestCellId;
1269 return 1;
1270 }
1271
1272 return 0;
1273 }
1274
1275 } // anonymous namespace
1276
1277 //------------------------------------------------------------------------------
1278 // Here is the VTK class proper.
1279
1280 //------------------------------------------------------------------------------
vtkStaticCellLocator()1281 vtkStaticCellLocator::vtkStaticCellLocator()
1282 {
1283 this->CacheCellBounds = 1; // always cached
1284
1285 this->Binner = nullptr;
1286 this->Processor = nullptr;
1287
1288 this->NumberOfCellsPerNode = 10;
1289 this->Divisions[0] = this->Divisions[1] = this->Divisions[2] = 100;
1290 this->H[0] = this->H[1] = this->H[2] = 0.0;
1291 for (int i = 0; i < 6; i++)
1292 {
1293 this->Bounds[i] = 0;
1294 }
1295
1296 this->MaxNumberOfBuckets = VTK_INT_MAX;
1297 this->LargeIds = false;
1298 }
1299
1300 //------------------------------------------------------------------------------
~vtkStaticCellLocator()1301 vtkStaticCellLocator::~vtkStaticCellLocator()
1302 {
1303 this->FreeSearchStructure();
1304 }
1305
1306 //------------------------------------------------------------------------------
FreeSearchStructure()1307 void vtkStaticCellLocator::FreeSearchStructure()
1308 {
1309 if (this->Binner)
1310 {
1311 delete this->Binner;
1312 this->Binner = nullptr;
1313 }
1314 if (this->Processor)
1315 {
1316 delete this->Processor;
1317 this->Processor = nullptr;
1318 }
1319 }
1320
1321 //------------------------------------------------------------------------------
FindCell(double pos[3],double,vtkGenericCell * cell,double pcoords[3],double * weights)1322 vtkIdType vtkStaticCellLocator::FindCell(
1323 double pos[3], double, vtkGenericCell* cell, double pcoords[3], double* weights)
1324 {
1325 this->BuildLocator();
1326 if (!this->Processor)
1327 {
1328 return -1;
1329 }
1330 return this->Processor->FindCell(pos, cell, pcoords, weights);
1331 }
1332
1333 //------------------------------------------------------------------------------
FindClosestPoint(const double x[3],double closestPoint[3],vtkGenericCell * cell,vtkIdType & cellId,int & subId,double & dist2)1334 void vtkStaticCellLocator::FindClosestPoint(const double x[3], double closestPoint[3],
1335 vtkGenericCell* cell, vtkIdType& cellId, int& subId, double& dist2)
1336 {
1337 int inside;
1338 double radius = vtkMath::Inf();
1339 double point[3] = { x[0], x[1], x[2] };
1340 this->FindClosestPointWithinRadius(
1341 point, radius, closestPoint, cell, cellId, subId, dist2, inside);
1342 }
1343
1344 //------------------------------------------------------------------------------
FindClosestPointWithinRadius(double x[3],double radius,double closestPoint[3],vtkGenericCell * cell,vtkIdType & cellId,int & subId,double & dist2,int & inside)1345 vtkIdType vtkStaticCellLocator::FindClosestPointWithinRadius(double x[3], double radius,
1346 double closestPoint[3], vtkGenericCell* cell, vtkIdType& cellId, int& subId, double& dist2,
1347 int& inside)
1348 {
1349 this->BuildLocator();
1350 if (!this->Processor)
1351 {
1352 return 0;
1353 }
1354 return this->Processor->FindClosestPointWithinRadius(
1355 x, radius, closestPoint, cell, cellId, subId, dist2, inside);
1356 }
1357
1358 //------------------------------------------------------------------------------
FindCellsWithinBounds(double * bbox,vtkIdList * cells)1359 void vtkStaticCellLocator::FindCellsWithinBounds(double* bbox, vtkIdList* cells)
1360 {
1361 this->BuildLocator();
1362 if (!this->Processor)
1363 {
1364 return;
1365 }
1366 return this->Processor->FindCellsWithinBounds(bbox, cells);
1367 }
1368
1369 //------------------------------------------------------------------------------
FindCellsAlongLine(const double p1[3],const double p2[3],double tol,vtkIdList * cells)1370 void vtkStaticCellLocator::FindCellsAlongLine(
1371 const double p1[3], const double p2[3], double tol, vtkIdList* cells)
1372 {
1373 this->BuildLocator();
1374 if (!this->Processor)
1375 {
1376 return;
1377 }
1378 return this->Processor->FindCellsAlongLine(p1, p2, tol, cells);
1379 }
1380
1381 //------------------------------------------------------------------------------
FindCellsAlongPlane(const double o[3],const double n[3],double tol,vtkIdList * cells)1382 void vtkStaticCellLocator::FindCellsAlongPlane(
1383 const double o[3], const double n[3], double tol, vtkIdList* cells)
1384 {
1385 this->BuildLocator();
1386 if (!this->Processor)
1387 {
1388 return;
1389 }
1390 return this->Processor->FindCellsAlongPlane(o, n, tol, cells);
1391 }
1392
1393 //------------------------------------------------------------------------------
IntersectWithLine(const double p1[3],const double p2[3],double tol,double & t,double x[3],double pcoords[3],int & subId,vtkIdType & cellId,vtkGenericCell * cell)1394 int vtkStaticCellLocator::IntersectWithLine(const double p1[3], const double p2[3], double tol,
1395 double& t, double x[3], double pcoords[3], int& subId, vtkIdType& cellId, vtkGenericCell* cell)
1396 {
1397 this->BuildLocator();
1398 if (!this->Processor)
1399 {
1400 return 0;
1401 }
1402 return this->Processor->IntersectWithLine(p1, p2, tol, t, x, pcoords, subId, cellId, cell);
1403 }
1404
1405 //------------------------------------------------------------------------------
BuildLocator()1406 void vtkStaticCellLocator::BuildLocator()
1407 {
1408 vtkDebugMacro(<< "Building static cell locator");
1409
1410 // Do we need to build?
1411 if ((this->Binner != nullptr) && (this->BuildTime > this->MTime) &&
1412 (this->BuildTime > this->DataSet->GetMTime()))
1413 {
1414 return;
1415 }
1416
1417 vtkIdType numCells;
1418 if (!this->DataSet || (numCells = this->DataSet->GetNumberOfCells()) < 1)
1419 {
1420 vtkErrorMacro(<< "No cells to build");
1421 return;
1422 }
1423
1424 // Prepare
1425 if (this->Binner)
1426 {
1427 this->FreeSearchStructure();
1428 }
1429
1430 // The bounding box can be slow
1431 int i, ndivs[3];
1432 const double* bounds = this->DataSet->GetBounds();
1433 vtkIdType numBins = static_cast<vtkIdType>(
1434 static_cast<double>(numCells) / static_cast<double>(this->NumberOfCellsPerNode));
1435 numBins = (numBins > this->MaxNumberOfBuckets ? this->MaxNumberOfBuckets : numBins);
1436
1437 vtkBoundingBox bbox(bounds);
1438 if (this->Automatic)
1439 {
1440 bbox.ComputeDivisions(numBins, this->Bounds, ndivs);
1441 }
1442 else
1443 {
1444 bbox.Inflate(); // make sure non-zero volume
1445 bbox.GetBounds(this->Bounds);
1446 for (i = 0; i < 3; i++)
1447 {
1448 ndivs[i] = (this->Divisions[i] < 1 ? 1 : this->Divisions[i]);
1449 }
1450 }
1451
1452 this->Divisions[0] = ndivs[0];
1453 this->Divisions[1] = ndivs[1];
1454 this->Divisions[2] = ndivs[2];
1455 numBins = static_cast<vtkIdType>(ndivs[0]) * static_cast<vtkIdType>(ndivs[1]) *
1456 static_cast<vtkIdType>(ndivs[2]);
1457
1458 // Compute bin/bucket widths
1459 for (i = 0; i < 3; i++)
1460 {
1461 this->H[i] = (this->Bounds[2 * i + 1] - this->Bounds[2 * i]) / this->Divisions[i];
1462 }
1463
1464 // Actually do the hard work of creating the locator. Clear out old stuff.
1465 delete this->Binner;
1466 delete this->Processor;
1467 this->Binner = new vtkCellBinner(this, numCells, numBins);
1468 vtkSMPTools::For(0, numCells, *(this->Binner));
1469
1470 // Create sorted cell fragments tuples of (cellId,binId). Depending
1471 // on problem size, different types are used.
1472 vtkIdType numFragments = this->Binner->NumFragments;
1473 if (numFragments >= VTK_INT_MAX)
1474 {
1475 this->LargeIds = true;
1476 CellProcessor<vtkIdType>* processor = new CellProcessor<vtkIdType>(this->Binner);
1477 vtkSMPTools::For(0, numCells, *processor);
1478 vtkSMPTools::Sort(processor->Map, processor->Map + numFragments);
1479 MapOffsets<vtkIdType> mapOffsets(processor);
1480 vtkSMPTools::For(0, processor->NumBatches, mapOffsets);
1481 this->Processor = processor;
1482 }
1483 else
1484 {
1485 this->LargeIds = false;
1486 CellProcessor<int>* processor = new CellProcessor<int>(this->Binner);
1487 vtkSMPTools::For(0, numCells, *processor);
1488 vtkSMPTools::Sort(processor->Map, processor->Map + numFragments);
1489 MapOffsets<int> mapOffsets(processor);
1490 vtkSMPTools::For(0, processor->NumBatches, mapOffsets);
1491 this->Processor = processor;
1492 }
1493
1494 this->BuildTime.Modified();
1495 }
1496
1497 //------------------------------------------------------------------------------
1498 // Produce a polygonal representation of the locator. Each bin which contains
1499 // a potential cell candidate contributes to the representation. Note that
1500 // since the locator has only a single level, the level method parameter is
1501 // ignored.
GenerateRepresentation(int vtkNotUsed (level),vtkPolyData * pd)1502 void vtkStaticCellLocator::GenerateRepresentation(int vtkNotUsed(level), vtkPolyData* pd)
1503 {
1504 // Make sure locator has been built successfully
1505 this->BuildLocator();
1506 if (!this->Processor)
1507 {
1508 return;
1509 }
1510
1511 vtkPoints* pts = vtkPoints::New();
1512 pts->SetDataTypeToFloat();
1513 vtkCellArray* polys = vtkCellArray::New();
1514 pd->SetPoints(pts);
1515 pd->SetPolys(polys);
1516
1517 int* dims = this->Divisions;
1518 int i, j, k;
1519 vtkIdType idx, kOffset, jOffset, kSlice = dims[0] * dims[1], pIds[8];
1520 double* s = this->H;
1521 double x[3], xT[3], origin[3];
1522 origin[0] = this->Bounds[0];
1523 origin[1] = this->Bounds[2];
1524 origin[2] = this->Bounds[4];
1525
1526 // A locator is used to avoid duplicate points
1527 vtkMergePoints* locator = vtkMergePoints::New();
1528 locator->InitPointInsertion(pts, this->Bounds, dims[0] * dims[1] * dims[2]);
1529
1530 for (k = 0; k < dims[2]; k++)
1531 {
1532 x[2] = origin[2] + k * s[2];
1533 kOffset = k * kSlice;
1534 for (j = 0; j < dims[1]; j++)
1535 {
1536 x[1] = origin[1] + j * s[1];
1537 jOffset = j * dims[0];
1538 for (i = 0; i < dims[0]; i++)
1539 {
1540 x[0] = origin[0] + i * s[0];
1541 idx = i + jOffset + kOffset;
1542
1543 // Check to see if bin contains anything
1544 // If so, insert up to eight points and six quad faces (depending on
1545 // local topology).
1546 if (!this->Processor->IsEmpty(idx))
1547 {
1548 // Points in (i-j-k) order. A locator is used to avoid duplicate points.
1549 locator->InsertUniquePoint(x, pIds[0]);
1550 xT[0] = x[0] + s[0];
1551 xT[1] = x[1];
1552 xT[2] = x[2];
1553 locator->InsertUniquePoint(xT, pIds[1]);
1554 xT[0] = x[0];
1555 xT[1] = x[1] + s[1];
1556 xT[2] = x[2];
1557 locator->InsertUniquePoint(xT, pIds[2]);
1558 xT[0] = x[0] + s[0];
1559 xT[1] = x[1] + s[1];
1560 xT[2] = x[2];
1561 locator->InsertUniquePoint(xT, pIds[3]);
1562 xT[0] = x[0];
1563 xT[1] = x[1];
1564 xT[2] = x[2] + s[2];
1565 locator->InsertUniquePoint(xT, pIds[4]);
1566 xT[0] = x[0] + s[0];
1567 xT[1] = x[1];
1568 xT[2] = x[2] + s[2];
1569 locator->InsertUniquePoint(xT, pIds[5]);
1570 xT[0] = x[0];
1571 xT[1] = x[1] + s[1];
1572 xT[2] = x[2] + s[2];
1573 locator->InsertUniquePoint(xT, pIds[6]);
1574 xT[0] = x[0] + s[0];
1575 xT[1] = x[1] + s[1];
1576 xT[2] = x[2] + s[2];
1577 locator->InsertUniquePoint(xT, pIds[7]);
1578
1579 // Loop over all bins. Any bin containing cell candidates may
1580 // generate output. Faces are output if they are on the boundary of
1581 // the locator or if the bin neighbor contains no cells (i.e., there
1582 // are no face neighbors). This prevents duplicate faces.
1583
1584 // -x bin boundary face
1585 if (i == 0 || this->Processor->IsEmpty(idx - 1))
1586 {
1587 polys->InsertNextCell(4);
1588 polys->InsertCellPoint(pIds[0]);
1589 polys->InsertCellPoint(pIds[4]);
1590 polys->InsertCellPoint(pIds[6]);
1591 polys->InsertCellPoint(pIds[2]);
1592 }
1593
1594 // +x boundary face
1595 if (i == (dims[0] - 1) || this->Processor->IsEmpty(idx + 1))
1596 {
1597 polys->InsertNextCell(4);
1598 polys->InsertCellPoint(pIds[1]);
1599 polys->InsertCellPoint(pIds[3]);
1600 polys->InsertCellPoint(pIds[7]);
1601 polys->InsertCellPoint(pIds[5]);
1602 }
1603
1604 // -y boundary face
1605 if (j == 0 || this->Processor->IsEmpty(idx - dims[0]))
1606 {
1607 polys->InsertNextCell(4);
1608 polys->InsertCellPoint(pIds[0]);
1609 polys->InsertCellPoint(pIds[1]);
1610 polys->InsertCellPoint(pIds[5]);
1611 polys->InsertCellPoint(pIds[4]);
1612 }
1613
1614 // +y boundary face
1615 if (j == (dims[1] - 1) || this->Processor->IsEmpty(idx + dims[0]))
1616 {
1617 polys->InsertNextCell(4);
1618 polys->InsertCellPoint(pIds[2]);
1619 polys->InsertCellPoint(pIds[6]);
1620 polys->InsertCellPoint(pIds[7]);
1621 polys->InsertCellPoint(pIds[3]);
1622 }
1623
1624 // -z boundary face
1625 if (k == 0 || this->Processor->IsEmpty(idx - kSlice))
1626 {
1627 polys->InsertNextCell(4);
1628 polys->InsertCellPoint(pIds[0]);
1629 polys->InsertCellPoint(pIds[2]);
1630 polys->InsertCellPoint(pIds[3]);
1631 polys->InsertCellPoint(pIds[1]);
1632 }
1633
1634 // +z boundary face
1635 if (k == (dims[2] - 1) || this->Processor->IsEmpty(idx + kSlice))
1636 {
1637 polys->InsertNextCell(4);
1638 polys->InsertCellPoint(pIds[4]);
1639 polys->InsertCellPoint(pIds[5]);
1640 polys->InsertCellPoint(pIds[7]);
1641 polys->InsertCellPoint(pIds[6]);
1642 }
1643 } // if not empty
1644 } // x
1645 } // y
1646 } // z
1647
1648 // Clean up
1649 locator->Delete();
1650 polys->Delete();
1651 pts->Delete();
1652 }
1653
1654 //------------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)1655 void vtkStaticCellLocator::PrintSelf(ostream& os, vtkIndent indent)
1656 {
1657 // Cell bounds are always cached
1658 this->CacheCellBounds = 1;
1659 this->Superclass::PrintSelf(os, indent);
1660
1661 os << indent << "Max Number Of Buckets: " << this->MaxNumberOfBuckets << "\n";
1662
1663 os << indent << "Large IDs: " << this->LargeIds << "\n";
1664 }
1665