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

  Program:   Visualization Toolkit
  Module:    vtkPointLocator.cxx
  Language:  C++
  Date:      $Date$
  Version:   $Revision$


Copyright (c) 1993-2001 Ken Martin, Will Schroeder, Bill Lorensen 
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

 * Redistributions of source code must retain the above copyright notice,
   this list of conditions and the following disclaimer.

 * Redistributions in binary form must reproduce the above copyright notice,
   this list of conditions and the following disclaimer in the documentation
   and/or other materials provided with the distribution.

 * Neither name of Ken Martin, Will Schroeder, or Bill Lorensen nor the names
   of any contributors may be used to endorse or promote products derived
   from this software without specific prior written permission.

 * Modified source versions must be plainly marked as such, and must not be
   misrepresented as being the original software.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ``AS IS''
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE AUTHORS OR CONTRIBUTORS BE LIABLE FOR
ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

=========================================================================*/
#include "vtkPointLocator.h"
#include "vtkMath.h"
#include "vtkIntArray.h"
#include "vtkPolyData.h"
#include "vtkObjectFactory.h"
#include "vtkIdList.h"

//------------------------------------------------------------------------------
vtkPointLocator* vtkPointLocator::New()
{
  // First try to create the object from the vtkObjectFactory
  vtkObject* ret = vtkObjectFactory::CreateInstance("vtkPointLocator");
  if(ret)
    {
    return (vtkPointLocator*)ret;
    }
  // If the factory was unable to create the object, then create it here.
  return new 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(float x, float y, float z)
{
  float 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 float x[3])
{
  int i, j;
  float minDist2;
  float dist2 = VTK_LARGE_FLOAT;
  float *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_LARGE_FLOAT,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(float radius, 
                                                        const float x[3],
                                                        float& dist2)
  {
  return FindClosestPointWithinRadius(radius, x, this->DataSet->GetLength(),
                                      dist2);
  }


vtkIdType vtkPointLocator::FindClosestPointWithinRadius(float radius,
                                                        const float x[3], 
                                                        float inputDataLength,
                                                        float& dist2)
{
  int i, j;
  float *pt;
  vtkIdType ptId, closest = -1;
  vtkIdList *ptIds;
  int ijk[3], *nei;
  float minDist2;
  
  float refinedRadius, radius2, refinedRadius2;
  float currentRadius;
  float 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)((float)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;
  float dist;
};

static int idsortcompare(const void *arg1, const void *arg2)
{
  idsort *v1 = (idsort *)arg1;
  idsort *v2 = (idsort *)arg2;
  if (v1->dist < v2->dist)
    {
    return -1;
    }
  if (v1->dist > v2->dist)
    {
    return 1;
    }
  return 0;  
}

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

static int GetOctent(const float x[3], const float pt[3])
{
  float 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 float GetMax(const float foo[8])
{
  float 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 float x[3],
                                            vtkIdList *result, int M)
{
  int i, j;
  float dist2;
  float *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;
  float 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),idsortcompare);
              }
            }
          else if (dist2 < maxDistance[oct])
            {
            res[oct][N-1].dist = dist2;
            res[oct][N-1].id = ptId;
            qsort(res[oct], N, sizeof(idsort), idsortcompare);
            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), idsortcompare);
    }
  
  // 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), idsortcompare);
          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, float x,
                                         float y, float z,
                                         vtkIdList *result)
{
  float p[3];
  p[0] = x;
  p[1] = y;
  p[2] = z;
  this->FindClosestNPoints(N,p,result);
}

void vtkPointLocator::FindClosestNPoints(int N, const float x[3],
                                         vtkIdList *result)
{
  int i, j;
  float dist2;
  float *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;
  float 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), idsortcompare);
              }
            }
          else if (dist2 < maxDistance)
            {
            res[N-1].dist = dist2;
            res[N-1].id = ptId;
            qsort(res, N, sizeof(idsort), idsortcompare);
            maxDistance = res[N-1].dist;
            }
          }
        }
      }
    level++;
    this->GetBucketNeighbors (&buckets, ijk, this->Divisions, level);
    }

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

  // 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), idsortcompare);
          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(float R, float x,
                                             float y, float z,
                                             vtkIdList *result)
{
  float p[3];
  p[0] = x;
  p[1] = y;
  p[2] = z;
  this->FindPointsWithinRadius(R,p,result);
}



void vtkPointLocator::FindPointsWithinRadius(float R, const float x[3],
                                             vtkIdList *result)
{
  int i, j;
  float dist2;
  float *pt;
  vtkIdType ptId, cno;
  vtkIdList *ptIds;
  int ijk[3], *nei;
  float 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()
{
  float *bounds;
  vtkIdType numBuckets;
  float level;
  int ndivs[3], product;
  int i, j, ijk[3];
  vtkIdType idx;
  vtkIdList *bucket;
  vtkIdType numPts;
  float *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 = (float) 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) ((float) ((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 float x[3],
                                            const int ijk[3], 
                                            float 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) ((float) (((x[i]-dist) - this->Bounds[2*i]) / 
        (this->Bounds[2*i+1] - this->Bounds[2*i])) * this->Divisions[i]);
    maxLevel[i] = (int) ((float) (((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 float x[3], float 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) ((float) (((x[i]-dist) - this->Bounds[2*i])
                                  / this->H[i]));
    maxLevel[i] = (int) ((float) (((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 float 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 float bounds[6],
                                        vtkIdType estNumPts)
{
  int i;
  int maxDivs;
  typedef vtkIdList *vtkIdListPtr;
  float hmin;
  int ndivs[3];
  float 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 = (float) 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_LARGE_FLOAT, 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 float 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) ((float) ((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 float 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) ((float) ((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 float x[3])
{
  int i, j, ijk[3];
  vtkNeighborPoints buckets;

  //  Locate bucket that point is in.
  //
  for (i=0; i<3; i++)
    {
    ijk[i] = (int) ((float) ((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;
  float *pt;

  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);
          pt = this->Points->GetPoint(ptId);

          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 float 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 float x[3])
{
  int i;
  float minDist2, dist2;
  float *pt;
  int level;
  vtkIdType closest, j;
  vtkIdType ptId, cno;
  vtkIdList *ptIds;
  int ijk[3], *nei;
  int MULTIPLES;
  float 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_LARGE_FLOAT,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);
          pt = this->Points->GetPoint(ptId);
          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);
            pt = this->Points->GetPoint(ptId);
            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 float 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];
  float 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!!!!
//
float vtkPointLocator::Distance2ToBucket(const float x[3],
                                         const int nei[3])
{
  float 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!!!!
float vtkPointLocator::Distance2ToBounds(const float x[3],
                                         const float bounds[6])
{
  float distance;
  float 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)
{
  vtkLocator::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";

}

