Method FindNPointsInShell

Common/DataModel/vtkStaticPointLocator.cxx:937–1039

// are processed.
template <typename TIds>
double BucketList<TIds>::FindNPointsInShell(int N, const double x[3], vtkDist2TupleArray& results,
  double minR2, bool sort, vtkDoubleArray* spheres)
{
  // Clear out any previous results.
  results.clear();

  // Find the bucket/bin the point is in. This is the center of the request
  // footprint.
  int center[3];
  this->GetBucketIndices(x, center);

  // Traverse and gather points in the bucket/bins contained in the shell
  // request (minR,maxR].
  double minR = sqrt(minR2);
  double maxR = 0, maxR2 = 0;

  // Determine absolute limits of iteration (based on possible number of points).
  vtkIdType numPts = this->DataSet->GetNumberOfPoints();
  N = (numPts < N ? numPts : N);

  // Gather N points if possible and determine maxR2. Make sure all points within maxR2
  // have been found. We use a shell iterator to grow a "rectangular" shell from the center bin.
  // Skip over bins inside the inner radius minR2. The bin index is updated during iteration,
  // a binIdx<0 is returned when the iteration is exhausted.
  vtkIdType binIdx;
  int currentLevel = (minR <= 0 ? 0 : static_cast<int>(std::floor(minR / (2.0 * this->BinRadius))));
  ShellIterator<TIds> siter(
    this->DataSet, this, this->Divisions, this->BinRadius, spheres, x, center);

  // Loop across levels, accruing points as we go. This will determine the maxR2. It also
  // carves out some inner levels (based on current level) that do not have to be revisited.
  while (static_cast<int>(results.size()) < N && currentLevel < this->MaxLevel)
  {
    int i, j, k;
    binIdx = siter.Initialize(currentLevel, i, j, k);
    // Basically iterating over a "rectangular" footprint defined from
    // the current level.
    while (binIdx >= 0)
    {
      maxR2 = siter.GatherPoints(i, j, k, binIdx, currentLevel, N, minR2, maxR2, results);
      binIdx = siter.NextBin(i, j, k);
    }
    ++currentLevel;
  }
  int level = currentLevel - 1; // reset to the last level processed

  // We have determined maxR2 and ~N points in the request annulus
  // (minR2 < p_d2 <= maxR2). Now gather any other remaining points
  // within the request. It's typical that the number of points
  // returned is >N.
  maxR = sqrt(maxR2);

  // Determine the range of indices in each direction based on radius maxR.
  // This block of bins is processed to gather any additional points with
  // radius <=maxR2.
  double xMin[3], xMax[3];
  xMin[0] = x[0] - maxR;
  xMin[1] = x[1] - maxR;