/*=========================================================================

  Program:   Visualization Toolkit
  Module:    vtkPointLocator.cxx

  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.

     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
     PURPOSE.  See the above copyright notice for more information.

=========================================================================*/
#include "vtkPointLocator.h"

#include "vtkBoundingBox.h"
#include "vtkCellArray.h"
#include "vtkIdList.h"
#include "vtkIntArray.h"
#include "vtkMath.h"
#include "vtkObjectFactory.h"
#include "vtkPolyData.h"

#include <algorithm> //std::sort

vtkStandardNewMacro(vtkPointLocator);

static const int VTK_INITIAL_SIZE = 1000;

//------------------------------------------------------------------------------
// Utility class to store an array of ijk values
class vtkNeighborPoints
{
public:
  vtkNeighborPoints()
  {
    this->Count = 0;
    this->P = this->InitialBuffer;
    this->MaxSize = VTK_INITIAL_SIZE;
  }
  ~vtkNeighborPoints()
  {
    this->Count = 0;
    if (this->P != this->InitialBuffer)
    {
      delete[] this->P;
    }
  }
  int GetNumberOfNeighbors() { return this->Count; }
  void Reset() { this->Count = 0; }

  int* GetPoint(int i) { return (this->Count > i ? (this->P + 3 * i) : nullptr); }

  int InsertNextPoint(const int x[3])
  {
    int* tmp;

    // Re-allocate if beyond the current max size.
    // (Increase by VTK_INITIAL_SIZE)
    if (this->Count == this->MaxSize)
    {
      tmp = this->P;

      this->MaxSize += VTK_INITIAL_SIZE;
      this->P = new int[this->MaxSize * 3];

      for (int i = 0; i < 3 * this->Count; i++)
      {
        this->P[i] = tmp[i];
      }
      if (tmp != this->InitialBuffer)
      {
        delete[] tmp;
      }
    }

    this->P[3 * this->Count] = x[0];
    this->P[3 * this->Count + 1] = x[1];
    this->P[3 * this->Count + 2] = x[2];
    this->Count++;
    return this->Count - 1;
  }

protected:
  // Start with an array to avoid memory allocation overhead
  int InitialBuffer[VTK_INITIAL_SIZE * 3];
  int* P;
  int Count;
  int MaxSize;
};

//------------------------------------------------------------------------------
// Construct with automatic computation of divisions, averaging
// 3 points per bucket.
vtkPointLocator::vtkPointLocator()
{
  this->Points = nullptr;
  this->Divisions[0] = this->Divisions[1] = this->Divisions[2] = 50;
  this->NumberOfPointsPerBucket = 3;
  this->HashTable = nullptr;
  this->H[0] = this->H[1] = this->H[2] = 0.0;
  this->InsertionPointId = 0;
  this->InsertionTol2 = 0.0001;
  this->InsertionLevel = 0;
}

//------------------------------------------------------------------------------
vtkPointLocator::~vtkPointLocator()
{
  if (this->Points)
  {
    this->Points->UnRegister(this);
    this->Points = nullptr;
  }
  this->FreeSearchStructure();
}

//------------------------------------------------------------------------------
void vtkPointLocator::Initialize()
{
  if (this->Points)
  {
    this->Points->UnRegister(this);
    this->Points = nullptr;
  }
  this->FreeSearchStructure();
}

//------------------------------------------------------------------------------
void vtkPointLocator::FreeSearchStructure()
{
  vtkIdList* ptIds;
  vtkIdType i;

  if (this->HashTable)
  {
    for (i = 0; i < this->NumberOfBuckets; i++)
    {
      if ((ptIds = this->HashTable[i]))
      {
        ptIds->Delete();
      }
    }
    delete[] this->HashTable;
    this->HashTable = nullptr;
  }
}

//------------------------------------------------------------------------------
// Initialize internal performance variables
void vtkPointLocator::ComputePerformanceFactors()
{
  this->HX = this->H[0];
  this->HY = this->H[1];
  this->HZ = this->H[2];
  this->FX = 1.0 / this->H[0];
  this->FY = 1.0 / this->H[1];
  this->FZ = 1.0 / this->H[2];
  this->BX = this->Bounds[0];
  this->BY = this->Bounds[2];
  this->BZ = this->Bounds[4];
  this->XD = this->Divisions[0];
  this->YD = this->Divisions[1];
  this->ZD = this->Divisions[2];
  this->SliceSize = this->Divisions[0] * this->Divisions[1];
}

//------------------------------------------------------------------------------
// Given a position x, return the id of the point closest to it.
vtkIdType vtkPointLocator::FindClosestPoint(const double x[3])
{
  this->BuildLocator();
  if (this->HashTable == nullptr)
  {
    return -1;
  }
  int i, j;
  double minDist2;
  double dist2 = VTK_DOUBLE_MAX;
  double pt[3];
  int level;
  vtkIdType ptId, closest, cno, nids;
  vtkIdList* ptIds;
  int ijk[3], *nei;
  vtkNeighborPoints buckets;

  //  Find bucket point is in.
  this->GetBucketIndices(x, ijk);

  //
  //  Need to search this bucket for closest point.  If there are no
  //  points in this bucket, search 1st level neighbors, and so on,
  //  until closest point found.
  //
  for (closest = (-1), minDist2 = VTK_DOUBLE_MAX, level = 0; (closest == -1) &&
       (level < this->Divisions[0] || level < this->Divisions[1] || level < this->Divisions[2]);
       level++)
  {
    this->GetBucketNeighbors(&buckets, ijk, this->Divisions, level);

    for (i = 0; i < buckets.GetNumberOfNeighbors(); i++)
    {
      nei = buckets.GetPoint(i);
      cno = nei[0] + nei[1] * this->XD + nei[2] * this->SliceSize;

      if ((ptIds = this->HashTable[cno]) != nullptr)
      {
        nids = ptIds->GetNumberOfIds();
        for (j = 0; j < nids; j++)
        {
          ptId = ptIds->GetId(j);
          this->DataSet->GetPoint(ptId, pt);
          if ((dist2 = vtkMath::Distance2BetweenPoints(x, pt)) < minDist2)
          {
            closest = ptId;
            minDist2 = dist2;
          }
        }
      }
    }
  }
  //
  // Because of the relative location of the points in the buckets, the
  // point found previously may not be the closest point.  Have to
  // search those bucket neighbors that might also contain point.
  //
  if (minDist2 > 0.0)
  {
    this->GetOverlappingBuckets(&buckets, x, ijk, sqrt(minDist2), 0);
    for (i = 0; i < buckets.GetNumberOfNeighbors(); i++)
    {
      nei = buckets.GetPoint(i);
      cno = nei[0] + nei[1] * this->XD + nei[2] * this->SliceSize;

      if ((ptIds = this->HashTable[cno]) != nullptr)
      {
        nids = ptIds->GetNumberOfIds();
        for (j = 0; j < nids; j++)
        {
          ptId = ptIds->GetId(j);
          this->DataSet->GetPoint(ptId, pt);
          if ((dist2 = vtkMath::Distance2BetweenPoints(x, pt)) < minDist2)
          {
            closest = ptId;
            minDist2 = dist2;
          }
        } // for each point
      }   // if points in bucket
    }     // for each overlapping bucket
  }       // if not identical point

  return closest;
}

//------------------------------------------------------------------------------
vtkIdType vtkPointLocator::FindClosestPointWithinRadius(
  double radius, const double x[3], double& dist2)
{
  return this->FindClosestPointWithinRadius(radius, x, this->DataSet->GetLength(), dist2);
}

//------------------------------------------------------------------------------
vtkIdType vtkPointLocator::FindClosestPointWithinRadius(
  double radius, const double x[3], double inputDataLength, double& dist2)
{
  this->BuildLocator();
  if (this->HashTable == nullptr)
  {
    return -1;
  }
  int i, j;
  double pt[3];
  vtkIdType ptId, nids, closest = -1;
  vtkIdList* ptIds;
  int ijk[3], *nei;
  double minDist2;

  double refinedRadius, radius2, refinedRadius2;
  double currentRadius;
  double distance2ToDataBounds, maxDistance;
  int ii, radiusLevels[3], radiusLevel, prevMinLevel[3], prevMaxLevel[3];
  vtkNeighborPoints buckets;

  dist2 = -1.0;
  radius2 = radius * radius;
  minDist2 = 1.01 * radius2; // something slightly bigger....

  vtkPointSet* pointSet = static_cast<vtkPointSet*>(this->DataSet);
  if (!pointSet)
  {
    return closest;
  }
  vtkPoints* points = pointSet->GetPoints();
  if (!points)
  {
    return closest;
  }
  vtkDataArray* pointData = points->GetData();

  //
  //  Find bucket point is in.
  //
  this->GetBucketIndices(x, ijk);

  // Start by searching the bucket that the point is in.
  //
  if ((ptIds = this->HashTable[ijk[0] + ijk[1] * this->XD + ijk[2] * this->SliceSize]) != nullptr)
  {
    nids = ptIds->GetNumberOfIds();
    for (j = 0; j < nids; j++)
    {
      ptId = ptIds->GetId(j);
      pointData->GetTuple(ptId, pt);
      if ((dist2 = vtkMath::Distance2BetweenPoints(x, pt)) < minDist2)
      {
        closest = ptId;
        minDist2 = dist2;
      }
    }
  }

  // Now, search only those buckets that are within a radius. The radius used
  // is the smaller of sqrt(dist2) and the radius that is passed in. To avoid
  // checking a large number of buckets unnecessarily, if the radius is
  // larger than the dimensions of a bucket, we search outward using a
  // simple heuristic of rings.  This heuristic ends up collecting inner
  // buckets multiple times, but this only happens in the case where these
  // buckets are empty, so they are discarded quickly.
  //
  if (minDist2 < radius2)
  {
    refinedRadius = sqrt(dist2);
    refinedRadius2 = dist2;
  }
  else
  {
    refinedRadius = radius;
    refinedRadius2 = radius2;
  }

  if (inputDataLength != 0.0)
  {
    distance2ToDataBounds = this->Distance2ToBounds(x, this->Bounds);
    maxDistance = sqrt(distance2ToDataBounds) + inputDataLength;
    if (refinedRadius > maxDistance)
    {
      refinedRadius = maxDistance;
      refinedRadius2 = maxDistance * maxDistance;
    }
  }

  for (i = 0; i < 3; i++)
  {
    radiusLevels[i] = static_cast<int>(refinedRadius / this->H[i]);
    if (radiusLevels[i] > this->Divisions[i] / 2)
    {
      radiusLevels[i] = this->Divisions[i] / 2;
    }
  }

  radiusLevel = radiusLevels[0];
  radiusLevel = radiusLevels[1] > radiusLevel ? radiusLevels[1] : radiusLevel;
  radiusLevel = radiusLevels[2] > radiusLevel ? radiusLevels[2] : radiusLevel;
  if (radiusLevel == 0)
  {
    radiusLevel = 1;
  }

  // radius schedule increases the radius each iteration, this is currently
  // implemented by decreasing ii by 1 each iteration.  another alternative
  // is to double the radius each iteration, i.e. ii = ii >> 1
  // In practice, reducing ii by one has been found to be more efficient.
  prevMinLevel[0] = prevMaxLevel[0] = ijk[0];
  prevMinLevel[1] = prevMaxLevel[1] = ijk[1];
  prevMinLevel[2] = prevMaxLevel[2] = ijk[2];
  for (ii = radiusLevel; ii >= 1; ii--)
  {
    currentRadius = refinedRadius; // used in if at bottom of this for loop

    // Build up a list of buckets that are arranged in rings
    this->GetOverlappingBuckets(&buckets, x, refinedRadius / ii, prevMinLevel, prevMaxLevel);

    for (i = 0; i < buckets.GetNumberOfNeighbors(); i++)
    {
      nei = buckets.GetPoint(i);

      // do we still need to test this bucket?
      if (this->Distance2ToBucket(x, nei) < refinedRadius2)
      {
        ptIds = this->HashTable[nei[0] + nei[1] * this->XD + nei[2] * this->SliceSize];

        nids = ptIds->GetNumberOfIds();
        for (j = 0; j < nids; j++)
        {
          ptId = ptIds->GetId(j);
          pointData->GetTuple(ptId, pt);
          if ((dist2 = vtkMath::Distance2BetweenPoints(x, pt)) < minDist2)
          {
            closest = ptId;
            minDist2 = dist2;
            refinedRadius = sqrt(minDist2);
            refinedRadius2 = minDist2;
          }
        } // for each pt in bucket
      }   // if bucket is within the current best distance
    }     // for each overlapping bucket

    // don't want to checker a smaller radius than we just checked so update
    // ii appropriately
    if (refinedRadius < currentRadius && ii > 2) // always check ii==1
    {
      ii = static_cast<int>(static_cast<double>(ii) * (refinedRadius / currentRadius)) + 1;
      if (ii < 2)
      {
        ii = 2;
      }
    }
  } // for each radius in the radius schedule

  if ((closest != -1) && (minDist2 <= radius2))
  {
    dist2 = minDist2;
  }
  else
  {
    closest = -1;
  }

  return closest;
}

namespace
{
//------------------------------------------------------------------------------
// Sorting closest points
class IdTuple
{
public:
  vtkIdType PtId;
  double Dist2;

  bool operator<(const IdTuple& tuple) const { return Dist2 < tuple.Dist2; }
};
}

//------------------------------------------------------------------------------
void vtkPointLocator::FindDistributedPoints(
  int N, double x, double y, double z, vtkIdList* result, int M)
{
  double p[3];
  p[0] = x;
  p[1] = y;
  p[2] = z;
  this->FindDistributedPoints(N, p, result, M);
}

//------------------------------------------------------------------------------
static int GetOctant(const double x[3], const double pt[3])
{
  double tmp[3];
  int res = 0;

  tmp[0] = pt[0] - x[0];
  tmp[1] = pt[1] - x[1];
  tmp[2] = pt[2] - x[2];

  if (tmp[0] > 0.0)
  {
    res += 1;
  }
  if (tmp[1] > 0.0)
  {
    res += 2;
  }
  if (tmp[2] > 0.0)
  {
    res += 4;
  }

  return res;
}

//------------------------------------------------------------------------------
static int GetMin(const int foo[8])
{
  int result = foo[0];
  int i;

  for (i = 1; i < 8; i++)
  {
    if (foo[i] < result)
    {
      result = foo[i];
    }
  }
  return result;
}

//------------------------------------------------------------------------------
static double GetMax(const double foo[8])
{
  double result = foo[0];
  int i;

  for (i = 1; i < 8; i++)
  {
    if (foo[i] > result)
    {
      result = foo[i];
    }
  }
  return result;
}

//------------------------------------------------------------------------------
void vtkPointLocator::FindDistributedPoints(int N, const double x[3], vtkIdList* result, int M)
{
  this->BuildLocator();
  if (this->HashTable == nullptr)
  {
    return;
  }
  int i, j;
  double dist2;
  double pt[3];
  int level;
  vtkIdType ptId, nids, cno;
  vtkIdList* ptIds;
  int ijk[3], *nei;
  int oct;
  int pointsChecked = 0;
  vtkNeighborPoints buckets;

  // clear out the result
  result->Reset();
  //
  //  Make sure candidate point is in bounds.  If not, it is outside.
  //
  for (i = 0; i < 3; i++)
  {
    if (x[i] < this->Bounds[2 * i] || x[i] > this->Bounds[2 * i + 1])
    {
      return;
    }
  }
  //
  //  Find bucket point is in.
  //
  this->GetBucketIndices(x, ijk);

  // there are two steps, first a simple expanding wave of buckets until
  // we have enough points. Then a refinement to make sure we have the
  // N closest points.
  level = 0;
  double maxDistance[8] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
  int currentCount[8] = { 0, 0, 0, 0, 0, 0, 0, 0 };
  int minCurrentCount = 0;

  IdTuple* res[8];
  for (i = 0; i < 8; i++)
  {
    res[i] = new IdTuple[N];
  }

  this->GetBucketNeighbors(&buckets, ijk, this->Divisions, level);
  while (buckets.GetNumberOfNeighbors() && minCurrentCount < N && pointsChecked < M)
  {
    for (i = 0; i < buckets.GetNumberOfNeighbors(); i++)
    {
      nei = buckets.GetPoint(i);
      cno = nei[0] + nei[1] * this->XD + nei[2] * this->SliceSize;

      if ((ptIds = this->HashTable[cno]) != nullptr)
      {
        nids = ptIds->GetNumberOfIds();
        for (j = 0; j < nids; j++)
        {
          pointsChecked++;
          ptId = ptIds->GetId(j);
          this->DataSet->GetPoint(ptId, pt);
          dist2 = vtkMath::Distance2BetweenPoints(x, pt);
          oct = GetOctant(x, pt);
          if (currentCount[oct] < N)
          {
            res[oct][currentCount[oct]].Dist2 = dist2;
            res[oct][currentCount[oct]].PtId = ptId;
            if (dist2 > maxDistance[oct])
            {
              maxDistance[oct] = dist2;
            }
            currentCount[oct] = currentCount[oct] + 1;
            // compute new minCurrentCount
            minCurrentCount = GetMin(currentCount);
            if (currentCount[oct] == N)
            {
              std::sort(res[oct], res[oct] + currentCount[oct]);
            }
          }
          else if (dist2 < maxDistance[oct])
          {
            res[oct][N - 1].Dist2 = dist2;
            res[oct][N - 1].PtId = ptId;
            std::sort(res[oct], res[oct] + N);
            maxDistance[oct] = res[oct][N - 1].Dist2;
          }
        }
      }
    }
    level++;
    this->GetBucketNeighbors(&buckets, ijk, this->Divisions, level);
  }

  // do a sort
  for (i = 0; i < 8; i++)
  {
    std::sort(res[i], res[i] + currentCount[i]);
  }

  // Now do the refinement
  this->GetOverlappingBuckets(&buckets, x, ijk, sqrt(GetMax(maxDistance)), level - 1);

  for (i = 0; pointsChecked < M && i < buckets.GetNumberOfNeighbors(); i++)
  {
    nei = buckets.GetPoint(i);
    cno = nei[0] + nei[1] * this->XD + nei[2] * this->SliceSize;

    if ((ptIds = this->HashTable[cno]) != nullptr)
    {
      nids = ptIds->GetNumberOfIds();
      for (j = 0; j < nids; j++)
      {
        pointsChecked++;
        ptId = ptIds->GetId(j);
        this->DataSet->GetPoint(ptId, pt);
        dist2 = vtkMath::Distance2BetweenPoints(x, pt);
        oct = GetOctant(x, pt);
        if (dist2 < maxDistance[oct])
        {
          res[oct][N - 1].Dist2 = dist2;
          res[oct][N - 1].PtId = ptId;
          std::sort(res[oct], res[oct] + N);
          maxDistance[oct] = res[oct][N - 1].Dist2;
        }
      }
    }
  }

  // Fill in the IdList
  for (j = 0; j < 8; j++)
  {
    for (i = 0; i < currentCount[j]; i++)
    {
      result->InsertNextId(res[j][i].PtId);
    }
    delete[] res[j];
  }
}

//------------------------------------------------------------------------------
void vtkPointLocator::FindClosestNPoints(int N, const double x[3], vtkIdList* result)
{
  this->BuildLocator();
  if (this->HashTable == nullptr)
  {
    return;
  }
  int i, j;
  double dist2;
  double pt[3];
  int level;
  vtkIdType ptId, nids, cno;
  vtkIdList* ptIds;
  int ijk[3], *nei;
  vtkNeighborPoints buckets;

  // clear out the result
  result->Reset();

  //  Find bucket point is in.
  this->GetBucketIndices(x, ijk);

  // there are two steps, first a simple expanding wave of buckets until
  // we have enough points. Then a refinement to make sure we have the
  // N closest points.
  level = 0;
  double maxDistance = 0.0;
  int currentCount = 0;
  IdTuple* res = new IdTuple[N];

  this->GetBucketNeighbors(&buckets, ijk, this->Divisions, level);
  while (buckets.GetNumberOfNeighbors() && currentCount < N)
  {
    for (i = 0; i < buckets.GetNumberOfNeighbors(); i++)
    {
      nei = buckets.GetPoint(i);
      cno = nei[0] + nei[1] * this->XD + nei[2] * this->SliceSize;

      if ((ptIds = this->HashTable[cno]) != nullptr)
      {
        nids = ptIds->GetNumberOfIds();
        for (j = 0; j < nids; j++)
        {
          ptId = ptIds->GetId(j);
          this->DataSet->GetPoint(ptId, pt);
          dist2 = vtkMath::Distance2BetweenPoints(x, pt);
          if (currentCount < N)
          {
            res[currentCount].Dist2 = dist2;
            res[currentCount].PtId = ptId;
            if (dist2 > maxDistance)
            {
              maxDistance = dist2;
            }
            currentCount++;
            if (currentCount == N)
            {
              std::sort(res, res + currentCount);
            }
          }
          else if (dist2 < maxDistance)
          {
            res[N - 1].Dist2 = dist2;
            res[N - 1].PtId = ptId;
            std::sort(res, res + N);
            maxDistance = res[N - 1].Dist2;
          }
        }
      }
    }
    level++;
    this->GetBucketNeighbors(&buckets, ijk, this->Divisions, level);
  }

  // do a sort
  std::sort(res, res + currentCount);

  // Now do the refinement
  this->GetOverlappingBuckets(&buckets, x, ijk, sqrt(maxDistance), level - 1);

  for (i = 0; i < buckets.GetNumberOfNeighbors(); i++)
  {
    nei = buckets.GetPoint(i);
    cno = nei[0] + nei[1] * this->XD + nei[2] * this->SliceSize;

    if ((ptIds = this->HashTable[cno]) != nullptr)
    {
      nids = ptIds->GetNumberOfIds();
      for (j = 0; j < nids; j++)
      {
        ptId = ptIds->GetId(j);
        this->DataSet->GetPoint(ptId, pt);
        dist2 = vtkMath::Distance2BetweenPoints(x, pt);
        if (dist2 < maxDistance)
        {
          res[N - 1].Dist2 = dist2;
          res[N - 1].PtId = ptId;
          std::sort(res, res + N);
          maxDistance = res[N - 1].Dist2;
        }
      }
    }
  }

  // Fill in the IdList
  result->SetNumberOfIds(currentCount);
  for (i = 0; i < currentCount; i++)
  {
    result->SetId(i, res[i].PtId);
  }

  delete[] res;
}

//------------------------------------------------------------------------------
void vtkPointLocator::FindPointsWithinRadius(double R, const double x[3], vtkIdList* result)
{
  this->BuildLocator();
  if (this->HashTable == nullptr)
  {
    return;
  }
  int i, j;
  double dist2;
  double pt[3];
  vtkIdType ptId, nids, cno;
  vtkIdList* ptIds;
  int ijk[3], *nei;
  double R2 = R * R;
  vtkNeighborPoints buckets;

  //  Find bucket point is in.
  this->GetBucketIndices(x, ijk);

  // get all buckets within a distance
  this->GetOverlappingBuckets(&buckets, x, ijk, R, 0);
  // add the original bucket
  buckets.InsertNextPoint(ijk);

  // clear out the result
  result->Reset();

  for (i = 0; i < buckets.GetNumberOfNeighbors(); i++)
  {
    nei = buckets.GetPoint(i);
    cno = nei[0] + nei[1] * this->XD + nei[2] * this->SliceSize;

    if ((ptIds = this->HashTable[cno]) != nullptr)
    {
      nids = ptIds->GetNumberOfIds();
      for (j = 0; j < nids; j++)
      {
        ptId = ptIds->GetId(j);
        this->DataSet->GetPoint(ptId, pt);
        dist2 = vtkMath::Distance2BetweenPoints(x, pt);
        if (dist2 <= R2)
        {
          result->InsertNextId(ptId);
        }
      }
    }
  }
}

//------------------------------------------------------------------------------
void vtkPointLocator::BuildLocator()
{
  // don't rebuild if build time is newer than modified and dataset modified time
  if (this->HashTable && this->BuildTime > this->MTime &&
    this->BuildTime > this->DataSet->GetMTime())
  {
    return;
  }
  // don't rebuild if UseExistingSearchStructure is ON and a search structure already exists
  if (this->HashTable && this->UseExistingSearchStructure)
  {
    this->BuildTime.Modified();
    vtkDebugMacro(<< "BuildLocator exited - UseExistingSearchStructure");
    return;
  }
  this->BuildLocatorInternal();
}

//------------------------------------------------------------------------------
void vtkPointLocator::ForceBuildLocator()
{
  this->BuildLocatorInternal();
}

//------------------------------------------------------------------------------
//  Method to form subdivision of space based on the points provided and
//  subject to the constraints of levels and NumberOfPointsPerBucket.
//  The result is directly addressable and of uniform subdivision.
void vtkPointLocator::BuildLocatorInternal()
{
  int ndivs[3];
  vtkIdType idx;
  vtkIdList* bucket;
  vtkIdType numPts;
  double x[3];
  typedef vtkIdList* vtkIdListPtr;

  vtkDebugMacro(<< "Hashing points...");
  this->Level = 1; // only single lowest level

  if (!this->DataSet || (numPts = this->DataSet->GetNumberOfPoints()) < 1)
  {
    vtkErrorMacro(<< "No points to subdivide");
    return;
  }

  //  Make sure the appropriate data is available
  this->FreeSearchStructure();

  //  Size the root bucket.  Initialize bucket data structure, compute
  //  level and divisions.
  const double* bounds = this->DataSet->GetBounds();
  vtkIdType numBuckets = static_cast<vtkIdType>(
    static_cast<double>(numPts) / static_cast<double>(this->NumberOfPointsPerBucket));

  vtkBoundingBox bbox(bounds);
  if (this->Automatic)
  {
    bbox.ComputeDivisions(numBuckets, this->Bounds, ndivs);
  }
  else
  {
    bbox.Inflate(); // make sure non-zero volume
    bbox.GetBounds(this->Bounds);
    for (int i = 0; i < 3; ++i)
    {
      ndivs[i] = (this->Divisions[i] < 1 ? 1 : this->Divisions[i]);
    }
  }

  this->Divisions[0] = ndivs[0];
  this->Divisions[1] = ndivs[1];
  this->Divisions[2] = ndivs[2];
  this->NumberOfBuckets = numBuckets = static_cast<vtkIdType>(ndivs[0]) *
    static_cast<vtkIdType>(ndivs[1]) * static_cast<vtkIdType>(ndivs[2]);

  //  Compute width of bucket in three directions
  //
  for (int i = 0; i < 3; ++i)
  {
    this->H[i] = (this->Bounds[2 * i + 1] - this->Bounds[2 * i]) / static_cast<double>(ndivs[i]);
  }

  // Allocate the bins/buckets and initialize
  this->HashTable = new vtkIdListPtr[numBuckets];
  // NOLINTNEXTLINE(bugprone-sizeof-expression)
  memset(this->HashTable, 0, numBuckets * sizeof(*this->HashTable));

  // Compute local variables (for performance reasons)
  // Setup internal data members for more efficient processing.
  this->ComputePerformanceFactors();

  //  Insert each point into the appropriate bucket.  Make sure point
  //  falls within bucket.
  //
  for (vtkIdType i = 0; i < numPts; ++i)
  {
    this->DataSet->GetPoint(i, x);
    idx = this->GetBucketIndex(x);
    bucket = this->HashTable[idx];
    if (!bucket)
    {
      bucket = vtkIdList::New();
      bucket->Allocate(this->NumberOfPointsPerBucket, this->NumberOfPointsPerBucket / 3);
      this->HashTable[idx] = bucket;
    }
    bucket->InsertNextId(i);
  }

  // Okay we're done update mtime
  this->BuildTime.Modified();
}

//------------------------------------------------------------------------------
//  Internal function to get bucket neighbors at specified level
//
void vtkPointLocator::GetBucketNeighbors(
  vtkNeighborPoints* buckets, const int ijk[3], const int ndivs[3], int level)
{
  int i, j, k, min, max, minLevel[3], maxLevel[3];
  int nei[3];
  //
  //  Initialize
  //
  buckets->Reset();
  //
  //  If at this bucket, just place into list
  //
  if (level == 0)
  {
    buckets->InsertNextPoint(ijk);
    return;
  }
  //
  //  Create permutations of the ijk indices that are at the level
  //  required. If these are legal buckets, add to list for searching.
  //
  for (i = 0; i < 3; i++)
  {
    min = ijk[i] - level;
    max = ijk[i] + level;
    minLevel[i] = (min > 0 ? min : 0);
    maxLevel[i] = (max < (ndivs[i] - 1) ? max : (ndivs[i] - 1));
  }

  for (i = minLevel[0]; i <= maxLevel[0]; i++)
  {
    for (j = minLevel[1]; j <= maxLevel[1]; j++)
    {
      for (k = minLevel[2]; k <= maxLevel[2]; k++)
      {
        if (i == (ijk[0] + level) || i == (ijk[0] - level) || j == (ijk[1] + level) ||
          j == (ijk[1] - level) || k == (ijk[2] + level) || k == (ijk[2] - level))
        {
          nei[0] = i;
          nei[1] = j;
          nei[2] = k;
          buckets->InsertNextPoint(nei);
        }
      }
    }
  }
}

//------------------------------------------------------------------------------
// Internal method to find those buckets that are within distance specified
// only those buckets outside of level radiuses of ijk are returned
void vtkPointLocator::GetOverlappingBuckets(
  vtkNeighborPoints* buckets, const double x[3], const int ijk[3], double dist, int level)
{
  int i, j, k, nei[3], minLevel[3], maxLevel[3];
  double xMin[3], xMax[3];

  // Initialize
  buckets->Reset();

  // Determine the range of indices in each direction
  // Determine the range of indices in each direction
  xMin[0] = x[0] - dist;
  xMin[1] = x[1] - dist;
  xMin[2] = x[2] - dist;
  xMax[0] = x[0] + dist;
  xMax[1] = x[1] + dist;
  xMax[2] = x[2] + dist;

  this->GetBucketIndices(xMin, minLevel);
  this->GetBucketIndices(xMax, maxLevel);

  for (i = minLevel[0]; i <= maxLevel[0]; i++)
  {
    for (j = minLevel[1]; j <= maxLevel[1]; j++)
    {
      for (k = minLevel[2]; k <= maxLevel[2]; k++)
      {
        if (i < (ijk[0] - level) || i > (ijk[0] + level) || j < (ijk[1] - level) ||
          j > (ijk[1] + level) || k < (ijk[2] - level) || k > (ijk[2] + level))
        {
          nei[0] = i;
          nei[1] = j;
          nei[2] = k;
          buckets->InsertNextPoint(nei);
        }
      }
    }
  }
}

//------------------------------------------------------------------------------
// Internal method to find those buckets that are within distance specified
// only those buckets outside of level radiuses of ijk are returned
void vtkPointLocator::GetOverlappingBuckets(vtkNeighborPoints* buckets, const double x[3],
  double dist, int prevMinLevel[3], int prevMaxLevel[3])
{
  int i, j, k, nei[3], minLevel[3], maxLevel[3];
  vtkIdType kFactor, jFactor;
  int jkSkipFlag, kSkipFlag;
  double xMin[3], xMax[3];

  // Initialize
  buckets->Reset();

  // Determine the range of indices in each direction
  // Determine the range of indices in each direction
  xMin[0] = x[0] - dist;
  xMin[1] = x[1] - dist;
  xMin[2] = x[2] - dist;
  xMax[0] = x[0] + dist;
  xMax[1] = x[1] + dist;
  xMax[2] = x[2] + dist;

  this->GetBucketIndices(xMin, minLevel);
  this->GetBucketIndices(xMax, maxLevel);

  if (minLevel[0] == prevMinLevel[0] && maxLevel[0] == prevMaxLevel[0] &&
    minLevel[1] == prevMinLevel[1] && maxLevel[1] == prevMaxLevel[1] &&
    minLevel[2] == prevMinLevel[2] && maxLevel[2] == prevMaxLevel[2])
  {
    return;
  }

  for (k = minLevel[2]; k <= maxLevel[2]; k++)
  {
    kFactor = k * this->SliceSize;
    if (k >= prevMinLevel[2] && k <= prevMaxLevel[2])
    {
      kSkipFlag = 1;
    }
    else
    {
      kSkipFlag = 0;
    }
    for (j = minLevel[1]; j <= maxLevel[1]; j++)
    {
      if (kSkipFlag && j >= prevMinLevel[1] && j <= prevMaxLevel[1])
      {
        jkSkipFlag = 1;
      }
      else
      {
        jkSkipFlag = 0;
      }
      jFactor = j * this->XD;
      for (i = minLevel[0]; i <= maxLevel[0]; i++)
      {
        if (jkSkipFlag && i == prevMinLevel[0])
        {
          i = prevMaxLevel[0];
          continue;
        }
        // if this bucket has any cells, add it to the list
        if (this->HashTable[i + jFactor + kFactor])
        {
          nei[0] = i;
          nei[1] = j;
          nei[2] = k;
          buckets->InsertNextPoint(nei);
        }
      }
    }
  }

  prevMinLevel[0] = minLevel[0];
  prevMinLevel[1] = minLevel[1];
  prevMinLevel[2] = minLevel[2];
  prevMaxLevel[0] = maxLevel[0];
  prevMaxLevel[1] = maxLevel[1];
  prevMaxLevel[2] = maxLevel[2];
}

//------------------------------------------------------------------------------
// Initialize the point insertion process. The newPts is an object representing
// point coordinates into which incremental insertion methods place their
// data. Bounds are the box that the points lie in.
int vtkPointLocator::InitPointInsertion(vtkPoints* newPts, const double bounds[6])
{
  return this->InitPointInsertion(newPts, bounds, 0);
}

//------------------------------------------------------------------------------
// Initialize the point insertion process. The newPts is an object representing
// point coordinates into which incremental insertion methods place their
// data. Bounds are the box that the points lie in.
int vtkPointLocator::InitPointInsertion(
  vtkPoints* newPts, const double bounds[6], vtkIdType estNumPts)
{
  int i;
  int maxDivs;
  typedef vtkIdList* vtkIdListPtr;
  double hmin;
  int ndivs[3];

  this->InsertionPointId = 0;
  if (this->HashTable)
  {
    this->FreeSearchStructure();
  }
  if (newPts == nullptr)
  {
    vtkErrorMacro(<< "Must define points for point insertion");
    return 0;
  }
  if (this->Points != nullptr)
  {
    this->Points->UnRegister(this);
  }
  this->Points = newPts;
  this->Points->Register(this);

  // Configure the locator
  vtkIdType numBuckets;
  vtkBoundingBox bbox(bounds);
  if (this->Automatic && (estNumPts > 0))
  {
    numBuckets = static_cast<vtkIdType>(
      static_cast<double>(estNumPts) / static_cast<double>(this->NumberOfPointsPerBucket));
    bbox.ComputeDivisions(numBuckets, this->Bounds, ndivs);
  }
  else
  {
    bbox.Inflate(); // make sure non-zero volume
    bbox.GetBounds(this->Bounds);
    for (i = 0; i < 3; i++)
    {
      ndivs[i] = (this->Divisions[i] < 1 ? 1 : this->Divisions[i]);
    }
  }

  this->Divisions[0] = ndivs[0];
  this->Divisions[1] = ndivs[1];
  this->Divisions[2] = ndivs[2];
  this->NumberOfBuckets = numBuckets = static_cast<vtkIdType>(ndivs[0]) *
    static_cast<vtkIdType>(ndivs[1]) * static_cast<vtkIdType>(ndivs[2]);

  // Initialize bins/buckets
  this->HashTable = new vtkIdListPtr[this->NumberOfBuckets];
  // NOLINTNEXTLINE(bugprone-sizeof-expression)
  memset(this->HashTable, 0, this->NumberOfBuckets * sizeof(*this->HashTable));

  //  Compute width of bucket in three directions
  //
  for (i = 0; i < 3; i++)
  {
    this->H[i] = (this->Bounds[2 * i + 1] - this->Bounds[2 * i]) / ndivs[i];
  }

  this->InsertionTol2 = this->Tolerance * this->Tolerance;

  for (maxDivs = 0, hmin = VTK_DOUBLE_MAX, i = 0; i < 3; i++)
  {
    hmin = (this->H[i] < hmin ? this->H[i] : hmin);
    maxDivs = (maxDivs > this->Divisions[i] ? maxDivs : this->Divisions[i]);
  }
  this->InsertionLevel = ceil(static_cast<double>(this->Tolerance) / hmin);
  this->InsertionLevel = (this->InsertionLevel > maxDivs ? maxDivs : this->InsertionLevel);

  // Update internal performance variables
  this->ComputePerformanceFactors();

  return 1;
}

//------------------------------------------------------------------------------
// Incrementally insert a point into search structure. The method returns
// the insertion location (i.e., point id). You should use the method
// IsInsertedPoint() to see whether this point has already been
// inserted (that is, if you desire to prevent duplicate points).
// Before using this method you must make sure that newPts have been
// supplied, the bounds has been set properly, and that divs are
// properly set. (See InitPointInsertion().)
vtkIdType vtkPointLocator::InsertNextPoint(const double x[3])
{
  vtkIdType idx;
  vtkIdList* bucket;

  idx = this->GetBucketIndex(x);

  if (!(bucket = this->HashTable[idx]))
  {
    bucket = vtkIdList::New();
    bucket->Allocate(this->NumberOfPointsPerBucket / 2, this->NumberOfPointsPerBucket / 3);
    this->HashTable[idx] = bucket;
  }

  bucket->InsertNextId(this->InsertionPointId);
  this->Points->InsertPoint(this->InsertionPointId, x);
  return this->InsertionPointId++;
}

//------------------------------------------------------------------------------
// Incrementally insert a point into search structure with a particular
// index value. You should use the method IsInsertedPoint() to see whether
// this point has already been inserted (that is, if you desire to prevent
// duplicate points). Before using this method you must make sure that
// newPts have been supplied, the bounds has been set properly, and that
// divs are properly set. (See InitPointInsertion().)
void vtkPointLocator::InsertPoint(vtkIdType ptId, const double x[3])
{
  vtkIdType idx;
  vtkIdList* bucket;

  idx = this->GetBucketIndex(x);

  if (!(bucket = this->HashTable[idx]))
  {
    bucket = vtkIdList::New();
    bucket->Allocate(this->NumberOfPointsPerBucket, this->NumberOfPointsPerBucket / 3);
    this->HashTable[idx] = bucket;
  }

  bucket->InsertNextId(ptId);
  this->Points->InsertPoint(ptId, x);
}

//------------------------------------------------------------------------------
// Determine whether point given by x[3] has been inserted into points list.
// Return id of previously inserted point if this is true, otherwise return
// -1.
vtkIdType vtkPointLocator::IsInsertedPoint(const double x[3])
{
  int i, j, ijk[3];
  vtkNeighborPoints buckets;
  vtkIdType nids;

  //  Locate bucket that point is in.
  //
  this->GetBucketIndices(x, ijk);

  // Check the list of points in that bucket for merging.  Also need to
  // search all neighboring buckets within the tolerance.  The number
  // and level of neighbors to search depends upon the tolerance and
  // the bucket width.
  //
  int *nei, lvtk;
  vtkIdType ptId, cno;
  vtkIdList* ptIds;
  double pt[3];

  for (lvtk = 0; lvtk <= this->InsertionLevel; lvtk++)
  {
    this->GetBucketNeighbors(&buckets, ijk, this->Divisions, lvtk);

    for (i = 0; i < buckets.GetNumberOfNeighbors(); i++)
    {
      nei = buckets.GetPoint(i);
      cno = nei[0] + nei[1] * this->XD + nei[2] * this->SliceSize;

      if ((ptIds = this->HashTable[cno]) != nullptr)
      {
        nids = ptIds->GetNumberOfIds();
        for (j = 0; j < nids; j++)
        {
          ptId = ptIds->GetId(j);
          this->Points->GetPoint(ptId, pt);

          if (vtkMath::Distance2BetweenPoints(x, pt) <= this->InsertionTol2)
          {
            return ptId;
          }
        }
      } // if points in bucket
    }   // for each neighbor
  }     // for neighbors at this level

  return -1;
}

//------------------------------------------------------------------------------
int vtkPointLocator::InsertUniquePoint(const double x[3], vtkIdType& id)
{
  vtkIdType ptId;

  ptId = this->IsInsertedPoint(x);

  if (ptId > -1)
  {
    id = ptId;
    return 0;
  }
  else
  {
    id = this->InsertNextPoint(x);
    return 1;
  }
}

//------------------------------------------------------------------------------
// Given a position x, return the id of the point closest to it. This method
// is used when performing incremental point insertion.
vtkIdType vtkPointLocator::FindClosestInsertedPoint(const double x[3])
{
  int i;
  double minDist2, dist2;
  double pt[3];
  int level;
  vtkIdType closest, j;
  vtkIdType ptId, nids, cno;
  vtkIdList* ptIds;
  int ijk[3], *nei;
  int MULTIPLES;
  double diff;
  vtkNeighborPoints buckets;

  //
  //  Make sure candidate point is in bounds.  If not, it is outside.
  //
  for (i = 0; i < 3; i++)
  {
    if (x[i] < this->Bounds[2 * i] || x[i] > this->Bounds[2 * i + 1])
    {
      return -1;
    }
  }

  //  Find bucket point is in.
  //
  this->GetBucketIndices(x, ijk);

  //  Need to search this bucket for closest point.  If there are no
  //  points in this bucket, search 1st level neighbors, and so on,
  //  until closest point found.
  //
  for (closest = -1, minDist2 = VTK_DOUBLE_MAX, level = 0; (closest == -1) &&
       (level < this->Divisions[0] || level < this->Divisions[1] || level < this->Divisions[2]);
       level++)
  {
    this->GetBucketNeighbors(&buckets, ijk, this->Divisions, level);

    for (i = 0; i < buckets.GetNumberOfNeighbors(); i++)
    {
      nei = buckets.GetPoint(i);
      cno = nei[0] + nei[1] * this->XD + nei[2] * this->SliceSize;

      if ((ptIds = this->HashTable[cno]) != nullptr)
      {
        nids = ptIds->GetNumberOfIds();
        for (j = 0; j < nids; j++)
        {
          ptId = ptIds->GetId(j);
          this->Points->GetPoint(ptId, pt);
          if ((dist2 = vtkMath::Distance2BetweenPoints(x, pt)) < minDist2)
          {
            closest = ptId;
            minDist2 = dist2;
          }
        }
      }
    }
  }
  //
  //  Because of the relative location of the points in the spatial_hash, this
  //  may not be the closest point.  Have to search those bucket
  //  neighbors (one level further out) that might also contain point.
  //
  this->GetBucketNeighbors(&buckets, ijk, this->Divisions, level);
  //
  //  Don't want to search all the neighbors, only those that could
  //  possibly have points closer than the current closest.
  //
  for (i = 0; i < buckets.GetNumberOfNeighbors(); i++)
  {
    nei = buckets.GetPoint(i);

    for (dist2 = 0, j = 0; j < 3; j++)
    {
      if (ijk[j] != nei[j])
      {
        MULTIPLES = (ijk[j] > nei[j] ? (nei[j] + 1) : nei[j]);
        diff = (this->Bounds[2 * j] + MULTIPLES * this->H[j]) - x[j];
        dist2 += diff * diff;
      }
    }

    if (dist2 < minDist2)
    {
      cno = nei[0] + nei[1] * this->XD + nei[2] * this->SliceSize;

      if ((ptIds = this->HashTable[cno]))
      {
        nids = ptIds->GetNumberOfIds();
        for (j = 0; j < nids; j++)
        {
          ptId = ptIds->GetId(j);
          this->Points->GetPoint(ptId, pt);
          if ((dist2 = vtkMath::Distance2BetweenPoints(x, pt)) < minDist2)
          {
            closest = ptId;
            minDist2 = dist2;
          }
        }
      }
    }
  }

  return closest;
}

//------------------------------------------------------------------------------
// Return the list of points in the bucket containing x.
vtkIdList* vtkPointLocator::GetPointsInBucket(const double x[3], int ijk[3])
{
  int i;

  //  Make sure candidate point is in bounds.  If not, it is outside.
  //
  for (i = 0; i < 3; i++)
  {
    if (x[i] < this->Bounds[2 * i] || x[i] > this->Bounds[2 * i + 1])
    {
      return nullptr;
    }
  }

  this->GetBucketIndices(x, ijk);

  // Get the id list, if any
  if (this->HashTable)
  {
    vtkIdType idx = this->GetBucketIndex(x);
    return this->HashTable[idx];
  }

  return nullptr;
}

//------------------------------------------------------------------------------
// Build polygonal representation of locator. Create faces that separate
// inside/outside buckets, or separate inside/boundary of locator.
void vtkPointLocator::GenerateRepresentation(int vtkNotUsed(level), vtkPolyData* pd)
{
  vtkPoints* pts;
  vtkCellArray* polys;
  int ii, i, j, k, idx, offset[3], minusOffset[3], inside, sliceSize;

  if (this->HashTable == nullptr)
  {
    vtkErrorMacro(<< "Can't build representation...no data!");
    return;
  }

  pts = vtkPoints::New();
  pts->Allocate(5000);
  polys = vtkCellArray::New();
  polys->AllocateEstimate(2048, 3);

  // loop over all buckets, creating appropriate faces
  sliceSize = this->Divisions[0] * this->Divisions[1];
  for (k = 0; k < this->Divisions[2]; k++)
  {
    offset[2] = k * sliceSize;
    minusOffset[2] = (k - 1) * sliceSize;
    for (j = 0; j < this->Divisions[1]; j++)
    {
      offset[1] = j * this->Divisions[0];
      minusOffset[1] = (j - 1) * this->Divisions[0];
      for (i = 0; i < this->Divisions[0]; i++)
      {
        offset[0] = i;
        minusOffset[0] = i - 1;
        idx = offset[0] + offset[1] + offset[2];
        if (this->HashTable[idx] == nullptr)
        {
          inside = 0;
        }
        else
        {
          inside = 1;
        }

        // check "negative" neighbors
        for (ii = 0; ii < 3; ii++)
        {
          if (minusOffset[ii] < 0)
          {
            if (inside)
            {
              this->GenerateFace(ii, i, j, k, pts, polys);
            }
          }
          else
          {
            if (ii == 0)
            {
              idx = minusOffset[0] + offset[1] + offset[2];
            }
            else if (ii == 1)
            {
              idx = offset[0] + minusOffset[1] + offset[2];
            }
            else
            {
              idx = offset[0] + offset[1] + minusOffset[2];
            }

            if ((this->HashTable[idx] == nullptr && inside) ||
              (this->HashTable[idx] != nullptr && !inside))
            {
              this->GenerateFace(ii, i, j, k, pts, polys);
            }
          }
          // those buckets on "positive" boundaries can generate faces specially
          if ((i + 1) >= this->Divisions[0] && inside)
          {
            this->GenerateFace(0, i + 1, j, k, pts, polys);
          }
          if ((j + 1) >= this->Divisions[1] && inside)
          {
            this->GenerateFace(1, i, j + 1, k, pts, polys);
          }
          if ((k + 1) >= this->Divisions[2] && inside)
          {
            this->GenerateFace(2, i, j, k + 1, pts, polys);
          }

        } // over negative faces
      }   // over i divisions
    }     // over j divisions
  }       // over k divisions

  pd->SetPoints(pts);
  pts->Delete();
  pd->SetPolys(polys);
  polys->Delete();
  pd->Squeeze();
}

//------------------------------------------------------------------------------
void vtkPointLocator::GenerateFace(
  int face, int i, int j, int k, vtkPoints* pts, vtkCellArray* polys)
{
  vtkIdType ids[4];
  double origin[3], x[3];

  // define first corner
  origin[0] = this->Bounds[0] + i * this->HX;
  origin[1] = this->Bounds[2] + j * this->HY;
  origin[2] = this->Bounds[4] + k * this->HZ;
  ids[0] = pts->InsertNextPoint(origin);

  if (face == 0) // x face
  {
    x[0] = origin[0];
    x[1] = origin[1] + this->HY;
    x[2] = origin[2];
    ids[1] = pts->InsertNextPoint(x);

    x[0] = origin[0];
    x[1] = origin[1] + this->HY;
    x[2] = origin[2] + this->HZ;
    ids[2] = pts->InsertNextPoint(x);

    x[0] = origin[0];
    x[1] = origin[1];
    x[2] = origin[2] + this->HZ;
    ids[3] = pts->InsertNextPoint(x);
  }

  else if (face == 1) // y face
  {
    x[0] = origin[0] + this->HX;
    x[1] = origin[1];
    x[2] = origin[2];
    ids[1] = pts->InsertNextPoint(x);

    x[0] = origin[0] + this->HX;
    x[1] = origin[1];
    x[2] = origin[2] + this->HZ;
    ids[2] = pts->InsertNextPoint(x);

    x[0] = origin[0];
    x[1] = origin[1];
    x[2] = origin[2] + this->HZ;
    ids[3] = pts->InsertNextPoint(x);
  }

  else // z face
  {
    x[0] = origin[0] + this->HX;
    x[1] = origin[1];
    x[2] = origin[2];
    ids[1] = pts->InsertNextPoint(x);

    x[0] = origin[0] + this->HX;
    x[1] = origin[1] + this->HY;
    x[2] = origin[2];
    ids[2] = pts->InsertNextPoint(x);

    x[0] = origin[0];
    x[1] = origin[1] + this->HY;
    x[2] = origin[2];
    ids[3] = pts->InsertNextPoint(x);
  }

  polys->InsertNextCell(4, ids);
}

//------------------------------------------------------------------------------
// Calculate the distance between the point x to the bucket "nei".
//
// WARNING!!!!! Be very careful altering this routine.  Simple changes to this
// routine can make is 25% slower!!!!
//
double vtkPointLocator::Distance2ToBucket(const double x[3], const int nei[3])
{
  double bounds[6];

  bounds[0] = nei[0] * this->HX + this->BX;
  bounds[1] = (nei[0] + 1) * this->HX + this->BX;
  bounds[2] = nei[1] * this->HY + this->BY;
  bounds[3] = (nei[1] + 1) * this->HY + this->BY;
  bounds[4] = nei[2] * this->HZ + this->BZ;
  bounds[5] = (nei[2] + 1) * this->HZ + this->BZ;

  return this->Distance2ToBounds(x, bounds);
}

//------------------------------------------------------------------------------
// Calculate the distance between the point x and the specified bounds
//
// WARNING!!!!! Be very careful altering this routine.  Simple changes to this
// routine can make is 25% slower!!!!
double vtkPointLocator::Distance2ToBounds(const double x[3], const double bounds[6])
{
  double distance;
  double deltas[3];

  // Are we within the bounds?
  if (x[0] >= bounds[0] && x[0] <= bounds[1] && x[1] >= bounds[2] && x[1] <= bounds[3] &&
    x[2] >= bounds[4] && x[2] <= bounds[5])
  {
    return 0.0;
  }

  deltas[0] = deltas[1] = deltas[2] = 0.0;

  // dx
  //
  if (x[0] < bounds[0])
  {
    deltas[0] = bounds[0] - x[0];
  }
  else if (x[0] > bounds[1])
  {
    deltas[0] = x[0] - bounds[1];
  }

  // dy
  //
  if (x[1] < bounds[2])
  {
    deltas[1] = bounds[2] - x[1];
  }
  else if (x[1] > bounds[3])
  {
    deltas[1] = x[1] - bounds[3];
  }

  // dz
  //
  if (x[2] < bounds[4])
  {
    deltas[2] = bounds[4] - x[2];
  }
  else if (x[2] > bounds[5])
  {
    deltas[2] = x[2] - bounds[5];
  }

  distance = vtkMath::Dot(deltas, deltas);
  return distance;
}

//------------------------------------------------------------------------------
void vtkPointLocator::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os, indent);

  os << indent << "Number of Points Per Bucket: " << this->NumberOfPointsPerBucket << "\n";
  os << indent << "Divisions: (" << this->Divisions[0] << ", " << this->Divisions[1] << ", "
     << this->Divisions[2] << ")\n";
  if (this->Points)
  {
    os << indent << "Points:\n";
    this->Points->PrintSelf(os, indent.GetNextIndent());
  }
  else
  {
    os << indent << "Points: (none)\n";
  }
}