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

  Program:   Visualization Toolkit
  Module:    vtkSphereTree.cxx

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

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

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

#include "vtkCellData.h"
#include "vtkDataArray.h"
#include "vtkDataSet.h"
#include "vtkDebugLeaks.h"
#include "vtkDoubleArray.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkLine.h"
#include "vtkMath.h"
#include "vtkNew.h"
#include "vtkObjectFactory.h"
#include "vtkPlane.h"
#include "vtkPointData.h"
#include "vtkSMPThreadLocal.h"
#include "vtkSMPTools.h"
#include "vtkSphere.h"
#include "vtkStructuredGrid.h"
#include "vtkUnstructuredGrid.h"

#include <vector>

vtkStandardNewMacro(vtkSphereTree);
vtkCxxSetObjectMacro(vtkSphereTree, DataSet, vtkDataSet);

// Implementation notes:
// Currently only two levels of the sphere tree are being built: the leaf
// spheres (one sphere per cell) and then the next level groupings of the
// leaf spheres. This is done because it is easier to thread, and the
// benefits of additional sphere tree hierarchy diminish quickly in a
// threaded environment. Future work may want to revisit this. In particular,
// huge datasets probably would benefit from more levels.
//
// Further room for improvement: while the leaf spheres are built in
// parallel, the hierarchy is built serially. The hierarchy could also
// be built in parallel.
//
// Note the sphere generation uses Ritter's algorithm. While fast, it can
// overestimate the sphere size by 5-20%. Tighter spheres would improve
// performance.

// Type of sphere tree hierarchy generated
#define VTK_SPHERE_TREE_HIERARCHY_NONE 0
#define VTK_SPHERE_TREE_HIERARCHY_STRUCTURED 1
#define VTK_SPHERE_TREE_HIERARCHY_UNSTRUCTURED 2

// Different types of sphere tree hierarchies can be created. These are
// basically data structures for different types of dataset (structured
// and unstructured)
struct vtkSphereTreeHierarchy
{
  virtual ~vtkSphereTreeHierarchy() = default;
};

struct vtkStructuredHierarchy : public vtkSphereTreeHierarchy
{
  vtkIdType NumCells;
  vtkDoubleArray* H;
  vtkIdType Dims[3];
  int Resolution;

  vtkIdType GridSize;
  vtkIdType GridDims[3];
  double* GridSpheres;

  vtkStructuredHierarchy(vtkIdType numCells, vtkIdType size)
    : NumCells(numCells)
  {
    this->Resolution = 0;
    this->Dims[0] = this->Dims[1] = this->Dims[2] = 0;
    this->GridSize = 0;
    this->GridDims[0] = this->GridDims[1] = this->GridDims[2] = 0;
    this->GridSpheres = nullptr;
    this->H = vtkDoubleArray::New();
    this->H->SetNumberOfComponents(1);
    this->H->SetNumberOfTuples(size);
  }

  ~vtkStructuredHierarchy() override
  {
    this->H->Delete();
    this->H = nullptr;
  }
};

// Currently the unstructured hierarchy is one level deep (to keep it
// simple). In the future a full blown hierarchy could be created. Note that
// there is significant cost to memory allocation/deletion etc. so the
// benefits run out quickly.
struct vtkUnstructuredHierarchy : public vtkSphereTreeHierarchy
{
  vtkIdType NumCells;
  int Dims[3];
  double Bounds[6], Spacing[3];
  vtkIdType GridSize;
  vtkIdType* NumSpheres;
  vtkIdType* Offsets;
  vtkIdType* CellLoc;
  vtkIdType* CellMap;
  double* GridSpheres;

  vtkUnstructuredHierarchy(int dims[3], double bounds[6], double spacing[3], vtkIdType numCells)
    : NumCells(numCells)
    , NumSpheres(nullptr)
    , Offsets(nullptr)
    , CellLoc(nullptr)
    , CellMap(nullptr)
    , GridSpheres(nullptr)
  {
    this->GridSize = static_cast<vtkIdType>(dims[0]) * dims[1] * dims[2];
    for (int i = 0; i < 3; ++i)
    {
      this->Dims[i] = dims[i];
      this->Spacing[i] = spacing[i];
      this->Bounds[2 * i] = bounds[2 * i];
      this->Bounds[2 * i + 1] = bounds[2 * i + 1];
    }

    // Create high-level meta structure that points to grid cells
    this->NumSpheres = new vtkIdType[this->GridSize];
    this->Offsets = new vtkIdType[this->GridSize + 1];
    std::fill_n(this->NumSpheres, this->GridSize, 0);
    this->CellLoc = new vtkIdType[numCells];
    this->CellMap = new vtkIdType[numCells];
  }
  ~vtkUnstructuredHierarchy() override
  {
    delete[] this->NumSpheres;
    this->NumSpheres = nullptr;
    delete[] this->Offsets;
    this->Offsets = nullptr;
    delete[] this->CellLoc;
    this->CellLoc = nullptr;
    delete[] this->CellMap;
    this->CellMap = nullptr;
    delete[] this->GridSpheres;
    this->GridSpheres = nullptr;
  }
};

// Threaded helper functions placed in anonymous namespace
//------------------------------------------------------------------------------
namespace
{

//------------------------------------------------------------------------------
// Compute bounds for each cell in any type of dataset
struct DataSetSpheres
{
  vtkDataSet* DataSet;
  double* Spheres;
  bool ComputeBoundsAndRadius;
  double AverageRadius;
  double Bounds[6];
  vtkSMPThreadLocal<double> Radius;
  vtkSMPThreadLocal<vtkIdType> Count;
  vtkSMPThreadLocal<double> XMin;
  vtkSMPThreadLocal<double> XMax;
  vtkSMPThreadLocal<double> YMin;
  vtkSMPThreadLocal<double> YMax;
  vtkSMPThreadLocal<double> ZMin;
  vtkSMPThreadLocal<double> ZMax;

  DataSetSpheres(vtkDataSet* ds, double* s)
    : DataSet(ds)
    , Spheres(s)
    , ComputeBoundsAndRadius(true)
    , AverageRadius(0.0)
  {
    this->Bounds[0] = this->Bounds[2] = this->Bounds[4] = 0.0;
    this->Bounds[1] = this->Bounds[3] = this->Bounds[5] = 0.0;
  }

  void Initialize()
  {
    double& radius = this->Radius.Local();
    radius = 0.0;

    vtkIdType& count = this->Count.Local();
    count = 0;

    double& xmin = this->XMin.Local();
    xmin = VTK_DOUBLE_MAX;
    double& ymin = this->YMin.Local();
    ymin = VTK_DOUBLE_MAX;
    double& zmin = this->ZMin.Local();
    zmin = VTK_DOUBLE_MAX;

    double& xmax = this->XMax.Local();
    xmax = VTK_DOUBLE_MIN;
    double& ymax = this->YMax.Local();
    ymax = VTK_DOUBLE_MIN;
    double& zmax = this->ZMax.Local();
    zmax = VTK_DOUBLE_MIN;
  }

  void operator()(vtkIdType cellId, vtkIdType endCellId)
  {
    vtkDataSet* ds = this->DataSet;
    double* sphere = this->Spheres + 4 * cellId;
    double r, bounds[6];
    double& radius = this->Radius.Local();
    vtkIdType& count = this->Count.Local();
    double& xmin = this->XMin.Local();
    double& ymin = this->YMin.Local();
    double& zmin = this->ZMin.Local();
    double& xmax = this->XMax.Local();
    double& ymax = this->YMax.Local();
    double& zmax = this->ZMax.Local();

    for (; cellId < endCellId; ++cellId, sphere += 4)
    {
      ds->GetCellBounds(cellId, bounds);
      sphere[0] = (bounds[0] + bounds[1]) / 2.0;
      sphere[1] = (bounds[2] + bounds[3]) / 2.0;
      sphere[2] = (bounds[4] + bounds[5]) / 2.0;
      sphere[3] = sqrt((bounds[1] - sphere[0]) * (bounds[1] - sphere[0]) +
        (bounds[3] - sphere[1]) * (bounds[3] - sphere[1]) +
        (bounds[5] - sphere[2]) * (bounds[5] - sphere[2]));

      if (this->ComputeBoundsAndRadius)
      {
        // Keep a bounds for the dataset
        r = sphere[3];
        xmin = std::min(xmin, (sphere[0] - r));
        xmax = std::max(xmax, (sphere[0] + r));
        ymin = std::min(ymin, (sphere[1] - r));
        ymax = std::max(ymax, (sphere[1] + r));
        zmin = std::min(zmin, (sphere[2] - r));
        zmax = std::max(zmax, (sphere[2] + r));

        // Keep a running average of the radius
        count++;
        radius = radius + (r - radius) / static_cast<double>(count);
      }
    }
  }

  // Compute approximation to the average radius, compute bounds
  void Reduce()
  {
    vtkSMPThreadLocal<double>::iterator iter;
    vtkSMPThreadLocal<double>::iterator end = this->Radius.end();

    double aveRadius = 0.0;
    int numThreads = 0;
    for (iter = this->Radius.begin(); iter != end; ++iter)
    {
      numThreads++;
      aveRadius += *iter;
    }
    if (numThreads < 1)
    {
      this->AverageRadius = 1.0;
    }
    else
    {
      this->AverageRadius = aveRadius / static_cast<double>(numThreads);
    }

    // Reduce bounds from all threads
    double xmin = VTK_DOUBLE_MAX;
    for (iter = this->XMin.begin(); iter != this->XMin.end(); ++iter)
    {
      xmin = (*iter < xmin ? *iter : xmin);
    }
    double ymin = VTK_DOUBLE_MAX;
    for (iter = this->YMin.begin(); iter != this->YMin.end(); ++iter)
    {
      ymin = (*iter < ymin ? *iter : ymin);
    }
    double zmin = VTK_DOUBLE_MAX;
    for (iter = this->ZMin.begin(); iter != this->ZMin.end(); ++iter)
    {
      zmin = (*iter < zmin ? *iter : zmin);
    }

    double xmax = VTK_DOUBLE_MIN;
    for (iter = this->XMax.begin(); iter != this->XMax.end(); ++iter)
    {
      xmax = (*iter > xmax ? *iter : xmax);
    }
    double ymax = VTK_DOUBLE_MIN;
    for (iter = this->YMax.begin(); iter != this->YMax.end(); ++iter)
    {
      ymax = (*iter > ymax ? *iter : ymax);
    }
    double zmax = VTK_DOUBLE_MIN;
    for (iter = this->ZMax.begin(); iter != this->ZMax.end(); ++iter)
    {
      zmax = (*iter > zmax ? *iter : zmax);
    }

    this->Bounds[0] = xmin;
    this->Bounds[1] = xmax;
    this->Bounds[2] = ymin;
    this->Bounds[3] = ymax;
    this->Bounds[4] = zmin;
    this->Bounds[5] = zmax;
  }

  void GetBounds(double bounds[6])
  {
    bounds[0] = this->Bounds[0];
    bounds[1] = this->Bounds[1];
    bounds[2] = this->Bounds[2];
    bounds[3] = this->Bounds[3];
    bounds[4] = this->Bounds[4];
    bounds[5] = this->Bounds[5];
  }

  static void Execute(vtkIdType numCells, vtkDataSet* ds, double* s,
    bool vtkNotUsed(computeBoundsAndRadius), double& aveRadius, double sphereBounds[6])
  {
    if (ds->GetNumberOfCells() > 0 && numCells <= ds->GetNumberOfCells())
    {
      // Dummy call to GetCellBounds to enable its uses in the threaded code
      double dummy[6];
      ds->GetCellBounds(0, dummy);

      DataSetSpheres spheres(ds, s);
      vtkSMPTools::For(0, numCells, spheres);
      aveRadius = spheres.AverageRadius;
      spheres.GetBounds(sphereBounds);
    }
  }
};

//------------------------------------------------------------------------------
// Compute bounds for each cell in an unstructured grid
struct UnstructuredSpheres : public DataSetSpheres
{
  UnstructuredSpheres(vtkUnstructuredGrid* grid, double* s)
    : DataSetSpheres(grid, s)
  {
  }

  void Initialize() { DataSetSpheres::Initialize(); }

  void operator()(vtkIdType cellId, vtkIdType endCellId)
  {
    double* sphere = this->Spheres + 4 * cellId;
    vtkUnstructuredGrid* grid = static_cast<vtkUnstructuredGrid*>(this->DataSet);
    double cellPts[120], *p, r;
    vtkIdType ptNum;
    vtkNew<vtkIdList> cellIds;
    vtkIdType numCellPts;
    double& radius = this->Radius.Local();
    vtkIdType& count = this->Count.Local();
    double& xmin = this->XMin.Local();
    double& ymin = this->YMin.Local();
    double& zmin = this->ZMin.Local();
    double& xmax = this->XMax.Local();
    double& ymax = this->YMax.Local();
    double& zmax = this->ZMax.Local();

    for (; cellId < endCellId; ++cellId, sphere += 4)
    {
      grid->GetCellPoints(cellId, cellIds);
      numCellPts = std::min(cellIds->GetNumberOfIds(), static_cast<vtkIdType>(40));
      p = cellPts;
      for (ptNum = 0; ptNum < numCellPts; ++ptNum, p += 3)
      {
        grid->GetPoint(cellIds->GetId(ptNum), p);
      }
      vtkSphere::ComputeBoundingSphere(cellPts, numCellPts, sphere, nullptr);

      if (this->ComputeBoundsAndRadius)
      {
        // Keep a bounds for the grid
        r = sphere[3];
        xmin = std::min(xmin, (sphere[0] - r));
        xmax = std::max(xmax, (sphere[0] + r));
        ymin = std::min(ymin, (sphere[1] - r));
        ymax = std::max(ymax, (sphere[1] + r));
        zmin = std::min(zmin, (sphere[2] - r));
        zmax = std::max(zmax, (sphere[2] + r));
        // Keep a running average of the radius
        count++;
        radius = radius + (r - radius) / static_cast<double>(count);
      }
    }
  }

  void Reduce() { DataSetSpheres::Reduce(); }

  static void Execute(vtkIdType numCells, vtkUnstructuredGrid* grid, double* s,
    bool vtkNotUsed(computeBoundsAndRadius), double& aveRadius, double sphereBounds[6])
  {
    if (grid->GetNumberOfCells() > 0 && numCells <= grid->GetNumberOfCells())
    {
      // Dummy call to GetCellPoints to enable its uses in the threaded code
      vtkNew<vtkIdList> dummy;
      grid->GetCellPoints(0, dummy);

      UnstructuredSpheres spheres(grid, s);
      vtkSMPTools::For(0, numCells, spheres);
      aveRadius = spheres.AverageRadius;
      spheres.GetBounds(sphereBounds);
    }
  }
};

//------------------------------------------------------------------------------
// Compute bounds for each cell in a structured grid
struct StructuredSpheres : public DataSetSpheres
{
  int Dims[3];
  vtkPoints* Points;

  StructuredSpheres(vtkStructuredGrid* grid, double* s)
    : DataSetSpheres(grid, s)
  {
    grid->GetDimensions(this->Dims);
    this->Points = grid->GetPoints();
  }

  void Initialize() { DataSetSpheres::Initialize(); }

  void operator()(vtkIdType slice, vtkIdType endSlice)
  {
    double *p, cellPts[24];
    vtkIdType cellIds[8], ptId, idx, i, j, jOffset, kOffset, sliceOffset, hint[2];
    hint[0] = 0;
    hint[1] = 6;
    int* dims = this->Dims;
    sliceOffset = static_cast<vtkIdType>(dims[0]) * dims[1];
    vtkPoints* inPts = this->Points;
    double* sphere = this->Spheres + slice * 4 * (dims[0] - 1) * (dims[1] - 1);
    for (; slice < endSlice; ++slice)
    {
      kOffset = slice * sliceOffset;
      for (j = 0; j < (dims[1] - 1); ++j)
      {
        jOffset = j * dims[0];
        for (i = 0; i < (dims[0] - 1); ++i)
        {
          ptId = i + jOffset + kOffset;
          cellIds[0] = ptId;
          cellIds[1] = ptId + 1;
          cellIds[2] = ptId + 1 + dims[0];
          cellIds[3] = ptId + dims[0];
          cellIds[4] = ptId + sliceOffset;
          cellIds[5] = ptId + 1 + sliceOffset;
          cellIds[6] = ptId + 1 + dims[0] + sliceOffset;
          cellIds[7] = ptId + dims[0] + sliceOffset;

          p = cellPts;
          for (idx = 0; idx < 8; ++idx, p += 3)
          {
            inPts->GetPoint(cellIds[idx], p);
          }

          vtkSphere::ComputeBoundingSphere(cellPts, 8, sphere, hint);
          sphere += 4;
        } // i
      }   // j
    }     // slices
  }

  void Reduce() { DataSetSpheres::Reduce(); }

  static void Execute(vtkStructuredGrid* grid, double* s)
  {
    StructuredSpheres spheres(grid, s);
    vtkSMPTools::For(0, spheres.Dims[2] - 1, spheres);
  }
};

//------------------------------------------------------------------------------
// Base class for selection of cells via geometric operations
struct BaseCellSelect
{
  vtkIdType NumberOfCells;
  vtkIdType NumberOfCellsSelected;
  vtkSMPThreadLocal<vtkIdType> NumberSelected;
  unsigned char* Selected;
  double* Spheres;
  double Point[3];

  BaseCellSelect(vtkIdType numCells, unsigned char* select, double* spheres, double p[3])
    : NumberOfCells(numCells)
    , NumberOfCellsSelected(0)
    , Selected(select)
    , Spheres(spheres)
  {
    for (int i = 0; i < 3; ++i)
    {
      this->Point[i] = p[i];
    }
    std::fill_n(this->Selected, numCells, 0);
  }

  void Initialize()
  {
    this->NumberOfCellsSelected = 0;
    vtkIdType& numSelected = this->NumberSelected.Local();
    numSelected = 0;
  }

  void Reduce()
  {
    vtkSMPThreadLocal<vtkIdType>::iterator iter;
    vtkSMPThreadLocal<vtkIdType>::iterator end = this->NumberSelected.end();
    this->NumberOfCellsSelected = 0;
    for (iter = this->NumberSelected.begin(); iter != end; ++iter)
    {
      this->NumberOfCellsSelected += *iter;
    }
  }
};

//------------------------------------------------------------------------------
// Select cells from point based on leaf-level spheres (default)
struct DefaultPointSelect : public BaseCellSelect
{
  DefaultPointSelect(vtkIdType numCells, unsigned char* select, double* spheres, double p[3])
    : BaseCellSelect(numCells, select, spheres, p)
  {
  }

  void Initialize() { BaseCellSelect::Initialize(); }

  void operator()(vtkIdType cellId, vtkIdType endCellId)
  {
    double* sphere = this->Spheres + 4 * cellId;
    double* p = this->Point;
    unsigned char* s = this->Selected + cellId;
    vtkIdType& numSelected = this->NumberSelected.Local();

    for (; cellId < endCellId; ++cellId, sphere += 4, ++s)
    {
      if (vtkMath::Distance2BetweenPoints(sphere, p) <= (sphere[3] * sphere[3]))
      {
        *s = 1;
        ++numSelected;
      }
    } // for cells
  }

  void Reduce() { BaseCellSelect::Reduce(); }
};

// Select cells with point from unstructured hierarchy
struct UnstructuredPointSelect : public DefaultPointSelect
{
  vtkUnstructuredHierarchy* H;

  UnstructuredPointSelect(vtkIdType numCells, unsigned char* select, double* spheres, double p[3],
    vtkUnstructuredHierarchy* h)
    : DefaultPointSelect(numCells, select, spheres, p)
    , H(h)
  {
  }

  void Initialize() { DefaultPointSelect::Initialize(); }

  void operator()(vtkIdType gridId, vtkIdType endGridId)
  {
    double* spheres = this->Spheres;
    double* sph;
    double* gs = this->H->GridSpheres + 4 * gridId;
    double* p = this->Point;
    unsigned char* s = this->Selected;
    const vtkIdType* cellMap = this->H->CellMap;
    const vtkIdType* offsets = this->H->Offsets;
    vtkIdType ii, numSph, cellId;
    vtkIdType& numSelected = this->NumberSelected.Local();

    // Loop over grid buckets. The cell spheres that are located in buckets
    // that intersect are processed further.
    for (; gridId < endGridId; ++gridId, gs += 4)
    {
      if (vtkMath::Distance2BetweenPoints(gs, p) <= (gs[3] * gs[3]))
      {
        numSph = offsets[gridId + 1] - offsets[gridId];
        for (ii = 0; ii < numSph; ++ii)
        {
          cellId = *(cellMap + offsets[gridId] + ii);
          sph = spheres + 4 * cellId;
          if (vtkMath::Distance2BetweenPoints(sph, p) <= (sph[3] * sph[3]))
          {
            s[cellId] = 1;
            ++numSelected;
          }
        } // for cells in bucket
      }   // if bucket sphere intersects point
    }     // for grid buckets
  }

  void Reduce() { DefaultPointSelect::Reduce(); }
};

// Select cells from unstructured hierarchy
struct StructuredPointSelect : public DefaultPointSelect
{
  vtkStructuredHierarchy* H;

  StructuredPointSelect(vtkIdType numCells, unsigned char* select, double* spheres, double p[3],
    vtkStructuredHierarchy* h)
    : DefaultPointSelect(numCells, select, spheres, p)
    , H(h)
  {
  }

  void Initialize() { DefaultPointSelect::Initialize(); }

  void operator()(vtkIdType gridId, vtkIdType endGridId)
  {
    double* p = this->Point;
    unsigned char* s = this->Selected;
    double* spheres = this->Spheres;
    double *sph, *gs = this->H->GridSpheres + 4 * gridId;
    const vtkIdType* gridDims = this->H->GridDims;
    int gridSliceOffset = gridDims[0] * gridDims[1];
    const vtkIdType* dims = this->H->Dims;
    vtkIdType i0, j0, k0, i, j, k, jOffset, kOffset, sliceOffset = dims[0] * dims[1];
    vtkIdType iEnd, jEnd, kEnd, cellId;
    int resolution = this->H->Resolution;
    vtkIdType& numSelected = this->NumberSelected.Local();

    // Loop over grid buckets. The cell spheres that are located in buckets
    // that intersect the point are processed further.
    for (; gridId < endGridId; ++gridId, gs += 4)
    {
      if (vtkMath::Distance2BetweenPoints(gs, p) <= (gs[3] * gs[3]))
      {
        // i-j-k coordinates in grid space
        i0 = (gridId % gridDims[0]) * resolution;
        j0 = ((gridId / gridDims[0]) % gridDims[1]) * resolution;
        k0 = (gridId / gridSliceOffset) * resolution;

        iEnd = ((i0 + resolution) < dims[0] ? i0 + resolution : dims[0]);
        jEnd = ((j0 + resolution) < dims[1] ? j0 + resolution : dims[1]);
        kEnd = ((k0 + resolution) < dims[2] ? k0 + resolution : dims[2]);

        // Now loop over resolution*resolution*resolution block of leaf cells
        for (k = k0; k < kEnd; ++k)
        {
          kOffset = k * sliceOffset;
          for (j = j0; j < jEnd; ++j)
          {
            jOffset = j * dims[0];
            for (i = i0; i < iEnd; ++i)
            {
              cellId = (i + jOffset + kOffset);
              sph = spheres + 4 * cellId;
              if (vtkMath::Distance2BetweenPoints(sph, p) <= (sph[3] * sph[3]))
              {
                s[cellId] = 1; // mark as candidate
                ++numSelected;
              }
            }
          }
        }

      } // if bucket sphere contains point
    }   // for grid buckets
  }

  void Reduce() { DefaultPointSelect::Reduce(); }
};

//------------------------------------------------------------------------------
// Select cells from line based on leaf-level spheres (default)
struct DefaultLineSelect : public BaseCellSelect
{
  double P1[3];

  DefaultLineSelect(
    vtkIdType numCells, unsigned char* select, double* spheres, double p[3], double ray[3])
    : BaseCellSelect(numCells, select, spheres, p)
  {
    for (int i = 0; i < 3; ++i)
    {
      this->P1[i] = this->Point[i] + ray[i];
    }
  }

  void Initialize() { BaseCellSelect::Initialize(); }

  void operator()(vtkIdType cellId, vtkIdType endCellId)
  {
    double* sph = this->Spheres + 4 * cellId;
    double *p0 = this->Point, *p1 = this->P1;
    unsigned char* s = this->Selected + cellId;
    vtkIdType& numSelected = this->NumberSelected.Local();

    for (; cellId < endCellId; ++cellId, sph += 4, ++s)
    {
      if (vtkLine::DistanceToLine(sph, p0, p1) <= (sph[3] * sph[3]))
      {
        *s = 1;
        ++numSelected;
      }
    } // for cells
  }

  void Reduce() { BaseCellSelect::Reduce(); }
};

// Select cells with line from unstructured hierarchy
struct UnstructuredLineSelect : public DefaultLineSelect
{
  vtkUnstructuredHierarchy* H;

  UnstructuredLineSelect(vtkIdType numCells, unsigned char* select, double* spheres,
    vtkUnstructuredHierarchy* h, double o[3], double ray[3])
    : DefaultLineSelect(numCells, select, spheres, o, ray)
    , H(h)
  {
  }

  void Initialize() { DefaultLineSelect::Initialize(); }

  void operator()(vtkIdType gridId, vtkIdType endGridId)
  {
    double* spheres = this->Spheres;
    double* sph;
    double* gs = this->H->GridSpheres + 4 * gridId;
    double *p0 = this->Point, *p1 = this->P1;
    unsigned char* s = this->Selected;
    const vtkIdType* cellMap = this->H->CellMap;
    const vtkIdType* offsets = this->H->Offsets;
    vtkIdType ii, numSph, cellId;
    vtkIdType& numSelected = this->NumberSelected.Local();

    // Loop over grid buckets. The cell spheres that are located in buckets
    // that intersect are processed further.
    for (; gridId < endGridId; ++gridId, gs += 4)
    {
      if (vtkLine::DistanceToLine(gs, p0, p1) <= gs[3])
      {
        numSph = offsets[gridId + 1] - offsets[gridId];
        for (ii = 0; ii < numSph; ++ii)
        {
          cellId = *(cellMap + offsets[gridId] + ii);
          sph = spheres + 4 * cellId;
          if (vtkLine::DistanceToLine(sph, p0, p1) <= (sph[3] * sph[3]))
          {
            s[cellId] = 1;
            ++numSelected;
          }
        } // for cells in bucket
      }   // if bucket sphere intersects line
    }     // for grid buckets
  }

  void Reduce() { DefaultLineSelect::Reduce(); }
};

// Select cells from unstructured hierarchy
struct StructuredLineSelect : public DefaultLineSelect
{
  vtkStructuredHierarchy* H;

  StructuredLineSelect(vtkIdType numCells, unsigned char* select, double* spheres,
    vtkStructuredHierarchy* h, double o[3], double ray[3])
    : DefaultLineSelect(numCells, select, spheres, o, ray)
    , H(h)
  {
  }

  void Initialize() { DefaultLineSelect::Initialize(); }

  void operator()(vtkIdType gridId, vtkIdType endGridId)
  {
    double *p0 = this->Point, *p1 = this->P1;
    unsigned char* s = this->Selected;
    double* spheres = this->Spheres;
    double *sph, *gs = this->H->GridSpheres + 4 * gridId;
    const vtkIdType* gridDims = this->H->GridDims;
    int gridSliceOffset = gridDims[0] * gridDims[1];
    const vtkIdType* dims = this->H->Dims;
    vtkIdType i0, j0, k0, i, j, k, jOffset, kOffset, sliceOffset = dims[0] * dims[1];
    vtkIdType iEnd, jEnd, kEnd, cellId;
    int resolution = this->H->Resolution;
    vtkIdType& numSelected = this->NumberSelected.Local();

    // Loop over grid buckets. The cell spheres that are located in buckets
    // that intersect the line are processed further.
    for (; gridId < endGridId; ++gridId, gs += 4)
    {
      if (vtkLine::DistanceToLine(gs, p0, p1) <= gs[3])
      {
        // i-j-k coordinates in grid space
        i0 = (gridId % gridDims[0]) * resolution;
        j0 = ((gridId / gridDims[0]) % gridDims[1]) * resolution;
        k0 = (gridId / gridSliceOffset) * resolution;

        iEnd = ((i0 + resolution) < dims[0] ? i0 + resolution : dims[0]);
        jEnd = ((j0 + resolution) < dims[1] ? j0 + resolution : dims[1]);
        kEnd = ((k0 + resolution) < dims[2] ? k0 + resolution : dims[2]);

        // Now loop over resolution*resolution*resolution block of leaf cells
        for (k = k0; k < kEnd; ++k)
        {
          kOffset = k * sliceOffset;
          for (j = j0; j < jEnd; ++j)
          {
            jOffset = j * dims[0];
            for (i = i0; i < iEnd; ++i)
            {
              cellId = (i + jOffset + kOffset);
              sph = spheres + 4 * cellId;
              if (vtkLine::DistanceToLine(sph, p0, p1) <= (sph[3] * sph[3]))
              {
                s[cellId] = 1; // mark as candidate
                ++numSelected;
              }
            }
          }
        }

      } // if bucket sphere intersects line
    }   // for grid buckets
  }

  void Reduce() { DefaultLineSelect::Reduce(); }
};

//------------------------------------------------------------------------------
// Select cells from plane based on leaf-level spheres (default)
struct DefaultPlaneSelect : public BaseCellSelect
{
  double Normal[3];

  DefaultPlaneSelect(
    vtkIdType numCells, unsigned char* select, double* spheres, double o[3], double n[3])
    : BaseCellSelect(numCells, select, spheres, o)
  {
    for (int i = 0; i < 3; ++i)
    {
      this->Normal[i] = n[i];
    }
    vtkMath::Normalize(this->Normal);
  }

  void Initialize() { BaseCellSelect::Initialize(); }

  void operator()(vtkIdType cellId, vtkIdType endCellId)
  {
    double* sphere = this->Spheres + 4 * cellId;
    double *o = this->Point, *n = this->Normal;
    unsigned char* s = this->Selected + cellId;
    vtkIdType& numSelected = this->NumberSelected.Local();

    for (; cellId < endCellId; ++cellId, sphere += 4, ++s)
    {
      if (vtkPlane::DistanceToPlane(sphere, n, o) <= sphere[3])
      {
        *s = 1;
        ++numSelected;
      }
    } // for cells
  }

  void Reduce() { BaseCellSelect::Reduce(); }
};

// Select cells with plane from unstructured hierarchy
struct UnstructuredPlaneSelect : public DefaultPlaneSelect
{
  vtkUnstructuredHierarchy* H;

  UnstructuredPlaneSelect(vtkIdType numCells, unsigned char* select, double* spheres,
    vtkUnstructuredHierarchy* h, double o[3], double n[3])
    : DefaultPlaneSelect(numCells, select, spheres, o, n)
    , H(h)
  {
  }

  void Initialize() { DefaultPlaneSelect::Initialize(); }

  void operator()(vtkIdType gridId, vtkIdType endGridId)
  {
    double* spheres = this->Spheres;
    double* sph;
    double* gs = this->H->GridSpheres + 4 * gridId;
    double *o = this->Point, *n = this->Normal;
    unsigned char* s = this->Selected;
    const vtkIdType* cellMap = this->H->CellMap;
    const vtkIdType* offsets = this->H->Offsets;
    vtkIdType ii, numSph, cellId;
    vtkIdType& numSelected = this->NumberSelected.Local();

    // Loop over grid buckets. The cell spheres that are located in buckets
    // that intersect are processed further.
    for (; gridId < endGridId; ++gridId, gs += 4)
    {
      if (vtkPlane::DistanceToPlane(gs, n, o) <= gs[3])
      {
        numSph = offsets[gridId + 1] - offsets[gridId];
        for (ii = 0; ii < numSph; ++ii)
        {
          cellId = *(cellMap + offsets[gridId] + ii);
          sph = spheres + 4 * cellId;
          if (vtkPlane::DistanceToPlane(sph, n, o) <= sph[3])
          {
            s[cellId] = 1;
            ++numSelected;
          }
        } // for cells in bucket
      }   // if bucket sphere intersects plane
    }     // for grid buckets
  }

  void Reduce() { DefaultPlaneSelect::Reduce(); }
};

// Select cells from unstructured hierarchy
struct StructuredPlaneSelect : public DefaultPlaneSelect
{
  vtkStructuredHierarchy* H;

  StructuredPlaneSelect(vtkIdType numCells, unsigned char* select, double* spheres,
    vtkStructuredHierarchy* h, double o[3], double n[3])
    : DefaultPlaneSelect(numCells, select, spheres, o, n)
    , H(h)
  {
  }

  void Initialize() { DefaultPlaneSelect::Initialize(); }

  void operator()(vtkIdType gridId, vtkIdType endGridId)
  {
    double *o = this->Point, *n = this->Normal;
    unsigned char* s = this->Selected;
    double* spheres = this->Spheres;
    double *sph, *gs = this->H->GridSpheres + 4 * gridId;
    const vtkIdType* gridDims = this->H->GridDims;
    int gridSliceOffset = gridDims[0] * gridDims[1];
    const vtkIdType* dims = this->H->Dims;
    vtkIdType i0, j0, k0, i, j, k, jOffset, kOffset, sliceOffset = dims[0] * dims[1];
    vtkIdType iEnd, jEnd, kEnd, cellId;
    int resolution = this->H->Resolution;
    vtkIdType& numSelected = this->NumberSelected.Local();

    // Loop over grid buckets. The cell spheres that are located in buckets
    // that intersect the plane are processed further.
    for (; gridId < endGridId; ++gridId, gs += 4)
    {
      if (vtkPlane::DistanceToPlane(gs, n, o) <= gs[3])
      {
        // i-j-k coordinates in grid space
        i0 = (gridId % gridDims[0]) * resolution;
        j0 = ((gridId / gridDims[0]) % gridDims[1]) * resolution;
        k0 = (gridId / gridSliceOffset) * resolution;

        iEnd = ((i0 + resolution) < dims[0] ? i0 + resolution : dims[0]);
        jEnd = ((j0 + resolution) < dims[1] ? j0 + resolution : dims[1]);
        kEnd = ((k0 + resolution) < dims[2] ? k0 + resolution : dims[2]);

        // Now loop over resolution*resolution*resolution block of leaf cells
        for (k = k0; k < kEnd; ++k)
        {
          kOffset = k * sliceOffset;
          for (j = j0; j < jEnd; ++j)
          {
            jOffset = j * dims[0];
            for (i = i0; i < iEnd; ++i)
            {
              cellId = (i + jOffset + kOffset);
              sph = spheres + 4 * cellId;
              if (vtkPlane::DistanceToPlane(sph, n, o) <= sph[3])
              {
                s[cellId] = 1; // mark as candidate
                ++numSelected;
              }
            }
          }
        }

      } // if bucket sphere intersects plane
    }   // for grid buckets
  }

  void Reduce() { DefaultPlaneSelect::Reduce(); }
};

} // anonymous namespace

//================================Sphere Tree class proper===================
//------------------------------------------------------------------------------
// Construct object.
vtkSphereTree::vtkSphereTree()
{
  this->DataSet = nullptr;
  this->Selected = nullptr;
  this->Resolution = 3;
  this->MaxLevel = 10;
  this->NumberOfLevels = 0;
  this->Tree = nullptr;
  this->Hierarchy = nullptr;
  this->BuildHierarchy = true;
  this->SphereTreeType = VTK_SPHERE_TREE_HIERARCHY_NONE;
  this->AverageRadius = 0.0;
  this->SphereBounds[0] = this->SphereBounds[1] = this->SphereBounds[2] = 0.0;
  this->SphereBounds[3] = this->SphereBounds[4] = this->SphereBounds[5] = 0.0;
}

//------------------------------------------------------------------------------
// Destroy object.
vtkSphereTree::~vtkSphereTree()
{
  this->SetDataSet(nullptr);
  delete[] this->Selected;
  delete this->Hierarchy;
  if (this->Tree)
  {
    this->Tree->Delete();
    this->Tree = nullptr;
  }
}

//================General tree methods========================================
//------------------------------------------------------------------------------
void vtkSphereTree::Build()
{
  if (!this->DataSet)
  {
    return;
  }
  else
  {
    this->Build(this->DataSet);
  }
}

//------------------------------------------------------------------------------
void vtkSphereTree::Build(vtkDataSet* input)
{
  this->SetDataSet(input);

  if (this->Tree != nullptr && this->Hierarchy != nullptr && this->BuildTime > this->MTime &&
    (this->BuildTime > this->DataSet->GetMTime()))
  {
    return;
  }

  this->SphereTreeType = VTK_SPHERE_TREE_HIERARCHY_NONE;
  this->BuildTreeSpheres(input);
  if (this->BuildHierarchy)
  {
    this->BuildTreeHierarchy(input);
  }

  this->BuildTime.Modified();
}

//------------------------------------------------------------------------------
// Compute the sphere tree leafs (i.e., spheres around each cell)
void vtkSphereTree::BuildTreeSpheres(vtkDataSet* input)
{
  // See if anything has to be done
  if (this->Tree != nullptr && this->BuildTime > this->MTime)
  {
    return;
  }
  else if (this->Tree != nullptr)
  {
    this->Tree->Delete();
    delete[] this->Selected;
  }

  // Allocate
  //
  vtkIdType numCells = input->GetNumberOfCells();
  vtkDoubleArray* newScalars = vtkDoubleArray::New();
  newScalars->SetNumberOfComponents(4);
  newScalars->SetNumberOfTuples(input->GetNumberOfCells());
  this->Tree = newScalars;
  this->TreePtr = newScalars->GetPointer(0);

  this->Selected = new unsigned char[numCells];

  if (input->GetDataObjectType() == VTK_STRUCTURED_GRID)
  {
    StructuredSpheres::Execute(vtkStructuredGrid::SafeDownCast(input), this->TreePtr);
  }

  else if (input->GetDataObjectType() == VTK_UNSTRUCTURED_GRID)
  {
    UnstructuredSpheres::Execute(numCells, vtkUnstructuredGrid::SafeDownCast(input), this->TreePtr,
      this->BuildHierarchy, this->AverageRadius, this->SphereBounds);
  }

  else // default algorithm
  {
    DataSetSpheres::Execute(numCells, input, this->TreePtr, this->BuildHierarchy,
      this->AverageRadius, this->SphereBounds);
  }

  this->BuildTime.Modified();
}

//------------------------------------------------------------------------------
void vtkSphereTree::BuildTreeHierarchy(vtkDataSet* input)
{
  if (input->GetDataObjectType() == VTK_STRUCTURED_GRID)
  {
    vtkSphereTree::BuildStructuredHierarchy(vtkStructuredGrid::SafeDownCast(input), this->TreePtr);
  }

  else if (input->GetDataObjectType() == VTK_UNSTRUCTURED_GRID)
  {
    vtkSphereTree::BuildUnstructuredHierarchy(
      vtkUnstructuredGrid::SafeDownCast(input), this->TreePtr);
  }

  else // default hierarchy
  {
    vtkSphereTree::BuildUnstructuredHierarchy(input, this->TreePtr);
  }

  this->BuildTime.Modified();
}

//================Specialized methods for structured grids====================
//------------------------------------------------------------------------------
// From the leaf spheres, build a sphere tree. Use the structure of the grid
// to control how the sphere tree hierarchy is constructed.
void vtkSphereTree::BuildStructuredHierarchy(vtkStructuredGrid* input, double* tree)
{
  this->SphereTreeType = VTK_SPHERE_TREE_HIERARCHY_STRUCTURED;

  // Determine the lay of the land. Note that the code below can build more than
  // the two levels, but for now we clamp to just two levels (the tree leaf
  // spheres plus one level up).
  //  vtkIdType numLevels = this->NumberOfLevels = this->MaxLevel;
  vtkIdType numLevels = this->NumberOfLevels = 2;
  vtkIdType i;
  int lDims[VTK_MAX_SPHERE_TREE_LEVELS][3], size[VTK_MAX_SPHERE_TREE_LEVELS];
  int resolution = this->Resolution;

  // Configure the various levels
  int curLevel = numLevels - 1;
  input->GetDimensions(lDims[curLevel]);
  lDims[curLevel][0] -= 1;
  lDims[curLevel][1] -= 1;
  lDims[curLevel][2] -= 1;
  size[curLevel] = lDims[curLevel][0] * lDims[curLevel][1] * lDims[curLevel][2];
  vtkIdType totalSize = 0;
  for (i = numLevels - 2; i >= 0; --i)
  {
    lDims[i][0] = (lDims[i + 1][0] - 1) / resolution + 1;
    lDims[i][1] = (lDims[i + 1][1] - 1) / resolution + 1;
    lDims[i][2] = (lDims[i + 1][2] - 1) / resolution + 1;
    size[i] = lDims[i][0] * lDims[i][1] * lDims[i][2];
    totalSize += size[i];
  }

  // Allocate space and set up storage.
  delete this->Hierarchy; // cleanup if necessary
  vtkStructuredHierarchy* sH =
    new vtkStructuredHierarchy(input->GetNumberOfCells(), 4 * totalSize + 2);
  this->Hierarchy = sH;

  double *sphere, *spheres[VTK_MAX_SPHERE_TREE_LEVELS];
  spheres[0] = sH->H->GetPointer(0);
  spheres[0][0] = numLevels;
  spheres[0][1] = resolution;
  spheres[0] += 2;
  // As long as numLevels=2; then this is dead code and hence commented out.
  // for (i=1; i < numLevels-1; ++i)
  // {
  //   spheres[i] = spheres[i-1] + 4*size[i-1];
  // }
  spheres[curLevel] = tree;

  // For now, we are going to do something really simple stupid. That is,
  // cull based on blocks of cells one level up from leaf spheres. In the
  // future this will be optimized.
  sH->Dims[0] = lDims[curLevel][0];
  sH->Dims[1] = lDims[curLevel][1];
  sH->Dims[2] = lDims[curLevel][2];
  sH->Resolution = resolution;
  sH->GridSize = size[curLevel - 1];
  sH->GridDims[0] = lDims[curLevel - 1][0];
  sH->GridDims[1] = lDims[curLevel - 1][1];
  sH->GridDims[2] = lDims[curLevel - 1][2];
  sH->GridSpheres = spheres[curLevel - 1];

  // Loop over all levels, from the bottom up, determining sphere tree from level below
  vtkIdType level, j, k, jOffset, kOffset, sliceOffset;
  vtkIdType ii, jj, kk, iStart, iEnd, jStart, jEnd, kStart, kEnd;
  vtkIdType jjOffset, kkOffset, blockSliceOffset, hints[2], numSpheres;
  double* blockSpheres[VTK_MAX_SPHERE_TREE_RESOLUTION * VTK_MAX_SPHERE_TREE_RESOLUTION *
    VTK_MAX_SPHERE_TREE_RESOLUTION];
  hints[0] = 0;

  for (level = numLevels - 2; level >= 0; --level)
  {
    sliceOffset = static_cast<vtkIdType>(lDims[level][0]) * lDims[level][1];
    for (k = 0; k < lDims[level][2]; ++k)
    {
      kOffset = k * sliceOffset;
      kStart = k * resolution;
      kEnd =
        (kStart + resolution < lDims[level + 1][2] ? kStart + resolution : lDims[level + 1][2]);
      for (j = 0; j < lDims[level][1]; ++j)
      {
        jOffset = j * lDims[level][0];
        jStart = j * resolution;
        jEnd =
          (jStart + resolution < lDims[level + 1][1] ? jStart + resolution : lDims[level + 1][1]);
        for (i = 0; i < lDims[level][0]; ++i)
        {
          iStart = i * resolution;
          iEnd =
            (iStart + resolution < lDims[level + 1][0] ? iStart + resolution : lDims[level + 1][0]);
          sphere = spheres[level] + 4 * (i + jOffset + kOffset);
          numSpheres = 0;

          // Now compute bounding sphere for this block of spheres
          hints[1] = (iEnd - iStart) * (jEnd - jStart) * (kEnd - kStart) - 1;

          blockSliceOffset = static_cast<vtkIdType>(lDims[level + 1][0]) * lDims[level + 1][1];
          for (kk = kStart; kk < kEnd; ++kk)
          {
            kkOffset = kk * blockSliceOffset;
            for (jj = jStart; jj < jEnd; ++jj)
            {
              jjOffset = jj * lDims[level + 1][0];
              for (ii = iStart; ii < iEnd; ++ii)
              {
                blockSpheres[numSpheres++] = spheres[level + 1] + 4 * (ii + jjOffset + kkOffset);
              } // for sub-block ii
            }   // for sub-block jj
          }     // for sub-block kk
          vtkSphere::ComputeBoundingSphere(blockSpheres, numSpheres, sphere, hints);
        } // for i
      }   // for j
    }     // for k
  }       // for all levels
}

//================Specialized methods for unstructured grids====================
// Here we create a pointerless binary sphere tree. The order of the spheres
// is implicit with the ordering of the cells. Note that the statistics
// gathered in the previous step are used to organize the grid. The average
// radius controls whether to create lots of spheres or less. Too many
// spheres is wasteful; too few and the computational benefit of the sphere
// tree is reduced.
//
// Based on the average radius and bounds, we'll grid a regular grid
// subdivided n x m x o in the x-y-z directions. We will attempt to
// make the grid buckets cubical. Once the grid is formed, cell
// spheres will be assigned to the grid buckets based on where the
// sphere's center is located. Finally, spheres will be associated
// with each grid bucket (which bound all spheres contained within the
// grid bucket).
void vtkSphereTree::BuildUnstructuredHierarchy(vtkDataSet* input, double* tree)
{
  this->SphereTreeType = VTK_SPHERE_TREE_HIERARCHY_UNSTRUCTURED;

  // Make sure we have something to do.
  vtkIdType numCells = input->GetNumberOfCells();
  if (this->AverageRadius <= 0.0 || numCells <= 0)
  {
    delete this->Hierarchy;
    this->Hierarchy = nullptr;
  }

  // Currently only two levels are being built (see implementation notes).
  this->NumberOfLevels = 2;

  // Compute the grid resolution in the x-y-z directions. Assume that
  // a grid cell should be this->Resolution times bigger than the average
  // radius (in each direction).
  double spacing[3], r = this->AverageRadius, *bds = this->SphereBounds;
  if (bds[1] <= bds[0] || bds[3] <= bds[2] || bds[5] <= bds[4])
  {
    vtkWarningMacro("Invalid bounds, cannot compute tree hierarchy");
    return;
  }
  int dims[3], res = this->Resolution;
  for (int i = 0; i < 3; ++i)
  {
    dims[i] = static_cast<int>((bds[2 * i + 1] - bds[2 * i]) / (res * r));
    dims[i] = (dims[i] < 1 ? 1 : dims[i]);
    spacing[i] = (bds[2 * i + 1] - bds[2 * i]) / dims[i];
  }

  // We are ready to create the hierarchy
  vtkUnstructuredHierarchy* h;
  delete this->Hierarchy; // cleanup if necessary
  this->Hierarchy = h = new vtkUnstructuredHierarchy(dims, bds, spacing, numCells);
  vtkIdType *cellLoc = h->CellLoc, *cellMap = h->CellMap;
  vtkIdType *numSpheres = h->NumSpheres, *offsets = h->Offsets;
  vtkIdType gridSize = h->GridSize;

  // Okay loop over all cell spheres and assign them to the grid cells.
  vtkIdType cellId, i, j, k, ii, idx, sliceOffset = static_cast<vtkIdType>(dims[0]) * dims[1];
  double* sphere = tree;
  for (cellId = 0; cellId < numCells; ++cellId, sphere += 4)
  {
    i = static_cast<int>(dims[0] * (sphere[0] - bds[0]) / (bds[1] - bds[0]));
    j = static_cast<int>(dims[1] * (sphere[1] - bds[2]) / (bds[3] - bds[2]));
    k = static_cast<int>(dims[2] * (sphere[2] - bds[4]) / (bds[5] - bds[4]));
    idx = i + j * dims[0] + k * sliceOffset;
    cellLoc[cellId] = idx;
    numSpheres[idx]++;
  }

  // Compute offsets into linear array. Also remember the max number of spheres
  // in any given bucket (for subsequent memory allocation).
  vtkIdType maxNumSpheres = numSpheres[0];
  offsets[0] = 0;
  for (idx = 1; idx < gridSize; ++idx)
  {
    offsets[idx] = offsets[idx - 1] + numSpheres[idx - 1];
    maxNumSpheres = (numSpheres[idx] > maxNumSpheres ? numSpheres[idx] : maxNumSpheres);
  }
  offsets[gridSize] = numCells;

  // Now associate cells with appropriate grid buckets.
  for (cellId = 0; cellId < numCells; ++cellId)
  {
    idx = cellLoc[cellId];
    *(cellMap + offsets[idx] + numSpheres[idx] - 1) = cellId;
    numSpheres[idx]--; // counting down towards offset
  }

  // Free extra data. What we have left is a grid with cells associated
  // with each bucket.
  delete[] h->NumSpheres;
  h->NumSpheres = nullptr;
  delete[] h->CellLoc;
  h->CellLoc = nullptr;

  // Now it's time to create a sphere per bucket, and adjust the spheres
  // to fit all of the cell spheres contained within it.
  std::vector<double*> tmpSpheres(maxNumSpheres);
  h->GridSpheres = new double[4 * gridSize];
  double* gs = h->GridSpheres;
  vtkIdType nSph;

  for (k = 0; k < dims[2]; ++k)
  {
    for (j = 0; j < dims[1]; ++j)
    {
      for (i = 0; i < dims[0]; ++i)
      {
        idx = i + j * dims[0] + k * sliceOffset;
        nSph = offsets[idx + 1] - offsets[idx];
        for (ii = 0; ii < nSph; ++ii)
        {
          cellId = *(cellMap + offsets[idx] + ii);
          tmpSpheres[ii] = tree + 4 * cellId;
        }
        vtkSphere::ComputeBoundingSphere(tmpSpheres.data(), nSph, gs, nullptr);
        gs += 4;
      } // i
    }   // j
  }     // k
}

//------------------------------------------------------------------------------
// Note that there is a long story behind these crude methods for selecting
// cells based on a sphere tree. Initially there was a complex hierarchy of
// iterators for different dataset types and geometric intersection entities
// (e.g., point, line or plane). However the performance of this approach was
// really poor and the code was excessively complex. To do it right requires
// extensive templating etc. Maybe someday.... In the mean time this approach
// (using a selection mask) is really simple and performs pretty well. It
// also suggests future approaches which use cell locators (and other
// classes) to produce selection masks as well.
const unsigned char* vtkSphereTree::SelectPoint(double x[3], vtkIdType& numSelected)
{
  // Check input
  if (this->DataSet == nullptr)
  {
    return nullptr;
  }

  vtkIdType numCells = this->DataSet->GetNumberOfCells();

  // Specialized for structured grids
  if (this->Hierarchy && this->DataSet->GetDataObjectType() == VTK_STRUCTURED_GRID)
  {
    vtkStructuredHierarchy* h = static_cast<vtkStructuredHierarchy*>(this->Hierarchy);
    vtkIdType gridSize = h->GridSize;
    StructuredPointSelect sPointSelect(numCells, this->Selected, this->TreePtr, x, h);
    vtkSMPTools::For(0, gridSize, sPointSelect);
    numSelected = sPointSelect.NumberOfCellsSelected;
  }

  // Specialized for unstructured grids
  else if (this->Hierarchy && this->DataSet->GetDataObjectType() == VTK_UNSTRUCTURED_GRID)
  {
    vtkUnstructuredHierarchy* h = static_cast<vtkUnstructuredHierarchy*>(this->Hierarchy);
    vtkIdType gridSize = h->GridSize;
    UnstructuredPointSelect uPointSelect(numCells, this->Selected, this->TreePtr, x, h);
    vtkSMPTools::For(0, gridSize, uPointSelect);
    numSelected = uPointSelect.NumberOfCellsSelected;
  }

  // default, process leaf spheres without hierarchy
  else
  {
    DefaultPointSelect defaultPointSelect(numCells, this->Selected, this->TreePtr, x);
    vtkSMPTools::For(0, numCells, defaultPointSelect);
    numSelected = defaultPointSelect.NumberOfCellsSelected;
  }

  return this->Selected;
}

//------------------------------------------------------------------------------
// Create selection mask based on intersection with an infinite line.
const unsigned char* vtkSphereTree::SelectLine(
  double origin[3], double ray[3], vtkIdType& numSelected)
{
  // Check input
  if (this->DataSet == nullptr)
  {
    return nullptr;
  }

  vtkIdType numCells = this->DataSet->GetNumberOfCells();

  // Specialized for structured grids
  if (this->Hierarchy && this->DataSet->GetDataObjectType() == VTK_STRUCTURED_GRID)
  {
    vtkStructuredHierarchy* h = static_cast<vtkStructuredHierarchy*>(this->Hierarchy);
    vtkIdType gridSize = h->GridSize;
    StructuredLineSelect sLineSelect(numCells, this->Selected, this->TreePtr, h, origin, ray);
    vtkSMPTools::For(0, gridSize, sLineSelect);
    numSelected = sLineSelect.NumberOfCellsSelected;
  }

  // Specialized for unstructured grids
  else if (this->Hierarchy && this->DataSet->GetDataObjectType() == VTK_UNSTRUCTURED_GRID)
  {
    vtkUnstructuredHierarchy* h = static_cast<vtkUnstructuredHierarchy*>(this->Hierarchy);
    vtkIdType gridSize = h->GridSize;
    UnstructuredLineSelect uLineSelect(numCells, this->Selected, this->TreePtr, h, origin, ray);
    vtkSMPTools::For(0, gridSize, uLineSelect);
    numSelected = uLineSelect.NumberOfCellsSelected;
  }

  // default, process leaf spheres without hierarchy
  else
  {
    DefaultLineSelect defaultLineSelect(numCells, this->Selected, this->TreePtr, origin, ray);
    vtkSMPTools::For(0, numCells, defaultLineSelect);
    numSelected = defaultLineSelect.NumberOfCellsSelected;
  }

  return this->Selected;
}

//------------------------------------------------------------------------------
// Create selection mask based on intersection with an infinite plane.
const unsigned char* vtkSphereTree::SelectPlane(
  double origin[3], double normal[3], vtkIdType& numSelected)
{
  // Check input
  if (this->DataSet == nullptr)
  {
    return nullptr;
  }

  vtkIdType numCells = this->DataSet->GetNumberOfCells();

  // Specialized for structured grids
  if (this->Hierarchy && this->DataSet->GetDataObjectType() == VTK_STRUCTURED_GRID)
  {
    vtkStructuredHierarchy* h = static_cast<vtkStructuredHierarchy*>(this->Hierarchy);
    vtkIdType gridSize = h->GridSize;
    StructuredPlaneSelect sPlaneSelect(numCells, this->Selected, this->TreePtr, h, origin, normal);
    vtkSMPTools::For(0, gridSize, sPlaneSelect);
    numSelected = sPlaneSelect.NumberOfCellsSelected;
  }

  // Specialized for unstructured grids
  else if (this->Hierarchy && this->DataSet->GetDataObjectType() == VTK_UNSTRUCTURED_GRID)
  {
    vtkUnstructuredHierarchy* h = static_cast<vtkUnstructuredHierarchy*>(this->Hierarchy);
    vtkIdType gridSize = h->GridSize;
    UnstructuredPlaneSelect uPlaneSelect(
      numCells, this->Selected, this->TreePtr, h, origin, normal);
    vtkSMPTools::For(0, gridSize, uPlaneSelect);
    numSelected = uPlaneSelect.NumberOfCellsSelected;
  }

  // default, process leaf spheres without hierarchy
  else
  {
    DefaultPlaneSelect defaultPlaneSelect(numCells, this->Selected, this->TreePtr, origin, normal);
    vtkSMPTools::For(0, numCells, defaultPlaneSelect);
    numSelected = defaultPlaneSelect.NumberOfCellsSelected;
  }

  return this->Selected;
}

//------------------------------------------------------------------------------
// Simply return the leaf spheres
const double* vtkSphereTree::GetCellSpheres()
{
  return this->TreePtr;
}

//------------------------------------------------------------------------------
// The number of levels is this->NumberOfLevels, with
// (this->NumberOfLevels-1) the cell (leaf) spheres, and level 0 the root
// level.
const double* vtkSphereTree::GetTreeSpheres(int level, vtkIdType& numSpheres)
{
  int numLevels = this->NumberOfLevels;

  // Check input for simple cases
  if (level == (numLevels - 1))
  {
    numSpheres = this->DataSet->GetNumberOfCells();
    return this->TreePtr; // just return leaf spheres
  }
  else if (level < 0 || level >= numLevels || this->DataSet == nullptr ||
    this->Hierarchy == nullptr)
  {
    numSpheres = 0;
    return nullptr;
  }

  // Asking for spheres within tree hierarchy
  if (this->SphereTreeType == VTK_SPHERE_TREE_HIERARCHY_STRUCTURED)
  {
    vtkStructuredHierarchy* h = static_cast<vtkStructuredHierarchy*>(this->Hierarchy);
    numSpheres = h->GridSize;
    return h->GridSpheres;
  }
  else if (this->SphereTreeType == VTK_SPHERE_TREE_HIERARCHY_UNSTRUCTURED)
  {
    vtkUnstructuredHierarchy* h = static_cast<vtkUnstructuredHierarchy*>(this->Hierarchy);
    numSpheres = h->GridSize;
    return h->GridSpheres;
  }

  // worst case shouldn't happen
  numSpheres = 0;
  return nullptr;
}

//------------------------------------------------------------------------------
void vtkSphereTree::SelectPoint(double point[3], vtkIdList* cellIds)
{
  vtkIdType numSelected;
  const unsigned char* selected = this->SelectPoint(point, numSelected);
  this->ExtractCellIds(selected, cellIds, numSelected);
}

//------------------------------------------------------------------------------
void vtkSphereTree::SelectLine(double origin[3], double ray[3], vtkIdList* cellIds)
{
  vtkIdType numSelected;
  const unsigned char* selected = this->SelectLine(origin, ray, numSelected);
  this->ExtractCellIds(selected, cellIds, numSelected);
}

//------------------------------------------------------------------------------
void vtkSphereTree::SelectPlane(double origin[3], double normal[3], vtkIdList* cellIds)
{
  vtkIdType numSelected;
  const unsigned char* selected = this->SelectPlane(origin, normal, numSelected);
  this->ExtractCellIds(selected, cellIds, numSelected);
}

//------------------------------------------------------------------------------
void vtkSphereTree::ExtractCellIds(
  const unsigned char* selected, vtkIdList* cellIds, vtkIdType numSelected)
{
  if (numSelected < 1 || selected == nullptr)
  {
    cellIds->Reset();
  }
  else
  {
    const unsigned char* s = selected;
    vtkIdType numCells = this->DataSet->GetNumberOfCells();
    vtkIdType numInserted = 0;
    cellIds->SetNumberOfIds(numSelected);
    for (vtkIdType cellId = 0; cellId < numCells; ++cellId)
    {
      if (*s++ > 0)
      {
        cellIds->SetId(numInserted, cellId);
        numInserted++;
      }
    }
  }
}

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

  os << indent << "Resolution: " << this->Resolution << "\n";
  os << indent << "Number Of Levels: " << this->NumberOfLevels << "\n";
  os << indent << "Maximum Number Of Levels: " << this->MaxLevel << "\n";
  os << indent << "Build Hierarchy: " << (this->BuildHierarchy ? "On\n" : "Off\n");
}