vtkOctreePointLocatorNode.cxx

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

  Program:   Visualization Toolkit
  Module:    vtkOctreePointLocatorNode.cxx

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

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

=========================================================================*/
/*----------------------------------------------------------------------------
 Copyright (c) Sandia Corporation
 See Copyright.txt or http://www.paraview.org/HTML/Copyright.html for details.
----------------------------------------------------------------------------*/

#include "vtkOctreePointLocatorNode.h"

#include "vtkBox.h"
#include "vtkCell.h"
#include "vtkMath.h"
#include "vtkObjectFactory.h"
#include "vtkPlanesIntersection.h"
#include "vtkPoints.h"

vtkStandardNewMacro(vtkOctreePointLocatorNode);

//----------------------------------------------------------------------------
vtkOctreePointLocatorNode::vtkOctreePointLocatorNode() : 
  NumberOfPoints(0), Children(NULL), ID(-1), MinID(-1)
{
  // set the min and max data value since we won't know it for a while
  for(int i=0;i<3;i++)
    {
    this->MinDataBounds[i] = VTK_DOUBLE_MAX;
    this->MaxDataBounds[i] = VTK_DOUBLE_MIN;
    }
}

//----------------------------------------------------------------------------
vtkOctreePointLocatorNode::~vtkOctreePointLocatorNode()
{
  this->DeleteChildNodes();
}

//----------------------------------------------------------------------------
void vtkOctreePointLocatorNode::CreateChildNodes()
{
  if(!this->Children)
    {
    int i;
    double midpoint[3];
    for(i=0;i<3;i++)
      {
      midpoint[i] = (this->MinBounds[i]+this->MaxBounds[i])*.5;
      }
    this->Children = new vtkOctreePointLocatorNode*[8];
    for(i=0;i<8;i++)
      {
      this->Children[i] = vtkOctreePointLocatorNode::New();
      double NewMin[3], NewMax[3];
      for(int j=0;j<3;j++)
        {
        if(!((i>>j) & 1))
          {
          NewMin[j] = this->MinBounds[j];
          NewMax[j] = midpoint[j];
          }
        else
          {
          NewMin[j] = midpoint[j];
          NewMax[j] = this->MaxBounds[j];
          }
        }
      this->Children[i]->SetMinBounds(NewMin);
      this->Children[i]->SetMaxBounds(NewMax);
      }
    }
}

//----------------------------------------------------------------------------
void vtkOctreePointLocatorNode::DeleteChildNodes()
{     
  if(this->Children)
    {
    for(int i=0;i<8;i++)
      {
      this->Children[i]->Delete();
      }
    delete []this->Children;
    this->Children = 0;
    }
}       

//----------------------------------------------------------------------------
vtkOctreePointLocatorNode* vtkOctreePointLocatorNode::GetChild(int i)
{
  if(this->Children)
    {
    return this->Children[i];
    }
  return 0;
}

//----------------------------------------------------------------------------
int vtkOctreePointLocatorNode::GetSubOctantIndex(
  double* point, int CheckContainment)
{
  int i, index = 0;
  if(CheckContainment)
    {
    for(i=0;i<3;i++)
      {
      if(point[i] <= this->MinBounds[i] || point[i] > this->MaxBounds[i])
        {
        return -1;
        }
      }
    }

  for(i=0;i<3;i++)
    {
    if(point[i] > (this->MinBounds[i]+this->MaxBounds[i])*.5)
      {
      index += 1 << i;
      }
    }
  return index;
}

//----------------------------------------------------------------------------
void vtkOctreePointLocatorNode::ComputeOctreeNodeInformation(
  vtkOctreePointLocatorNode* Parent, int & NextLeafId, int & NextMinId, 
  float* coordinates)
{
  this->MinID = NextMinId;
  if(this->Children)
    {
    int i;
    for(i=0;i<8;i++)
      {
      this->Children[i]->ComputeOctreeNodeInformation(
        this, NextLeafId, NextMinId, coordinates);
      }
    // a non-leaf region can get its data bounds from its children...
    this->SetMinDataBounds(this->Children[0]->GetMinDataBounds());
    this->SetMaxDataBounds(this->Children[0]->GetMaxDataBounds());

    for(i=1;i<8;i++)
      {
      double* min = this->Children[i]->GetMinDataBounds();
      double* max = this->Children[i]->GetMaxDataBounds();
      for(int j=0;j<3;j++)
        {
        if(min[j] < this->MinDataBounds[j])
          {
          this->MinDataBounds[j] = min[j];
          }
        if(max[j] > this->MaxDataBounds[j])
          {
          this->MaxDataBounds[j] = max[j];
          }
        }
      }
    }
  else
    {
    this->ID = NextLeafId;
    NextLeafId++;
    NextMinId = this->MinID+this->NumberOfPoints;
    if(this->NumberOfPoints == 0)
      {
      // since there are no points in this, set the data bounds
      // such that they won't affect anything else
      this->SetMinDataBounds(Parent->GetMaxBounds());
      this->SetMaxDataBounds(Parent->GetMinBounds());
      }
    else
      {
      int i;
      float* coordptr = coordinates+3*this->MinID;
      for(i=0;i<3;i++)
        {
        this->MinDataBounds[i] = this->MaxDataBounds[i] = coordptr[i];
        }
      for(i=1;i<this->NumberOfPoints;i++)
        {
        coordptr += 3;
        for(int j=0;j<3;j++)
          {
          if(coordptr[j] < this->MinDataBounds[j])
            {
            this->MinDataBounds[j] = coordptr[j];
            }
          else if(coordptr[j] > this->MaxDataBounds[j])
            {
            this->MaxDataBounds[j] = coordptr[j];
            }
          }
        }
      }
    }
}

//----------------------------------------------------------------------------
void vtkOctreePointLocatorNode::SetBounds(
  double xMin,double xMax,double yMin,double yMax,double zMin,double zMax)
{
   this->MinBounds[0] = xMin; this->MaxBounds[0] = xMax;
   this->MinBounds[1] = yMin; this->MaxBounds[1] = yMax;
   this->MinBounds[2] = zMin; this->MaxBounds[2] = zMax;
}

//----------------------------------------------------------------------------
void vtkOctreePointLocatorNode::GetBounds(double *b) const
{
   b[0] = this->MinBounds[0]; b[1] = this->MaxBounds[0];
   b[2] = this->MinBounds[1]; b[3] = this->MaxBounds[1];
   b[4] = this->MinBounds[2]; b[5] = this->MaxBounds[2];
}

//----------------------------------------------------------------------------
void vtkOctreePointLocatorNode::SetDataBounds(
  double xMin,double xMax,double yMin,double yMax,double zMin,double zMax)
{
   this->MinDataBounds[0] = xMin; this->MaxDataBounds[0] = xMax;
   this->MinDataBounds[1] = yMin; this->MaxDataBounds[1] = yMax;
   this->MinDataBounds[2] = zMin; this->MaxDataBounds[2] = zMax;
}

//----------------------------------------------------------------------------
void vtkOctreePointLocatorNode::GetDataBounds(double *b) const
{
   b[0] = this->MinDataBounds[0]; b[1] = this->MaxDataBounds[0];
   b[2] = this->MinDataBounds[1]; b[3] = this->MaxDataBounds[1];
   b[4] = this->MinDataBounds[2]; b[5] = this->MaxDataBounds[2];
}

//----------------------------------------------------------------------------
// Squared distance from any point anywhere to the boundary of spatial region
//   
double vtkOctreePointLocatorNode::GetDistance2ToBoundary(
  double x, double y, double z, vtkOctreePointLocatorNode* top,
  int useDataBounds=0)
{
  return this->_GetDistance2ToBoundary(x, y, z, NULL, 0, top, useDataBounds);
}
 
//----------------------------------------------------------------------------
// Squared distance from any point anywhere to the boundary of spatial region,
// and give me the point on the boundary closest to this point.
//
double vtkOctreePointLocatorNode::GetDistance2ToBoundary(
  double x, double y, double z, double *p,
  vtkOctreePointLocatorNode* top, int useDataBounds=0)
{
  return this->_GetDistance2ToBoundary(x, y, z, p, 0, top, useDataBounds);
}

//----------------------------------------------------------------------------
// The point is inside the region, and I want the distance (squared)
// to the closest "interior" wall, one that is not an outer boundary of
// the entire space.
//
double vtkOctreePointLocatorNode::GetDistance2ToInnerBoundary(
  double x, double y, double z, vtkOctreePointLocatorNode* top)
{
  return this->_GetDistance2ToBoundary(x, y, z, NULL, 1, top, 0);
}

//----------------------------------------------------------------------------
double vtkOctreePointLocatorNode::_GetDistance2ToBoundary( 
  double x, double y, double z,  // from this point
  double *p,              // set to point on boundary that is closest
  int innerBoundaryOnly, // ignore boundaries on "outside"
  vtkOctreePointLocatorNode* top,    // the top of the octree
  int useDataBounds=0)   // use bounds of data within region instead
{
  double minDistance, dist;
  double edgePt[3];
  double cornerPt[3];
  double pt3[3];

  double *min, *max;

  if (useDataBounds)
    {
    min = this->MinDataBounds; max = this->MaxDataBounds;  // data inside region
    }
  else
    {
    min = this->MinBounds; max = this->MaxBounds;   // region itself
    }

  double *outerBoundaryMin=NULL;
  double *outerBoundaryMax=NULL;

  if (innerBoundaryOnly)
    {
    // We want the distance to the nearest inner boundary, because we 
    // are only interested in boundaries such that there may be points on
    // the other side.  This option only makes sense when the point supplied
    // is inside this node (region).

    outerBoundaryMin = (useDataBounds? top->MinDataBounds : top->MinBounds);
    outerBoundaryMax = (useDataBounds? top->MaxDataBounds : top->MaxBounds);
    }

  double xmax = max[0]; double ymax = max[1]; double zmax = max[2];
  double xmin = min[0]; double ymin = min[1]; double zmin = min[2];

  int xless = (x < xmin);
  int xmore = (x > xmax);
  int yless = (y < ymin);
  int ymore = (y > ymax);
  int zless = (z < zmin);
  int zmore = (z > zmax);

  int withinX = !xless && !xmore;
  int withinY = !yless && !ymore;
  int withinZ = !zless && !zmore;

  int mindim=0;

  if (withinX && withinY && withinZ)  // point is inside the box
    {
    if (!innerBoundaryOnly)
      {
      minDistance = x - xmin;
      mindim = 0;
  
      if ((dist = xmax - x) < minDistance)
        {
        mindim = 1;
        minDistance = dist;
        }
      if ((dist = y - ymin) < minDistance)
        {
        mindim = 2;
        minDistance = dist;
        }
      if ((dist = ymax - y) < minDistance)
        {
        mindim = 3;
        minDistance = dist;
        }
      if ((dist = z - zmin) < minDistance)
        {
        mindim = 4;
        minDistance = dist;
        }
      if ((dist = zmax - z) < minDistance)
        {
        mindim = 5;
        minDistance = dist;
        }
      }
    else
      {
      int first = 1;
      minDistance = VTK_LARGE_FLOAT; // Suppresses warning message.

      if ((xmin != outerBoundaryMin[0]) && 
          (((dist = x - xmin) < minDistance) || first))
        {
        mindim = 0;
        minDistance = dist;
        first = 0;
        }
      if ((xmax != outerBoundaryMax[0]) && 
          (((dist = xmax - x) < minDistance) || first))
        {
        mindim = 1;
        minDistance = dist;
        first = 0;
        }
      if ((ymin != outerBoundaryMin[1]) && 
          (((dist = y - ymin) < minDistance) || first))
        {
        mindim = 2;
        minDistance = dist;
        first = 0;
        }
      if ((ymax != outerBoundaryMax[1]) && 
          (((dist = ymax - y) < minDistance) || first))
        {
        mindim = 3;
        minDistance = dist;
        first = 0;
        }
      if ((zmin != outerBoundaryMin[2]) && 
          (((dist = z - zmin) < minDistance) || first))
        {
        mindim = 4;
        minDistance = dist;
        first = 0;
        }
      if ((zmax != outerBoundaryMax[2]) && 
          (((dist = zmax - z) < minDistance) || first))
        {
        mindim = 5;
        minDistance = dist;
        }
      }

    // if there are no inner boundaries we dont want to square.
    if(minDistance != VTK_LARGE_FLOAT)
      {
      minDistance *= minDistance;
      }

    if (p)
      {
      p[0] = x; p[1] = y; p[2] = z;


      if (mindim == 0)      
        {
        p[0] = xmin;
        }
      else if (mindim == 1) 
        {
        p[0] = xmax;
        }
      else if (mindim == 2) 
        {
        p[1] = ymin;
        }
      else if (mindim == 3) 
        {
        p[1] = ymax;
        }
      else if (mindim == 4) 
        {
        p[2] = zmin;
        }
      else if (mindim == 5) 
        {
        p[2] = zmax;
        }
      }
    }
  else if (withinX && withinY)  // point projects orthogonally to a face
    {
    minDistance = (zless ? zmin - z : z - zmax);
    minDistance *= minDistance;

    if (p)
      {
      p[0] = x; p[1] = y;
      p[2] = (zless ? zmin : zmax);
      }
    }
  else if (withinX && withinZ)
    {
    minDistance = (yless ? ymin - y : y - ymax);
    minDistance *= minDistance;

    if (p)
      {
      p[0] = x; p[2] = z;
      p[1] = (yless ? ymin : ymax);
      }
    }
  else if (withinY && withinZ)
    {
    minDistance = (xless ? xmin - x : x - xmax);
    minDistance *= minDistance;

    if (p)
      {
      p[1] = y; p[2] = z;
      p[0] = (xless ? xmin: xmax);
      }
    }
  else if (withinX || withinY || withinZ)   // point is closest to an edge
    {
    edgePt[0] = (withinX ? x : (xless ? xmin : xmax));
    edgePt[1] = (withinY ? y : (yless ? ymin : ymax));
    edgePt[2] = (withinZ ? z : (zless ? zmin : zmax));

    pt3[0] = x; pt3[1] = y; pt3[2] = z;

    minDistance = vtkMath::Distance2BetweenPoints(pt3, edgePt);

    if (p)
      {
      p[0] = edgePt[0]; p[1] = edgePt[1]; p[2] = edgePt[2];
      }
    }
  else                        // point is closest to a corner
    {
    cornerPt[0] = (xless ? xmin : xmax);
    cornerPt[1] = (yless ? ymin : ymax);
    cornerPt[2] = (zless ? zmin : zmax);

    pt3[0] = x; pt3[1] = y; pt3[2] = z;

    minDistance = vtkMath::Distance2BetweenPoints(pt3, cornerPt);

    if (p)
      {
      p[0] = cornerPt[0]; p[1] = cornerPt[1]; p[2] = cornerPt[2];
      }
    }

  return minDistance;
}

//----------------------------------------------------------------------------
int vtkOctreePointLocatorNode::ContainsPoint(double x, double y, double z, 
                                 int useDataBounds=0)
{
  double *min, *max;

  if (useDataBounds)
    {
    min = this->MinDataBounds;
    max = this->MaxDataBounds;
    }
  else
    {
    min = this->MinBounds;
    max = this->MaxBounds;
    }

  if( (min[0] >= x) ||
      (max[0] < x) ||
      (min[1] >= y) ||
      (max[1] < y) ||
      (min[2] >= z) ||
      (max[2] < z))
    {
    return 0;
    }
  else
    {
    return 1;
    }
}

//----------------------------------------------------------------------------
int vtkOctreePointLocatorNode::IntersectsRegion(
  vtkPlanesIntersection *pi, int useDataBounds)
{
  double xMin, xMax, yMin, yMax, zMin, zMax;
  vtkPoints *box = vtkPoints::New();

  box->SetNumberOfPoints(8);

  double *min, *max;

  if (useDataBounds)
    {
    min = this->MinDataBounds;
    max = this->MaxDataBounds;
    }
  else
    {
    min = this->MinBounds;
    max = this->MaxBounds;
    }

  xMin = min[0]; xMax = max[0];
  yMin = min[1]; yMax = max[1];
  zMin = min[2]; zMax = max[2];

  box->SetPoint(0, xMax, yMin, zMax);
  box->SetPoint(1, xMax, yMin, zMin);
  box->SetPoint(2, xMax, yMax, zMin);
  box->SetPoint(3, xMax, yMax, zMax);
  box->SetPoint(4, xMin, yMin, zMax);
  box->SetPoint(5, xMin, yMin, zMin);
  box->SetPoint(6, xMin, yMax, zMin);
  box->SetPoint(7, xMin, yMax, zMax);

  int intersects = pi->IntersectsRegion(box);

  box->Delete();

  return intersects;
}

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

  os << indent << "NumberOfPoints: " << this->NumberOfPoints << endl;
  os << indent << "Children: " << this->Children << endl;
  os << indent << "ID: " << this->ID << endl;
  os << indent << "MinID: " << this->MinID << endl;
  os << indent << "MinBounds: " << this->MinBounds[0] << " " 
     << this->MinBounds[1] << " " << this->MinBounds[2] << endl;
  os << indent << "MaxBounds: " << this->MaxBounds[0] << " " 
     << this->MaxBounds[1] << " " << this->MaxBounds[2] << endl;
  os << indent << "MinDataBounds: " << this->MinDataBounds[0] << " " 
     << this->MinDataBounds[1] << " " << this->MinDataBounds[2] << endl;
  os << indent << "MaxDataBounds: " << this->MaxDataBounds[0] << " " 
     << this->MaxDataBounds[1] << " " << this->MaxDataBounds[2] << endl;
}