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

  Program:   Visualization Toolkit
  Module:    vtkIncrementalOctreeNode.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.

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

// Hide VTK_DEPRECATED_IN_9_1_0() warnings for this class.
#define VTK_DEPRECATION_LEVEL 0

#include "vtkIncrementalOctreeNode.h"
#include "vtkIdList.h"
#include "vtkMath.h"
#include "vtkObjectFactory.h"
#include "vtkPoints.h"

vtkStandardNewMacro(vtkIncrementalOctreeNode);

vtkCxxSetObjectMacro(vtkIncrementalOctreeNode, PointIdSet, vtkIdList);
vtkCxxSetObjectMacro(vtkIncrementalOctreeNode, Parent, vtkIncrementalOctreeNode);

//------------------------------------------------------------------------------
// ----------------------------- Helper functions ----------------------------
//------------------------------------------------------------------------------

//------------------------------------------------------------------------------
// This is an empty function which provides only the point index to the caller
// function vtkIncreemntalOctreeNode::InsertPoint(). The caller inserts the
// point index to the vtkIdList maintained by a leaf node, without inserting
// the point (coordinate) to vtkPoints at all.
static void OctreeNodeGetPointId(
  vtkPoints* vtkNotUsed(points), vtkIdType* vtkNotUsed(pntIdx), const double* vtkNotUsed(coords))
{
  // the 3D point coordinate is not inserted to vtkPoints at all
}

//------------------------------------------------------------------------------
// Insert a point, with a specified point index, to a vtkPoints object by
// calling vtkPoints::InsertPoint().
static void OctreeNodeInsertPoint(vtkPoints* points, vtkIdType* pntIdx, const double* coords)
{
  points->InsertPoint(*pntIdx, coords);
}

//------------------------------------------------------------------------------
// Insert a point to a vtkPoints by calling vtkPoints::InsertNextPoint().
static void OctreeNodeInsertNextPoint(vtkPoints* points, vtkIdType* pntIdx, const double* coords)
{
  *pntIdx = points->InsertNextPoint(coords);
}

//------------------------------------------------------------------------------
// Function pointers in support of three point insertion modes.
typedef void (*OCTREENODE_INSERTPOINT_FUNCTION)(
  vtkPoints* points, vtkIdType* pntIdx, const double* coords);

static OCTREENODE_INSERTPOINT_FUNCTION OCTREENODE_INSERTPOINT[3] = { OctreeNodeGetPointId,
  OctreeNodeInsertPoint, OctreeNodeInsertNextPoint };

//------------------------------------------------------------------------------
// ------------------------- vtkIncrementalOctreeNode ------------------------
//------------------------------------------------------------------------------

//------------------------------------------------------------------------------
vtkIncrementalOctreeNode::vtkIncrementalOctreeNode()
{
  this->Parent = nullptr;
  this->Children = nullptr;
  this->PointIdSet = nullptr;
  this->NumberOfPoints = 0;
  this->ID = 0;

  // unnecessary to initialize spatial and data bounding boxes here as
  // SetBounds() are always called by the user for the root node of an
  // octree and this function is also always called for any descendant
  // node upon node-subdivison.  for now though we'll set them to
  // something to avoid uninitialized variable warnings
  for (int i = 0; i < 3; i++)
  {
    this->MinBounds[i] = this->MinDataBounds[i] = VTK_DOUBLE_MIN;
    this->MaxBounds[i] = this->MaxDataBounds[i] = VTK_DOUBLE_MAX;
  }
}

//------------------------------------------------------------------------------
vtkIncrementalOctreeNode::~vtkIncrementalOctreeNode()
{
  if (this->Parent)
  {
    this->Parent->UnRegister(this);
    this->Parent = nullptr;
  }

  this->DeleteChildNodes();
  this->DeletePointIdSet();
}

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

    delete[] this->Children;
    this->Children = nullptr;
  }
}

//------------------------------------------------------------------------------
void vtkIncrementalOctreeNode::CreatePointIdSet(int initSize, int growSize)
{
  if (this->PointIdSet == nullptr)
  {
    this->PointIdSet = vtkIdList::New();
    this->PointIdSet->Allocate(initSize, growSize);
  }
}

//------------------------------------------------------------------------------
void vtkIncrementalOctreeNode::DeletePointIdSet()
{
  if (this->PointIdSet)
  {
    this->PointIdSet->Delete();
    this->PointIdSet = nullptr;
  }
}

//------------------------------------------------------------------------------
void vtkIncrementalOctreeNode::SetBounds(
  double x1, double x2, double y1, double y2, double z1, double z2)
{
  this->MinBounds[0] = x1;
  this->MaxBounds[0] = x2;
  this->MinBounds[1] = y1;
  this->MaxBounds[1] = y2;
  this->MinBounds[2] = z1;
  this->MaxBounds[2] = z2;

  this->MinDataBounds[0] = x2;
  this->MaxDataBounds[0] = x1;
  this->MinDataBounds[1] = y2;
  this->MaxDataBounds[1] = y1;
  this->MinDataBounds[2] = z2;
  this->MaxDataBounds[2] = z1;
}

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

//------------------------------------------------------------------------------
int vtkIncrementalOctreeNode::UpdateCounterAndDataBounds(
  const double point[3], int nHits, int updateData)
{
  this->NumberOfPoints += nHits;

  if (updateData == 0)
    return 0;

  int updated = 0;

  if (point[0] < this->MinDataBounds[0])
  {
    updated = 1;
    this->MinDataBounds[0] = point[0];
  }
  if (point[0] > this->MaxDataBounds[0])
  {
    updated = 1;
    this->MaxDataBounds[0] = point[0];
  }

  if (point[1] < this->MinDataBounds[1])
  {
    updated = 1;
    this->MinDataBounds[1] = point[1];
  }
  if (point[1] > this->MaxDataBounds[1])
  {
    updated = 1;
    this->MaxDataBounds[1] = point[1];
  }

  if (point[2] < this->MinDataBounds[2])
  {
    updated = 1;
    this->MinDataBounds[2] = point[2];
  }
  if (point[2] > this->MaxDataBounds[2])
  {
    updated = 1;
    this->MaxDataBounds[2] = point[2];
  }

  return updated;
}

//------------------------------------------------------------------------------
// Given the index (0 ~ 7) of a child node, the spatial bounding axis (0 ~ 2
// for x, y, and z), and the value (0 ~ 1 for min and max) to access, this LUT
// allows for rapid assignment of its spatial bounding box --- MinBounds[3]
// and MaxBounds[3], with each specific value or entry of this LUT pointing to
// MinBounds[3] for 0, center point for 1, or MaxBounds[3] for 2.
static int OCTREE_CHILD_BOUNDS_LUT[8][3][2] = {
  { { 0, 1 }, { 0, 1 }, { 0, 1 } },
  { { 1, 2 }, { 0, 1 }, { 0, 1 } },
  { { 0, 1 }, { 1, 2 }, { 0, 1 } },
  { { 1, 2 }, { 1, 2 }, { 0, 1 } },

  { { 0, 1 }, { 0, 1 }, { 1, 2 } },
  { { 1, 2 }, { 0, 1 }, { 1, 2 } },
  { { 0, 1 }, { 1, 2 }, { 1, 2 } },
  { { 1, 2 }, { 1, 2 }, { 1, 2 } },
};

//------------------------------------------------------------------------------
void vtkIncrementalOctreeNode::SeperateExactlyDuplicatePointsFromNewInsertion(vtkPoints* points,
  vtkIdList* pntIds, const double newPnt[3], vtkIdType* pntIdx, int maxPts, int ptMode)
{
  // the number of points already maintained in this leaf node
  // >= maxPts AND all of them are exactly duplicate with one another
  //           BUT the new point is not a duplicate of them any more

  int i;
  double dupPnt[3];
  double octMin[3];
  double octMid[3];
  double octMax[3];
  double* boxPtr[3] = { nullptr, nullptr, nullptr };
  vtkIncrementalOctreeNode* ocNode = nullptr;
  vtkIncrementalOctreeNode* duplic = this;
  vtkIncrementalOctreeNode* single = this;

  // the coordinate of the duplicate points: note pntIds == this->PointIdSet
  points->GetPoint(pntIds->GetId(0), dupPnt);

  while (duplic == single) // as long as separation has not been achieved
  {
    // update the current (in recursion) node and access the bounding box info
    ocNode = duplic;
    octMid[0] = (ocNode->MinBounds[0] + ocNode->MaxBounds[0]) * 0.5;
    octMid[1] = (ocNode->MinBounds[1] + ocNode->MaxBounds[1]) * 0.5;
    octMid[2] = (ocNode->MinBounds[2] + ocNode->MaxBounds[2]) * 0.5;
    boxPtr[0] = ocNode->MinBounds;
    boxPtr[1] = octMid;
    boxPtr[2] = ocNode->MaxBounds;

    // create eight child nodes
    ocNode->Children = new vtkIncrementalOctreeNode*[8];
    for (i = 0; i < 8; i++)
    {
      // x-bound: axis 0
      octMin[0] = boxPtr[OCTREE_CHILD_BOUNDS_LUT[i][0][0]][0];
      octMax[0] = boxPtr[OCTREE_CHILD_BOUNDS_LUT[i][0][1]][0];

      // y-bound: axis 1
      octMin[1] = boxPtr[OCTREE_CHILD_BOUNDS_LUT[i][1][0]][1];
      octMax[1] = boxPtr[OCTREE_CHILD_BOUNDS_LUT[i][1][1]][1];

      // z-bound: axis 2
      octMin[2] = boxPtr[OCTREE_CHILD_BOUNDS_LUT[i][2][0]][2];
      octMax[2] = boxPtr[OCTREE_CHILD_BOUNDS_LUT[i][2][1]][2];

      ocNode->Children[i] = vtkIncrementalOctreeNode::New();
      ocNode->Children[i]->SetParent(ocNode);
      ocNode->Children[i]->SetBounds(
        octMin[0], octMax[0], octMin[1], octMax[1], octMin[2], octMax[2]);
    }

    // determine the leaf node of the duplicate points & that of the new point
    duplic = ocNode->Children[ocNode->GetChildIndex(dupPnt)];
    single = ocNode->Children[ocNode->GetChildIndex(newPnt)];
  }
  boxPtr[0] = nullptr;
  boxPtr[1] = nullptr;
  boxPtr[2] = nullptr;

  // Now the duplicate points have been separated from the new point //

  // create a vtkIdList object for the new point
  // update the counter and the data bounding box until the root node
  // (including the root node)
  OCTREENODE_INSERTPOINT[ptMode](points, pntIdx, newPnt);
  single->CreatePointIdSet((maxPts >> 2), (maxPts >> 1));
  single->GetPointIdSet()->InsertNextId(*pntIdx);
  single->UpdateCounterAndDataBoundsRecursively(newPnt, 1, 1, nullptr);

  // We just need to reference pntIds while un-registering it from 'this'.
  // This avoids deep-copying point ids from pntIds to duplic's PointIdSet.
  // update the counter and the data bounding box, but until 'this' node
  // (excluding 'this' node)
  duplic->SetPointIdSet(pntIds);
  duplic->UpdateCounterAndDataBoundsRecursively(dupPnt, pntIds->GetNumberOfIds(), 1, this);

  // handle memory
  ocNode = nullptr;
  duplic = nullptr;
  single = nullptr;
}

//------------------------------------------------------------------------------
int vtkIncrementalOctreeNode::CreateChildNodes(vtkPoints* points, vtkIdList* pntIds,
  const double newPnt[3], vtkIdType* pntIdx, int maxPts, int ptMode, int& numberOfNodes)
{
  // There are two scenarios for which this function is invoked.
  //
  // (1) the number of points already maintained in this leaf node
  //     == maxPts AND not all of them are exactly duplicate
  //               AND the new point is not a duplicate of them all
  // (2) the number of points already maintained in this leaf node
  //     >= maxPts AND all of them are exactly duplicate with one another
  //               BUT the new point is not a duplicate of them any more

  // address case (2) first if necessary
  double sample[3];
  points->GetPoint(pntIds->GetId(0), sample);
  if (this->ContainsDuplicatePointsOnly(sample) == 1)
  {
    this->SeperateExactlyDuplicatePointsFromNewInsertion(
      points, pntIds, newPnt, pntIdx, maxPts, ptMode);

    // notify vtkIncrementalOctreeNode::InsertPoint() that pntIds just needs
    // to be unregistered from 'this', but must NOT be destroyed at all.
    return 0;
  }

  // then address case (1) below

  int i;
  int target;
  int dvidId = -1; // index of the sub-dividing octant, if any
  int fullId = -1; // index of the full octant, if any
  int numIds[8] = { 0, 0, 0, 0, 0, 0, 0, 0 };
  double octMin[3];
  double octMax[3];
  double tempPt[3];
  vtkIdType tempId;

  double octMid[3] = { (this->MinBounds[0] + this->MaxBounds[0]) * 0.5,
    (this->MinBounds[1] + this->MaxBounds[1]) * 0.5,
    (this->MinBounds[2] + this->MaxBounds[2]) * 0.5 };
  double* boxPtr[3];
  boxPtr[0] = this->MinBounds;
  boxPtr[1] = octMid;
  boxPtr[2] = this->MaxBounds;

  // create eight child nodes
  this->Children = new vtkIncrementalOctreeNode*[8];
  for (i = 0; i < 8; i++)
  {
    // x-bound: axis 0
    octMin[0] = boxPtr[OCTREE_CHILD_BOUNDS_LUT[i][0][0]][0];
    octMax[0] = boxPtr[OCTREE_CHILD_BOUNDS_LUT[i][0][1]][0];

    // y-bound: axis 1
    octMin[1] = boxPtr[OCTREE_CHILD_BOUNDS_LUT[i][1][0]][1];
    octMax[1] = boxPtr[OCTREE_CHILD_BOUNDS_LUT[i][1][1]][1];

    // z-bound: axis 2
    octMin[2] = boxPtr[OCTREE_CHILD_BOUNDS_LUT[i][2][0]][2];
    octMax[2] = boxPtr[OCTREE_CHILD_BOUNDS_LUT[i][2][1]][2];

    // This call internally sets the cener and default data bounding box, too.
    this->Children[i] = vtkIncrementalOctreeNode::New();
    this->Children[i]->ID = numberOfNodes++;
    this->Children[i]->SetParent(this);
    this->Children[i]->SetBounds(octMin[0], octMax[0], octMin[1], octMax[1], octMin[2], octMax[2]);

    // allocate a list of point-indices (size = 2^n) for index registration
    this->Children[i]->CreatePointIdSet((maxPts >> 2), (maxPts >> 1));
  }
  boxPtr[0] = nullptr;
  boxPtr[1] = nullptr;
  boxPtr[2] = nullptr;

  // distribute the available point-indices to the eight child nodes
  for (i = 0; i < maxPts; i++)
  {
    tempId = pntIds->GetId(i);
    points->GetPoint(tempId, tempPt);
    target = this->GetChildIndex(tempPt);
    this->Children[target]->GetPointIdSet()->InsertNextId(tempId);
    this->Children[target]->UpdateCounterAndDataBounds(tempPt);
    numIds[target]++;
  }

  // locate the full child, just if any
  for (i = 0; i < 8; i++)
  {
    if (numIds[i] == maxPts)
    {
      fullId = i;
      break;
    }
  }

  target = this->GetChildIndex(newPnt);
  if (fullId == target)
  {
    // The fact is that we are going to insert the new point to an already
    // full octant (child node). Thus we need to further divide this child
    // to avoid the overflow problem.
    this->Children[target]->CreateChildNodes(
      points, pntIds, newPnt, pntIdx, maxPts, ptMode, numberOfNodes);
    dvidId = fullId;
  }
  else
  {
    // the initial division is a success
    OCTREENODE_INSERTPOINT[ptMode](points, pntIdx, newPnt);
    this->Children[target]->GetPointIdSet()->InsertNextId(*pntIdx);
    this->Children[target]->UpdateCounterAndDataBoundsRecursively(newPnt, 1, 1, nullptr);

    // NOTE: The counter below might reach the threshold, though we delay the
    // sub-division of this child node until the next point insertion occurs.
    numIds[target]++;
  }

  // Now it is time to reclaim those un-used vtkIdList objects, of which each
  // either is empty or still needs to be deleted due to further division of
  // the child node. This post-deallocation of the un-used vtkIdList objects
  // (of some child nodes) is based on the assumption that retrieving the
  // 'maxPts' points from vtkPoints and the associated 'maxPts' point-indices
  // from vtkIdList is more expensive than reclaiming at most 8 vtkIdList
  // objects at hand.
  for (i = 0; i < 8; i++)
  {
    if (numIds[i] == 0 || i == dvidId)
    {
      this->Children[i]->DeletePointIdSet();
    }
  }

  // notify vtkIncrementalOctreeNode::InsertPoint() to destroy pntIds
  return 1;
}

//------------------------------------------------------------------------------
int vtkIncrementalOctreeNode::InsertPoint(
  vtkPoints* points, const double newPnt[3], int maxPts, vtkIdType* pntId, int ptMode)
{
  VTK_LEGACY_REPLACED_BODY(vtkIncrementalOctreeNode::InsertPoint(vtkPoints * points,
                             const double newPnt[3], int maxPts, vtkIdType* pntId, int ptMode),
    "VTK 9.1",
    vtkIncrementalOctreeNode::InsertPoint(vtkPoints * points, const double newPnt[3], int maxPts,
      vtkIdType* pntId, int ptMode, int& numberOfNodes));
  int numberOfNodes = 0;
  return InsertPoint(points, newPnt, maxPts, pntId, ptMode, numberOfNodes);
}

//------------------------------------------------------------------------------
int vtkIncrementalOctreeNode::InsertPoint(vtkPoints* points, const double newPnt[3], int maxPts,
  vtkIdType* pntId, int ptMode, int& numberOfNodes)
{
  if (this->PointIdSet)
  {
    // there has been at least one point index
    if (this->PointIdSet->GetNumberOfIds() < maxPts ||
      this->ContainsDuplicatePointsOnly(newPnt) == 1)
    {
      // this leaf node is not full or
      // this leaf node is full, but of all exactly duplicate points
      // and the point under check is another duplicate of these points
      OCTREENODE_INSERTPOINT[ptMode](points, pntId, newPnt);
      this->PointIdSet->InsertNextId(*pntId);
      this->UpdateCounterAndDataBoundsRecursively(newPnt, 1, 1, nullptr);
    }
    else
    {
      // overflow: divide this node and delete the list of point-indices.
      // Note that the number of exactly duplicate points might be greater
      // than or equal to maxPts.
      if (this->CreateChildNodes(
            points, this->PointIdSet, newPnt, pntId, maxPts, ptMode, numberOfNodes))
      {
        this->PointIdSet->Delete();
      }
      else
      {
        this->PointIdSet->UnRegister(this);
      }
      this->PointIdSet = nullptr;
    }
  }
  else
  {
    // there has been no any point index registered in this leaf node
    OCTREENODE_INSERTPOINT[ptMode](points, pntId, newPnt);
    this->PointIdSet = vtkIdList::New();
    this->PointIdSet->Allocate((maxPts >> 2), (maxPts >> 1));
    this->PointIdSet->InsertNextId(*pntId);
    this->UpdateCounterAndDataBoundsRecursively(newPnt, 1, 1, nullptr);
  }

  return 1;
}

//------------------------------------------------------------------------------
double vtkIncrementalOctreeNode::GetDistance2ToBoundary(const double point[3], double closest[3],
  int innerOnly, vtkIncrementalOctreeNode* rootNode, int checkData)
{
  // It is mandatory that GetMinDataBounds() and GetMaxDataBounds() be used.
  // Direct access to MinDataBounds and MaxDataBounds might incur problems.
  double* thisMin = nullptr;
  double* thisMax = nullptr;
  double* rootMin = nullptr;
  double* rootMax = nullptr;
  double minDist = VTK_DOUBLE_MAX; // minimum distance to the boundaries
  if (checkData)
  {
    thisMin = this->GetMinDataBounds();
    thisMax = this->GetMaxDataBounds();
    rootMin = rootNode->GetMinDataBounds();
    rootMax = rootNode->GetMaxDataBounds();
  }
  else
  {
    thisMin = this->MinBounds;
    thisMax = this->MaxBounds;
    rootMin = rootNode->GetMinBounds();
    rootMax = rootNode->GetMaxBounds();
  }

  int minFace = 0; // index of the face with min distance to the point
  int beXless = int(point[0] < thisMin[0]);
  int beXmore = int(point[0] > thisMax[0]);
  int beYless = int(point[1] < thisMin[1]);
  int beYmore = int(point[1] > thisMax[1]);
  int beZless = int(point[2] < thisMin[2]);
  int beZmore = int(point[2] > thisMax[2]);
  int withinX = int((!beXless) && (!beXmore));
  int withinY = int((!beYless) && (!beYmore));
  int withinZ = int((!beZless) && (!beZmore));
  int xyzFlag = (withinZ << 2) + (withinY << 1) + withinX;

  switch (xyzFlag)
  {
    case 0: // withinZ = 0; withinY = 0;  withinX = 0
            // closest to a corner

      closest[0] = (beXless ? thisMin[0] : thisMax[0]);
      closest[1] = (beYless ? thisMin[1] : thisMax[1]);
      closest[2] = (beZless ? thisMin[2] : thisMax[2]);
      minDist = vtkMath::Distance2BetweenPoints(point, closest);
      break;

    case 1: // withinZ = 0; withinY = 0; withinX = 1
            // closest to an x-aligned edge

      closest[0] = point[0];
      closest[1] = (beYless ? thisMin[1] : thisMax[1]);
      closest[2] = (beZless ? thisMin[2] : thisMax[2]);
      minDist = vtkMath::Distance2BetweenPoints(point, closest);
      break;

    case 2: // withinZ = 0; withinY = 1; withinX = 0
            // closest to a y-aligned edge

      closest[0] = (beXless ? thisMin[0] : thisMax[0]);
      closest[1] = point[1];
      closest[2] = (beZless ? thisMin[2] : thisMax[2]);
      minDist = vtkMath::Distance2BetweenPoints(point, closest);
      break;

    case 3: // withinZ = 0; withinY = 1; withinX = 1
            // closest to a z-face

      if (beZless)
      {
        minDist = thisMin[2] - point[2];
        closest[2] = thisMin[2];
      }
      else
      {
        minDist = point[2] - thisMax[2];
        closest[2] = thisMax[2];
      }

      minDist *= minDist;
      closest[0] = point[0];
      closest[1] = point[1];
      break;

    case 4: // withinZ = 1; withinY = 0; withinX = 0
            // cloest to a z-aligned edge

      closest[0] = (beXless ? thisMin[0] : thisMax[0]);
      closest[1] = (beYless ? thisMin[1] : thisMax[1]);
      closest[2] = point[2];
      minDist = vtkMath::Distance2BetweenPoints(point, closest);
      break;

    case 5: // withinZ = 1; withinY = 0; withinX = 1
            // closest to a y-face

      if (beYless)
      {
        minDist = thisMin[1] - point[1];
        closest[1] = thisMin[1];
      }
      else
      {
        minDist = point[1] - thisMax[1];
        closest[1] = thisMax[1];
      }

      minDist *= minDist;
      closest[0] = point[0];
      closest[2] = point[2];
      break;

    case 6: // withinZ = 1; withinY = 1; withinX = 0
            // closest to an x-face

      if (beXless)
      {
        minDist = thisMin[0] - point[0];
        closest[0] = thisMin[0];
      }
      else
      {
        minDist = point[0] - thisMax[0];
        closest[0] = thisMax[0];
      }

      minDist *= minDist;
      closest[1] = point[1];
      closest[2] = point[2];
      break;

    case 7: // withinZ = 1; withinY = 1;  withinZ = 1
            // point is inside the box

      if (innerOnly) // check only inner boundaries
      {
        double faceDst;

        faceDst = point[0] - thisMin[0]; // x-min face
        if (thisMin[0] != rootMin[0] && faceDst < minDist)
        {
          minFace = 0;
          minDist = faceDst;
        }

        faceDst = thisMax[0] - point[0]; // x-max face
        if (thisMax[0] != rootMax[0] && faceDst < minDist)
        {
          minFace = 1;
          minDist = faceDst;
        }

        faceDst = point[1] - thisMin[1]; // y-min face
        if (thisMin[1] != rootMin[1] && faceDst < minDist)
        {
          minFace = 2;
          minDist = faceDst;
        }

        faceDst = thisMax[1] - point[1]; // y-max face
        if (thisMax[1] != rootMax[1] && faceDst < minDist)
        {
          minFace = 3;
          minDist = faceDst;
        }

        faceDst = point[2] - thisMin[2]; // z-min face
        if (thisMin[2] != rootMin[2] && faceDst < minDist)
        {
          minFace = 4;
          minDist = faceDst;
        }

        faceDst = thisMax[2] - point[2]; // z-max face
        if (thisMax[2] != rootMax[2] && faceDst < minDist)
        {
          minFace = 5;
          minDist = faceDst;
        }
      }
      else // check all boundaries
      {
        double tmpDist[6];
        tmpDist[0] = point[0] - thisMin[0];
        tmpDist[1] = thisMax[0] - point[0];
        tmpDist[2] = point[1] - thisMin[1];
        tmpDist[3] = thisMax[1] - point[1];
        tmpDist[4] = point[2] - thisMin[2];
        tmpDist[5] = thisMax[2] - point[2];

        for (int i = 0; i < 6; i++)
        {
          if (tmpDist[i] < minDist)
          {
            minFace = i;
            minDist = tmpDist[i];
          }
        }
      }

      // no square operation if no any inner boundary
      if (minDist != VTK_DOUBLE_MAX)
      {
        minDist *= minDist;
      }

      closest[0] = point[0];
      closest[1] = point[1];
      closest[2] = point[2];

      // minFace: the quad with the min distance to the point
      // 0: x-min face  ===>  xyzIndx = 0:  x  and  minFace & 1 = 0:  thisMin
      // 1: x-max face  ===>  xyzIndx = 0:  x  and  minFace & 1 = 1:  thisMax
      // 2: y-min face  ===>  xyzIndx = 1:  y  and  minFace & 1 = 0:  thisMin
      // 3: y-max face  ===>  xyzIndx = 1:  y  and  minFace & 1 = 1:  thisMax
      // 4: z-min face  ===>  xyzIndx = 2:  z  and  minFace & 1 = 0:  thisMin
      // 5: z-max face  ===>  xyzIndx = 2:  z  and  minFace & 1 = 1:  thisMax
      double* pMinMax[2] = { thisMin, thisMax };
      int xyzIndx = (minFace >> 1);
      closest[xyzIndx] = pMinMax[minFace & 1][xyzIndx];
      pMinMax[0] = pMinMax[1] = nullptr;

      break;
  }

  return minDist;
}

//------------------------------------------------------------------------------
double vtkIncrementalOctreeNode::GetDistance2ToBoundary(
  const double point[3], vtkIncrementalOctreeNode* rootNode, int checkData)
{
  double dumbPnt[3];
  return (checkData == 1 && this->GetNumberOfPoints() == 0)
    ? VTK_DOUBLE_MAX
    : this->GetDistance2ToBoundary(point, dumbPnt, 0, rootNode, checkData);
}

//------------------------------------------------------------------------------
double vtkIncrementalOctreeNode::GetDistance2ToBoundary(
  const double point[3], double closest[3], vtkIncrementalOctreeNode* rootNode, int checkData)
{
  return (checkData == 1 && this->GetNumberOfPoints() == 0)
    ? VTK_DOUBLE_MAX
    : this->GetDistance2ToBoundary(point, closest, 0, rootNode, checkData);
}

//------------------------------------------------------------------------------
double vtkIncrementalOctreeNode::GetDistance2ToInnerBoundary(
  const double point[3], vtkIncrementalOctreeNode* rootNode)
{
  double dumbPnt[3];
  return this->GetDistance2ToBoundary(point, dumbPnt, 1, rootNode, 0);
}

//------------------------------------------------------------------------------
void vtkIncrementalOctreeNode::ExportAllPointIdsByInsertion(vtkIdList* idList)
{
  if (this->Children == nullptr)
  {
    for (vtkIdType localId = 0; localId < this->NumberOfPoints; localId++)
    {
      idList->InsertNextId(this->PointIdSet->GetId(localId));
    }
  }
  else
  {
    for (int i = 0; i < 8; i++)
    {
      this->Children[i]->ExportAllPointIdsByInsertion(idList);
    }
  }
}

//------------------------------------------------------------------------------
void vtkIncrementalOctreeNode::ExportAllPointIdsByDirectSet(vtkIdType* pntIdx, vtkIdList* idList)
{
  if (this->Children == nullptr)
  {
    for (vtkIdType localId = 0; localId < this->NumberOfPoints; localId++)
    {
      idList->SetId((*pntIdx), this->PointIdSet->GetId(localId));
      (*pntIdx)++;
    }
  }
  else
  {
    for (int i = 0; i < 8; i++)
    {
      this->Children[i]->ExportAllPointIdsByDirectSet(pntIdx, idList);
    }
  }
}

//------------------------------------------------------------------------------
void vtkIncrementalOctreeNode::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os, indent);
  os << indent << "Index: " << this->ID << endl;
  os << indent << "Parent: " << this->Parent << endl;
  os << indent << "Children: " << this->Children << endl;
  os << indent << "PointIdSet: " << this->PointIdSet << endl;
  os << indent << "NumberOfPoints: " << this->NumberOfPoints << 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;
}

//------------------------------------------------------------------------------
int vtkIncrementalOctreeNode::GetNumberOfLevels() const
{
  if (this->Children)
  {
    int maxLevel = 0;
    for (int i = 0; i < 8; i++)
    {
      maxLevel = std::max(maxLevel, this->Children[i]->GetNumberOfLevels());
    }
    return maxLevel + 1;
  }
  else
  {
    return 1;
  }
}