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

  Program:   Visualization Toolkit
  Module:    $RCSfile: vtkPointLocator.cxx,v $

  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 "vtkCellArray.h"
#include "vtkIdList.h"
#include "vtkIntArray.h"
#include "vtkMath.h"
#include "vtkObjectFactory.h"
#include "vtkPolyData.h"

vtkCxxRevisionMacro(vtkPointLocator, "$Revision: 1.1 $");
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[0]);
      this->MaxSize = VTK_INITIAL_SIZE;
    }
  ~vtkNeighborPoints() 
    { 
      this->Count = 0; 
      if ( this->P != &(this->InitialBuffer[0]) )
        {
        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]) : 0);
    }

  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[0]) )
          {
          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
// 25 points per bucket.
vtkPointLocator::vtkPointLocator()
{

  this->Points = NULL;
  this->Divisions[0] = this->Divisions[1] = this->Divisions[2] = 50;
  this->NumberOfPointsPerBucket = 3;
  this->HashTable = NULL;
  this->NumberOfBuckets = 0;
  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 = NULL;
    }
  this->FreeSearchStructure();
}

void vtkPointLocator::Initialize()
{
  if ( this->Points )
    {
    this->Points->UnRegister(this);
    this->Points = NULL;
    }
  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 = NULL;
    }
}

// Given a position x-y-z, return the id of the point closest to it.
vtkIdType vtkPointLocator::FindClosestPoint(double x, double y, double z)
{
  double xyz[3];

  xyz[0] = x; xyz[1] = y; xyz[2] = z;
  return this->FindClosestPoint(xyz);
}

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

  this->BuildLocator(); // will subdivide if modified; otherwise returns

  //
  //  Find bucket point is in.  
  //
  for (j=0; j<3; j++) 
    {
    ijk[j] = (int)(((x[j] - this->Bounds[2*j]) / 
        (this->Bounds[2*j+1] - this->Bounds[2*j])) * this->Divisions[j]);

    if (ijk[j] < 0)
      {
      ijk[j] = 0;
      }
    else if (ijk[j] >= this->Divisions[j])
      { 
      ijk[j] = this->Divisions[j] - 1;
      }
    }
  //
  //  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->Divisions[0] + 
            nei[2]*this->Divisions[0]*this->Divisions[1];

      if ( (ptIds = this->HashTable[cno]) != NULL )
        {
        for (j=0; j < ptIds->GetNumberOfIds(); j++) 
          {
          ptId = ptIds->GetId(j);
          pt = this->DataSet->GetPoint(ptId);
          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 ( dist2 > 0.0 )
    {
    this->GetOverlappingBuckets (&buckets, x, ijk, sqrt(dist2),0);
    for (i=0; i<buckets.GetNumberOfNeighbors(); i++) 
      {
      nei = buckets.GetPoint(i);
      cno = nei[0] + nei[1]*this->Divisions[0] + 
            nei[2]*this->Divisions[0]*this->Divisions[1];

      if ( (ptIds = this->HashTable[cno]) != NULL )
        {
        for (j=0; j < ptIds->GetNumberOfIds(); j++) 
          {
          ptId = ptIds->GetId(j);
          pt = this->DataSet->GetPoint(ptId);
          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 FindClosestPointWithinRadius(radius, x, this->DataSet->GetLength(),
                                      dist2);
  }


vtkIdType vtkPointLocator::FindClosestPointWithinRadius(double radius,
                                                        const double x[3], 
                                                        double inputDataLength,
                                                        double& dist2)
{
  int i, j;
  double *pt;
  vtkIdType ptId, 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;

  this->BuildLocator(); // will subdivide if modified; otherwise returns

  dist2 = -1.0;
  radius2 = radius*radius;
  minDist2 = 1.01*radius2;   // something slightly bigger....
  
  vtkDataArray *pointData = ((vtkPointSet *)this->DataSet)->GetPoints()->GetData();
  int flag = 1;

  //
  //  Find bucket point is in.  
  //
  for (j=0; j<3; j++) 
    {
    ijk[j] = (int)(((x[j] - this->Bounds[2*j]) / 
        (this->Bounds[2*j+1] - this->Bounds[2*j])) * this->Divisions[j]);

    if (ijk[j] < 0)
      {
      ijk[j] = 0;
      }
    else if (ijk[j] >= this->Divisions[j])
      { 
      ijk[j] = this->Divisions[j] - 1;
      }
    }

  // Start by searching the bucket that the point is in.
  //
  if ( (ptIds = this->HashTable[ijk[0] + ijk[1]*this->Divisions[0] + 
    ijk[2]*this->Divisions[0]*this->Divisions[1]]) != NULL )
    {
    for (j=0; j < ptIds->GetNumberOfIds(); j++) 
      {
      ptId = ptIds->GetId(j);
      if (flag)
        {
        pt = pointData->GetTuple(ptId);
        }
      else
        {
        pt = this->DataSet->GetPoint(ptId);
        }
      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 (dist2 < radius2 && dist2 >= 0.0)
    {
    refinedRadius = sqrt(dist2);
    refinedRadius2 = dist2;
    }
  else
    {
    refinedRadius = radius;
    refinedRadius2 = radius2;
    }

  if (inputDataLength)
    {
    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] = (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.
  int numberOfBucketsPerPlane;
  numberOfBucketsPerPlane = this->Divisions[0]*this->Divisions[1];
  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->Divisions[0] + 
          nei[2]*numberOfBucketsPerPlane];

        for (j=0; j < ptIds->GetNumberOfIds(); j++) 
          {
          ptId = ptIds->GetId(j);
          if (flag)
            {
            pt = pointData->GetTuple(ptId);
            }
          else
            {
            pt = this->DataSet->GetPoint(ptId);
            }
          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 = (int)((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;
  }



struct idsort
{
  vtkIdType id;
  double dist;
};

#ifdef _WIN32_WCE
static int __cdecl vtkidsortcompare(const void *arg1, const void *arg2)
#else
extern "C" 
{  
  int vtkidsortcompare(const void *arg1, const void *arg2)
#endif
{
  idsort *v1 = (idsort *)arg1;
  idsort *v2 = (idsort *)arg2;
  if (v1->dist < v2->dist)
    {
    return -1;
    }
  if (v1->dist > v2->dist)
    {
    return 1;
    }
  return 0;  
}
#ifndef _WIN32_WCE
} // close extern "C"
#endif

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 GetOctent(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)
{
  int i, j;
  double dist2;
  double *pt;
  int level;
  vtkIdType ptId, cno;
  vtkIdList *ptIds;
  int ijk[3], *nei;
  int oct;
  int pointsChecked = 0;
  vtkNeighborPoints buckets;
  
  // clear out the result
  result->Reset();
    
  this->BuildLocator(); // will subdivide if modified; otherwise returns
  //
  //  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.  
  //
  for (j=0; j<3; j++) 
    {
    ijk[j] = (int)(((x[j] - this->Bounds[2*j]) / 
          (this->Bounds[2*j+1] - this->Bounds[2*j])) * this->Divisions[j]);
    if (ijk[j] >= this->Divisions[j])
      {
      ijk[j] = this->Divisions[j] - 1;
      }
    }

  // 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;
  
  idsort *res[8];
  for (i = 0; i < 8; i++) 
    {
    res[i] = new idsort [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->Divisions[0] + 
            nei[2]*this->Divisions[0]*this->Divisions[1];

      if ( (ptIds = this->HashTable[cno]) != NULL )
        {
        for (j=0; j < ptIds->GetNumberOfIds(); j++) 
          {
          pointsChecked++;
          ptId = ptIds->GetId(j);
          pt = this->DataSet->GetPoint(ptId);
          dist2 = vtkMath::Distance2BetweenPoints(x,pt);
          oct = GetOctent(x,pt);
          if (currentCount[oct] < N)
            {
            res[oct][currentCount[oct]].dist = dist2;
            res[oct][currentCount[oct]].id = ptId;
            if (dist2 > maxDistance[oct])
              {
              maxDistance[oct] = dist2;
              }
            currentCount[oct] = currentCount[oct] + 1;
            // compute new minCurrentCount
            minCurrentCount = GetMin(currentCount);         
            if (currentCount[oct] == N)
              {
              qsort(res[oct], currentCount[oct], sizeof(idsort),vtkidsortcompare);
              }
            }
          else if (dist2 < maxDistance[oct])
            {
            res[oct][N-1].dist = dist2;
            res[oct][N-1].id = ptId;
            qsort(res[oct], N, sizeof(idsort), vtkidsortcompare);
            maxDistance[oct] = res[oct][N-1].dist;
            }
          }
        }
      }
    level++;
    this->GetBucketNeighbors (&buckets, ijk, this->Divisions, level);
    }

  // do a sort
  for (i = 0; i < 8; i++)
    {
    qsort(res[i], currentCount[i], sizeof(idsort), vtkidsortcompare);
    }
  
  // 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->Divisions[0] + 
      nei[2]*this->Divisions[0]*this->Divisions[1];
    
    if ( (ptIds = this->HashTable[cno]) != NULL )
      {
      for (j=0; j < ptIds->GetNumberOfIds(); j++) 
        {
        pointsChecked++;
        ptId = ptIds->GetId(j);
        pt = this->DataSet->GetPoint(ptId);
        dist2 = vtkMath::Distance2BetweenPoints(x,pt);
        oct = GetOctent(x,pt);
        if (dist2 < maxDistance[oct])
          {
          res[oct][N-1].dist = dist2;
          res[oct][N-1].id = ptId;
          qsort(res[oct], N, sizeof(idsort), vtkidsortcompare);
          maxDistance[oct] = res[oct][N-1].dist;
          }
        }
      }
    }
    
  // Fill in the IdList
  for (j = 0; j < 8; j++)
    {
    for (i = 0; i < currentCount[j]; i++)
      {
      result->InsertNextId(res[j][i].id);
      }
    delete [] res[j];
    }
}

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

void vtkPointLocator::FindClosestNPoints(int N, const double x[3],
                                         vtkIdList *result)
{
  int i, j;
  double dist2;
  double *pt;
  int level;
  vtkIdType ptId, cno;
  vtkIdList *ptIds;
  int ijk[3], *nei;
  vtkNeighborPoints buckets;
  
  // clear out the result
  result->Reset();
    
  this->BuildLocator(); // will subdivide if modified; otherwise returns

  //
  //  Find bucket point is in.  
  //
  for (j=0; j<3; j++) 
    {
    ijk[j] = (int)(((x[j] - this->Bounds[2*j]) / 
        (this->Bounds[2*j+1] - this->Bounds[2*j])) * this->Divisions[j]);

    if (ijk[j] < 0)
      {
      ijk[j] = 0;
      }
    else if (ijk[j] >= this->Divisions[j])
      { 
      ijk[j] = this->Divisions[j] - 1;
      }
    }

  // 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;
  idsort *res = new idsort [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->Divisions[0] + 
            nei[2]*this->Divisions[0]*this->Divisions[1];

      if ( (ptIds = this->HashTable[cno]) != NULL )
        {
        for (j=0; j < ptIds->GetNumberOfIds(); j++) 
          {
          ptId = ptIds->GetId(j);
          pt = this->DataSet->GetPoint(ptId);
          dist2 = vtkMath::Distance2BetweenPoints(x,pt);
          if (currentCount < N)
            {
            res[currentCount].dist = dist2;
            res[currentCount].id = ptId;
            if (dist2 > maxDistance)
              {
              maxDistance = dist2;
              }
            currentCount++;
            if (currentCount == N)
              {
              qsort(res, currentCount, sizeof(idsort), vtkidsortcompare);
              }
            }
          else if (dist2 < maxDistance)
            {
            res[N-1].dist = dist2;
            res[N-1].id = ptId;
            qsort(res, N, sizeof(idsort), vtkidsortcompare);
            maxDistance = res[N-1].dist;
            }
          }
        }
      }
    level++;
    this->GetBucketNeighbors (&buckets, ijk, this->Divisions, level);
    }

  // do a sort
  qsort(res, currentCount, sizeof(idsort), vtkidsortcompare);

  // 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->Divisions[0] + 
      nei[2]*this->Divisions[0]*this->Divisions[1];
    
    if ( (ptIds = this->HashTable[cno]) != NULL )
      {
      for (j=0; j < ptIds->GetNumberOfIds(); j++) 
        {
        ptId = ptIds->GetId(j);
        pt = this->DataSet->GetPoint(ptId);
        dist2 = vtkMath::Distance2BetweenPoints(x,pt);
        if (dist2 < maxDistance)
          {
          res[N-1].dist = dist2;
          res[N-1].id = ptId;
          qsort(res, N, sizeof(idsort), vtkidsortcompare);
          maxDistance = res[N-1].dist;
          }
        }
      }
    }
    
  // Fill in the IdList
  result->SetNumberOfIds(currentCount);
  for (i = 0; i < currentCount; i++)
    {
    result->SetId(i,res[i].id);
    }

  delete [] res;
}


void vtkPointLocator::FindPointsWithinRadius(double R, double x,
                                             double y, double z,
                                             vtkIdList *result)
{
  double p[3];
  p[0] = x;
  p[1] = y;
  p[2] = z;
  this->FindPointsWithinRadius(R,p,result);
}



void vtkPointLocator::FindPointsWithinRadius(double R, const double x[3],
                                             vtkIdList *result)
{
  int i, j;
  double dist2;
  double *pt;
  vtkIdType ptId, cno;
  vtkIdList *ptIds;
  int ijk[3], *nei;
  double R2 = R*R;
  vtkNeighborPoints buckets;
  
  this->BuildLocator(); // will subdivide if modified; otherwise returns
  //
  //  Find bucket point is in.  
  //
  for (j=0; j<3; j++) 
    {
    ijk[j] = (int)(((x[j] - this->Bounds[2*j]) / 
        (this->Bounds[2*j+1] - this->Bounds[2*j])) * this->Divisions[j]);

    if (ijk[j] < 0)
      {
      ijk[j] = 0;
      }
    else if (ijk[j] >= this->Divisions[j])
      { 
      ijk[j] = this->Divisions[j] - 1;
      }
    }

  // 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->Divisions[0] + 
      nei[2]*this->Divisions[0]*this->Divisions[1];
    
    if ( (ptIds = this->HashTable[cno]) != NULL )
      {
      for (j=0; j < ptIds->GetNumberOfIds(); j++) 
        {
        ptId = ptIds->GetId(j);
        pt = this->DataSet->GetPoint(ptId);
        dist2 = vtkMath::Distance2BetweenPoints(x,pt);
        if (dist2 <= R2)
          {
          result->InsertNextId(ptId);
          }
        }
      }
    }

}

//
//  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::BuildLocator()
{
  double *bounds;
  vtkIdType numBuckets;
  double level;
  int ndivs[3], product;
  int i, j, ijk[3];
  vtkIdType idx;
  vtkIdList *bucket;
  vtkIdType numPts;
  double *x;
  typedef vtkIdList *vtkIdListPtr;

  if ( (this->HashTable != NULL) && (this->BuildTime > this->MTime)
       && (this->BuildTime > this->DataSet->GetMTime()) )
    {
    return;
    }

  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
  //
  if ( this->HashTable )
    {
    this->FreeSearchStructure();
    }
  //
  //  Size the root bucket.  Initialize bucket data structure, compute 
  //  level and divisions.
  //
  bounds = this->DataSet->GetBounds();
  for (i=0; i<3; i++)
    {
    this->Bounds[2*i] = bounds[2*i];
    this->Bounds[2*i+1] = bounds[2*i+1];
    if ( this->Bounds[2*i+1] <= this->Bounds[2*i] ) //prevent zero width
      {
      this->Bounds[2*i+1] = this->Bounds[2*i] + 1.0;
      }
    }

  if ( this->Automatic ) 
    {
    level = (double) numPts / this->NumberOfPointsPerBucket;
    level = ceil( pow((double)level,(double)0.33333333) );
    for (i=0; i<3; i++)
      {
      ndivs[i] = (int) level;
      }
    } 
  else 
    {
    for (i=0; i<3; i++)
      {
      ndivs[i] = (int) this->Divisions[i];
      }
    }

  for (i=0; i<3; i++) 
    {
    ndivs[i] = (ndivs[i] > 0 ? ndivs[i] : 1);
    this->Divisions[i] = ndivs[i];
    }

  this->NumberOfBuckets = numBuckets = ndivs[0]*ndivs[1]*ndivs[2];
  this->HashTable = new vtkIdListPtr[numBuckets];
  memset (this->HashTable, 0, numBuckets*sizeof(vtkIdListPtr));
  //
  //  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] ;
    }
  //
  //  Insert each point into the appropriate bucket.  Make sure point
  //  falls within bucket.
  //
  product = ndivs[0]*ndivs[1];
  for (i=0; i<numPts; i++) 
    {
    x = this->DataSet->GetPoint(i);
    for (j=0; j<3; j++) 
      {
      ijk[j] = (int) ((double) ((x[j] - this->Bounds[2*j]) / 
                        (this->Bounds[2*j+1] - this->Bounds[2*j])) * ndivs[j]);
      if (ijk[j] >= this->Divisions[j])
        {
        ijk[j] = this->Divisions[j] - 1;
        }
      }

    idx = ijk[0] + ijk[1]*ndivs[0] + ijk[2]*product;
    bucket = this->HashTable[idx];
    if ( ! bucket )
      {
      bucket = vtkIdList::New();
      bucket->Allocate(this->NumberOfPointsPerBucket,
                       this->NumberOfPointsPerBucket/3);
      this->HashTable[idx] = bucket;
      }
    bucket->InsertNextId(i);
    }

  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);
          }
        }
      }
    }

  return;
}


//
// 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];

  // Initialize
  buckets->Reset();

  // Determine the range of indices in each direction
  for (i=0; i < 3; i++)
    {
    minLevel[i] = (int) ((double) (((x[i]-dist) - this->Bounds[2*i]) / 
        (this->Bounds[2*i+1] - this->Bounds[2*i])) * this->Divisions[i]);
    maxLevel[i] = (int) ((double) (((x[i]+dist) - this->Bounds[2*i]) / 
        (this->Bounds[2*i+1] - this->Bounds[2*i])) * this->Divisions[i]);

    if ( minLevel[i] < 0 )
      {
      minLevel[i] = 0;
      }
    if ( maxLevel[i] >= this->Divisions[i] )
      {
      maxLevel[i] = this->Divisions[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], double dist, 
                                            int prevMinLevel[3],
                                            int prevMaxLevel[3])
{
  int i, j, k, nei[3], minLevel[3], maxLevel[3];
  int kFactor, jFactor;
  int jkSkipFlag, kSkipFlag;

  // Initialize
  buckets->Reset();

  // Determine the range of indices in each direction
  for (i=0; i < 3; i++)
    {
    minLevel[i] = (int) ((double) (((x[i]-dist) - this->Bounds[2*i])
                                  / this->H[i]));
    maxLevel[i] = (int) ((double) (((x[i]+dist) - this->Bounds[2*i])
                                  / this->H[i]));
    
    if ( minLevel[i] < 0 )
      {
      minLevel[i] = 0;
      }
    else if (minLevel[i] >= this->Divisions[i] )
      {
      minLevel[i] = this->Divisions[i] - 1;
      }
    if ( maxLevel[i] >= this->Divisions[i] )
      {
      maxLevel[i] = this->Divisions[i] - 1;
      }
    else if ( maxLevel[i] < 0 )
      {
      maxLevel[i] = 0;
      }
    }


  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->Divisions[0]*this->Divisions[1];
    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->Divisions[0];
      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];
  double level;

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

  for (i=0; i<3; i++)
    {
    this->Bounds[2*i] = bounds[2*i];
    this->Bounds[2*i+1] = bounds[2*i+1];
    if ( this->Bounds[2*i+1] <= this->Bounds[2*i] )
      {
      this->Bounds[2*i+1] = this->Bounds[2*i] + 1.0;
      }
    }

  if ( this->Automatic && (estNumPts > 0) )
    {
    level = (double) estNumPts / this->NumberOfPointsPerBucket;
    level = ceil( pow((double)level,(double)0.33333333) );
    for (i=0; i<3; i++)
      {
      ndivs[i] = (int) level;
      }
    } 
  else 
    {
    for (i=0; i<3; i++)
      {
      ndivs[i] = (int) this->Divisions[i];
      }
    }

  for (i=0; i<3; i++) 
    {
    ndivs[i] = (ndivs[i] > 0 ? ndivs[i] : 1);
    this->Divisions[i] = ndivs[i];
    }

  this->NumberOfBuckets = ndivs[0]*ndivs[1]*ndivs[2];
  this->HashTable = new vtkIdListPtr[this->NumberOfBuckets];
  memset (this->HashTable, 0, this->NumberOfBuckets*
          sizeof(vtkIdListPtr));
  //
  //  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 ((double) this->Tolerance / hmin);
  this->InsertionLevel = (this->InsertionLevel > maxDivs ? maxDivs : this->InsertionLevel);
  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 dulicate 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])
{
  int i, ijk[3];
  vtkIdType idx;
  vtkIdList *bucket;
  //
  //  Locate bucket that point is in.
  //
  for (i=0; i<3; i++)
    {
    ijk[i] = (int) ((double) ((x[i] - this->Bounds[2*i]) / 
        (this->Bounds[2*i+1] - this->Bounds[2*i])) * this->Divisions[i]);
    if (ijk[i] >= this->Divisions[i])
      {
      ijk[i] = this->Divisions[i] - 1;
      }
    }

  idx = ijk[0] + ijk[1]*this->Divisions[0] + 
        ijk[2]*this->Divisions[0]*this->Divisions[1];

  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
// dulicate 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])
{
  int i, ijk[3];
  vtkIdType idx;
  vtkIdList *bucket;
  //
  //  Locate bucket that point is in.
  //
  for (i=0; i<3; i++)
    {
    ijk[i] = (int) ((double) ((x[i] - this->Bounds[2*i]) / 
       (this->Bounds[2*i+1] - this->Bounds[2*i])) * this->Divisions[i]);
    if (ijk[i] >= this->Divisions[i])
      {
      ijk[i] = this->Divisions[i] - 1;
      }
    }

  idx = ijk[0] + ijk[1]*this->Divisions[0] + 
        ijk[2]*this->Divisions[0]*this->Divisions[1];

  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;

  //  Locate bucket that point is in.
  //
  for (i=0; i<3; i++)
    {
    ijk[i] = (int) ((double) ((x[i] - this->Bounds[2*i]) / 
        (this->Bounds[2*i+1] - this->Bounds[2*i])) * this->Divisions[i]);
    if (ijk[i] >= this->Divisions[i])
      {
      ijk[i] = this->Divisions[i] - 1;
      }
    }

  // 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->Divisions[0] + 
        nei[2]*this->Divisions[0]*this->Divisions[1];

      if ( (ptIds = this->HashTable[cno]) != NULL )
        {
        for (j=0; j < ptIds->GetNumberOfIds(); 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, 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.  
  //
  for (j=0; j<3; j++) 
    {
    ijk[j] = (int)(((x[j] - this->Bounds[2*j]) / 
        (this->Bounds[2*j+1] - this->Bounds[2*j])) * this->Divisions[j]);
    if (ijk[j] >= this->Divisions[j])
      {
      ijk[j] = this->Divisions[j] - 1;
      }
    }

  //  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=0,minDist2=VTK_DOUBLE_MAX,level=0; (closest == 0) && 
  (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->Divisions[0] + 
            nei[2]*this->Divisions[0]*this->Divisions[1];

      if ( (ptIds = this->HashTable[cno]) != NULL )
        {
        for (j=0; j < ptIds->GetNumberOfIds(); 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->Divisions[0] + nei[2]*this->Divisions[0]*this->Divisions[1];

        if ( (ptIds = this->HashTable[cno]) )
          {
          for (j=0; j < ptIds->GetNumberOfIds(); 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 NULL;
      }
    }

  //  Find bucket point is in.  
  //
  for (i=0; i<3; i++) 
    {
    ijk[i] = (int)(((x[i] - this->Bounds[2*i]) / 
        (this->Bounds[2*i+1] - this->Bounds[2*i])) * this->Divisions[i]);
    if (ijk[i] >= this->Divisions[i])
      {
      ijk[i] = this->Divisions[i] - 1;
      }
    }
  
  // Get the id list, if any
  //
  if ( this->HashTable )
    {
    int idx = ijk[0] + ijk[1]*this->Divisions[0] + 
              ijk[2]*this->Divisions[0]*this->Divisions[1];

    return this->HashTable[idx];
    }

  return NULL;
}


// 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 == NULL ) 
    {
    vtkErrorMacro(<<"Can't build representation...no data!");
    return;
    }

  pts = vtkPoints::New();
  pts->Allocate(5000);
  polys = vtkCellArray::New();
  polys->Allocate(10000);

  // 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] == NULL )
          {
          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] == NULL && inside) ||
            (this->HashTable[idx] != NULL && !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->H[0];
  origin[1] = this->Bounds[2] + j * this->H[1];
  origin[2] = this->Bounds[4] + k * this->H[2];
  ids[0] = pts->InsertNextPoint(origin);

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

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

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

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

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

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

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

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

    x[0] = origin[0];
    x[1] = origin[1] + this->H[1];
    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->H[0] + this->Bounds[0];
  bounds[1] = (nei[0]+1)*this->H[0] + this->Bounds[0];
  bounds[2] =     nei[1]*this->H[1] + this->Bounds[2];
  bounds[3] = (nei[1]+1)*this->H[1] + this->Bounds[2];
  bounds[4] =     nei[2]*this->H[2] + this->Bounds[4];
  bounds[5] = (nei[2]+1)*this->H[2] + this->Bounds[4];

  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";
    }
}