// SPDX-FileCopyrightText: Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
// SPDX-License-Identifier: BSD-3-Clause
#include "vtkStructuredData.h"

#include "vtkCellType.h"
#include "vtkConstantArray.h"
#include "vtkDataSetAttributes.h"
#include "vtkIdList.h"
#include "vtkObjectFactory.h"
#include "vtkPoints.h"
#include "vtkStructuredCellArray.h"
#include "vtkStructuredExtent.h"
#include "vtkStructuredPointArray.h"
#include "vtkUnsignedCharArray.h"

#include <algorithm>
#include <cassert>

VTK_ABI_NAMESPACE_BEGIN
namespace
{
// FIXME
// I don't know if it makes sense to add `REFINEDCELL` in this mask.
// Needs further investigation.
constexpr unsigned char MASKED_CELL_VALUE =
  vtkDataSetAttributes::HIDDENCELL | vtkDataSetAttributes::REFINEDCELL;
} // anonymous namespace

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

// Return the topological dimension of the data (e.g., 0, 1, 2, or 3D).
int vtkStructuredData::GetDataDimension(int dataDescription)
{
  switch (dataDescription)
  {
    case VTK_STRUCTURED_EMPTY:
      return 0; // Should I put -1?
    case VTK_STRUCTURED_SINGLE_POINT:
      return 0;
    case VTK_STRUCTURED_X_LINE:
    case VTK_STRUCTURED_Y_LINE:
    case VTK_STRUCTURED_Z_LINE:
      return 1;
    case VTK_STRUCTURED_XY_PLANE:
    case VTK_STRUCTURED_YZ_PLANE:
    case VTK_STRUCTURED_XZ_PLANE:
      return 2;
    case VTK_STRUCTURED_XYZ_GRID:
      return 3;
    default:
      return -1;
  }
}

//------------------------------------------------------------------------------
int vtkStructuredData::GetDataDimension(int ext[6])
{
  int dataDescription = vtkStructuredData::GetDataDescriptionFromExtent(ext);
  return (vtkStructuredData::GetDataDimension(dataDescription));
}

//------------------------------------------------------------------------------
// Returns the data description given the dimensions (eg. VTK_STRUCTURED_SINGLE_POINT,
// VTK_STRUCTURED_X_LINE, VTK_STRUCTURED_XY_PLANE etc.)
int vtkStructuredData::GetDataDescription(int dims[3])
{
  int tempDims[3];
  // It is essential that dims != tempDims, then alone will SetDimensions()
  // return the correct data description.
  tempDims[0] = dims[0] + 1;
  tempDims[1] = dims[1] + 1;
  tempDims[2] = dims[2] + 1;

  return vtkStructuredData::SetDimensions(dims, tempDims);
}

//------------------------------------------------------------------------------
// Given the extent, returns the data description given the dimensions
// (eg. VTK_STRUCTURED_SINGLE_POINT,VTK_STRUCTURED_X_LINE, VTK_STRUCTURED_XY_PLANE etc.)
int vtkStructuredData::GetDataDescriptionFromExtent(int ext[6])
{
  int dims[3];
  vtkStructuredExtent::GetDimensions(ext, dims);
  return (vtkStructuredData::GetDataDescription(dims));
}

//------------------------------------------------------------------------------
// Specify the dimensions of a regular, rectangular dataset. The input is
// the new dimensions (inDim) and the current dimensions (dim). The function
// returns the dimension of the dataset (0-3D). If the dimensions are
// improperly specified a -1 is returned. If the dimensions are unchanged, a
// value of 100 is returned.
int vtkStructuredData::SetDimensions(VTK_FUTURE_CONST int inDim[3], int dim[3])
{
  int dataDim, i;
  int dataDescription = vtkStructuredData::VTK_STRUCTURED_UNCHANGED;

  if (inDim[0] != dim[0] || inDim[1] != dim[1] || inDim[2] != dim[2])
  {
    for (dataDim = 0, i = 0; i < 3; i++)
    {
      dim[i] = inDim[i];
      if (inDim[i] > 1)
      {
        dataDim++;
      }
    }

    if (inDim[0] < 1 || inDim[1] < 1 || inDim[2] < 1)
    {
      return VTK_STRUCTURED_EMPTY;
    }

    if (dataDim == 3)
    {
      dataDescription = VTK_STRUCTURED_XYZ_GRID;
    }
    else if (dataDim == 2)
    {
      if (inDim[0] == 1)
      {
        dataDescription = VTK_STRUCTURED_YZ_PLANE;
      }
      else if (inDim[1] == 1)
      {
        dataDescription = VTK_STRUCTURED_XZ_PLANE;
      }
      else
      {
        dataDescription = VTK_STRUCTURED_XY_PLANE;
      }
    }
    else if (dataDim == 1)
    {
      if (inDim[0] != 1)
      {
        dataDescription = VTK_STRUCTURED_X_LINE;
      }
      else if (inDim[1] != 1)
      {
        dataDescription = VTK_STRUCTURED_Y_LINE;
      }
      else
      {
        dataDescription = VTK_STRUCTURED_Z_LINE;
      }
    }
    else
    {
      dataDescription = VTK_STRUCTURED_SINGLE_POINT;
    }
  }

  return dataDescription;
}

//------------------------------------------------------------------------------
// Specify the extent of a regular, rectangular dataset. The input is
// the new extent (inExt) and the current extent (ext). The function
// returns the dimension of the dataset (0-3D). If the extents are
// improperly specified a -1 is returned. If the dimensions are unchanged, a
// value of 100 is returned.
int vtkStructuredData::SetExtent(VTK_FUTURE_CONST int inExt[6], int ext[6])
{
  int dataDim, i;
  int dataDescription;

  if (inExt[0] == ext[0] && inExt[1] == ext[1] && inExt[2] == ext[2] && inExt[3] == ext[3] &&
    inExt[4] == ext[4] && inExt[5] == ext[5])
  {
    return vtkStructuredData::VTK_STRUCTURED_UNCHANGED;
  }

  dataDim = 0;
  for (i = 0; i < 3; ++i)
  {
    ext[i * 2] = inExt[i * 2];
    ext[i * 2 + 1] = inExt[i * 2 + 1];
    if (inExt[i * 2] < inExt[i * 2 + 1])
    {
      dataDim++;
    }
  }

  if (inExt[0] > inExt[1] || inExt[2] > inExt[3] || inExt[4] > inExt[5])
  {
    return VTK_STRUCTURED_EMPTY;
  }

  if (dataDim == 3)
  {
    dataDescription = VTK_STRUCTURED_XYZ_GRID;
  }
  else if (dataDim == 2)
  {
    if (inExt[0] == inExt[1])
    {
      dataDescription = VTK_STRUCTURED_YZ_PLANE;
    }
    else if (inExt[2] == inExt[3])
    {
      dataDescription = VTK_STRUCTURED_XZ_PLANE;
    }
    else
    {
      dataDescription = VTK_STRUCTURED_XY_PLANE;
    }
  }
  else if (dataDim == 1)
  {
    if (inExt[0] < inExt[1])
    {
      dataDescription = VTK_STRUCTURED_X_LINE;
    }
    else if (inExt[2] < inExt[3])
    {
      dataDescription = VTK_STRUCTURED_Y_LINE;
    }
    else
    {
      dataDescription = VTK_STRUCTURED_Z_LINE;
    }
  }
  else
  {
    dataDescription = VTK_STRUCTURED_SINGLE_POINT;
  }

  return dataDescription;
}

//------------------------------------------------------------------------------
void vtkStructuredData::ComputeCellStructuredMinMaxCoords(
  vtkIdType cellId, const int dim[3], int ijkMin[3], int ijkMax[3], int dataDescription)
{
  switch (dataDescription)
  {
    case VTK_STRUCTURED_EMPTY:
      ijkMin[0] = ijkMin[1] = ijkMin[2] = 0;
      ijkMax[0] = ijkMax[1] = ijkMax[2] = 0;
      return;

    case VTK_STRUCTURED_SINGLE_POINT:
      ijkMin[0] = ijkMin[1] = ijkMin[2] = 0;
      ijkMax[0] = ijkMax[1] = ijkMax[2] = 0;
      break;

    case VTK_STRUCTURED_X_LINE:
      ijkMin[0] = cellId;
      ijkMin[1] = 0;
      ijkMin[2] = 0;
      ijkMax[0] = ijkMin[0] + 1;
      ijkMax[1] = 0;
      ijkMax[2] = 0;
      break;

    case VTK_STRUCTURED_Y_LINE:
      ijkMin[0] = 0;
      ijkMin[1] = cellId;
      ijkMin[2] = 0;
      ijkMax[0] = 0;
      ijkMax[1] = ijkMin[1] + 1;
      ijkMax[2] = 0;
      break;

    case VTK_STRUCTURED_Z_LINE:
      ijkMin[0] = 0;
      ijkMax[1] = 0;
      ijkMin[2] = cellId;
      ijkMax[0] = 0;
      ijkMin[1] = 0;
      ijkMax[2] = ijkMin[2] + 1;
      break;

    case VTK_STRUCTURED_XY_PLANE:
    {
      const auto div = std::div(cellId, (vtkIdType)(dim[0] - 1));
      ijkMin[0] = div.rem;
      ijkMin[1] = div.quot;
      ijkMin[2] = 0;
      ijkMax[0] = ijkMin[0] + 1;
      ijkMax[1] = ijkMin[1] + 1;
      ijkMax[2] = 0;
      break;
    }

    case VTK_STRUCTURED_YZ_PLANE:
    {
      const auto div = std::div(cellId, (vtkIdType)(dim[1] - 1));
      ijkMin[0] = 0;
      ijkMin[1] = div.rem;
      ijkMin[2] = div.quot;
      ijkMax[0] = 0;
      ijkMax[1] = ijkMin[1] + 1;
      ijkMax[2] = ijkMin[2] + 1;
      break;
    }

    case VTK_STRUCTURED_XZ_PLANE:
    {
      const auto div = std::div(cellId, (vtkIdType)(dim[0] - 1));
      ijkMin[0] = div.rem;
      ijkMin[1] = 0;
      ijkMin[2] = div.quot;
      ijkMax[0] = ijkMin[0] + 1;
      ijkMax[1] = 0;
      ijkMax[2] = ijkMin[2] + 1;
      break;
    }

    case VTK_STRUCTURED_XYZ_GRID:
    {
      const auto div1 = std::div(cellId, (vtkIdType)(dim[0] - 1));
      const auto div2 = std::div(div1.quot, (vtkIdType)(dim[1] - 1));
      ijkMin[0] = div1.rem;
      ijkMin[1] = div2.rem;
      ijkMin[2] = div2.quot;
      ijkMax[0] = ijkMin[0] + 1;
      ijkMax[1] = ijkMin[1] + 1;
      ijkMax[2] = ijkMin[2] + 1;
      break;
    }

    default:
      vtkErrorWithObjectMacro(nullptr, "Invalid DataDescription.");
      return;
  }
}

//------------------------------------------------------------------------------
vtkSmartPointer<vtkStructuredCellArray> vtkStructuredData::GetCellArray(
  int extent[6], bool usePixelVoxelOrientation)
{
  auto implicitCellArray = vtkSmartPointer<vtkStructuredCellArray>::New();
  implicitCellArray->SetData(extent, usePixelVoxelOrientation);
  return implicitCellArray;
}

//------------------------------------------------------------------------------
vtkSmartPointer<vtkConstantArray<unsigned char>> vtkStructuredData::GetCellTypes(
  int extent[6], bool usePixelVoxelOrientation)
{
  const int dataDescription = vtkStructuredData::GetDataDescriptionFromExtent(extent);
  const int dimension = vtkStructuredData::GetDataDimension(dataDescription);
  // NOLINTNEXTLINE(readability-avoid-nested-conditional-operator)
  const int cellType = dimension == 3 ? (usePixelVoxelOrientation ? VTK_VOXEL : VTK_HEXAHEDRON)
    // NOLINTNEXTLINE(readability-avoid-nested-conditional-operator)
    : dimension == 2 ? (usePixelVoxelOrientation ? VTK_PIXEL : VTK_QUAD)
    // NOLINTNEXTLINE(readability-avoid-nested-conditional-operator)
    : dimension == 1 ? VTK_LINE
    // NOLINTNEXTLINE(readability-avoid-nested-conditional-operator)
    : dimension == 0 && dataDescription == VTK_STRUCTURED_SINGLE_POINT ? VTK_VERTEX
                                                                       : VTK_EMPTY_CELL;
  auto cellTypesArray = vtkSmartPointer<vtkConstantArray<unsigned char>>::New();
  cellTypesArray->ConstructBackend(cellType);
  cellTypesArray->SetNumberOfComponents(1);
  cellTypesArray->SetNumberOfTuples(vtkStructuredData::GetNumberOfCells(extent));
  return cellTypesArray;
}

//------------------------------------------------------------------------------
vtkSmartPointer<vtkConstantArray<int>> vtkStructuredData::GetCellTypesArray(
  int extent[6], bool usePixelVoxelOrientation)
{
  auto resultUC = vtkStructuredData::GetCellTypes(extent, usePixelVoxelOrientation);
  auto resultInt = vtkSmartPointer<vtkConstantArray<int>>::New();
  resultInt->ConstructBackend(static_cast<int>(resultUC->GetBackend()->Value));
  resultInt->SetNumberOfComponents(1);
  resultInt->SetNumberOfTuples(resultUC->GetNumberOfTuples());
  return resultInt;
}

//------------------------------------------------------------------------------
// Get the points defining a cell. (See vtkDataSet for more info.)
void vtkStructuredData::GetCellPoints(
  vtkIdType cellId, vtkIdList* ptIds, int dataDescription, int dim[3])
{
  int loc[3];
  vtkIdType idx, npts;
  int iMin, iMax, jMin, jMax, kMin, kMax;
  vtkIdType d01 = static_cast<vtkIdType>(dim[0]) * dim[1];

  ptIds->Reset();
  iMin = iMax = jMin = jMax = kMin = kMax = 0;

  switch (dataDescription)
  {
    case VTK_STRUCTURED_EMPTY:
      return;

    case VTK_STRUCTURED_SINGLE_POINT: // cellId can only be = 0
      break;

    case VTK_STRUCTURED_X_LINE:
      iMin = cellId;
      iMax = cellId + 1;
      break;

    case VTK_STRUCTURED_Y_LINE:
      jMin = cellId;
      jMax = cellId + 1;
      break;

    case VTK_STRUCTURED_Z_LINE:
      kMin = cellId;
      kMax = cellId + 1;
      break;

    case VTK_STRUCTURED_XY_PLANE:
      iMin = cellId % (dim[0] - 1);
      iMax = iMin + 1;
      jMin = cellId / (dim[0] - 1);
      jMax = jMin + 1;
      break;

    case VTK_STRUCTURED_YZ_PLANE:
      jMin = cellId % (dim[1] - 1);
      jMax = jMin + 1;
      kMin = cellId / (dim[1] - 1);
      kMax = kMin + 1;
      break;

    case VTK_STRUCTURED_XZ_PLANE:
      iMin = cellId % (dim[0] - 1);
      iMax = iMin + 1;
      kMin = cellId / (dim[0] - 1);
      kMax = kMin + 1;
      break;

    case VTK_STRUCTURED_XYZ_GRID:
      iMin = cellId % (dim[0] - 1);
      iMax = iMin + 1;
      jMin = (cellId / (dim[0] - 1)) % (dim[1] - 1);
      jMax = jMin + 1;
      kMin = cellId / (static_cast<vtkIdType>(dim[0] - 1) * (dim[1] - 1));
      kMax = kMin + 1;
      break;
    default:
      assert("check: impossible case." && 0); // reaching this line is a bug.
      break;
  }

  // Extract point ids
  for (npts = 0, loc[2] = kMin; loc[2] <= kMax; loc[2]++)
  {
    for (loc[1] = jMin; loc[1] <= jMax; loc[1]++)
    {
      for (loc[0] = iMin; loc[0] <= iMax; loc[0]++)
      {
        idx = loc[0] + loc[1] * static_cast<vtkIdType>(dim[0]) + loc[2] * d01;
        ptIds->InsertId(npts++, idx);
      }
    }
  }
}

//------------------------------------------------------------------------------
vtkSmartPointer<vtkPoints> vtkStructuredData::GetPoints(vtkDataArray* xCoords,
  vtkDataArray* yCoords, vtkDataArray* zCoords, int extent[6], double dirMatrix[9])
{
  auto points = vtkSmartPointer<vtkPoints>::New();
  const auto implicitPointArray = vtk::CreateStructuredPointArray<double>(xCoords, yCoords, zCoords,
    extent, vtkStructuredData::GetDataDescriptionFromExtent(extent), dirMatrix);
  points->SetData(implicitPointArray);
  return points;
}

//------------------------------------------------------------------------------
// Get the cells using a point. (See vtkDataSet for more info.)
void vtkStructuredData::GetPointCells(
  vtkIdType ptId, vtkIdList* cellIds, VTK_FUTURE_CONST int dim[3])
{
  vtkIdType cellDim[3];
  int ptLoc[3], cellLoc[3];
  int i, j;
  vtkIdType cellId;
  static const int offset[8][3] = { { -1, 0, 0 }, { -1, -1, 0 }, { -1, -1, -1 }, { -1, 0, -1 },
    { 0, 0, 0 }, { 0, -1, 0 }, { 0, -1, -1 }, { 0, 0, -1 } };

  for (i = 0; i < 3; i++)
  {
    cellDim[i] = dim[i] - 1;
    if (cellDim[i] == 0)
    {
      cellDim[i] = 1;
    }
  }

  //  Get the location of the point
  //
  ptLoc[0] = ptId % dim[0];
  ptLoc[1] = (ptId / dim[0]) % dim[1];
  ptLoc[2] = ptId / (static_cast<vtkIdType>(dim[0]) * dim[1]);

  //  From the point location, compute the cell locations.  There are at
  //  most eight possible.
  //
  cellIds->Reset();

  for (j = 0; j < 8; j++)
  {
    for (i = 0; i < 3; i++)
    {
      cellLoc[i] = ptLoc[i] + offset[j][i];
      if (cellLoc[i] < 0 || cellLoc[i] >= cellDim[i])
      {
        break;
      }
    }
    if (i >= 3) // add cell
    {
      cellId = cellLoc[0] + cellLoc[1] * cellDim[0] + cellLoc[2] * cellDim[0] * cellDim[1];
      cellIds->InsertNextId(cellId);
    }
  }
}

//------------------------------------------------------------------------------
bool vtkStructuredData::IsPointVisible(vtkIdType pointId, vtkUnsignedCharArray* ghosts)
{
  return !(ghosts && (ghosts->GetValue(pointId) & vtkDataSetAttributes::HIDDENPOINT));
}

//------------------------------------------------------------------------------
bool vtkStructuredData::IsCellVisible(vtkIdType cellId, VTK_FUTURE_CONST int dimensions[3],
  int dataDescription, vtkUnsignedCharArray* cellGhostArray, vtkUnsignedCharArray* pointGhostArray)
{
  if (cellGhostArray && (cellGhostArray->GetValue(cellId) & MASKED_CELL_VALUE))
  {
    return false;
  }

  if (!pointGhostArray)
  {
    return dataDescription != VTK_STRUCTURED_EMPTY;
  }

  int numIds = 0;
  vtkIdType ptIds[8];
  int iMin, iMax, jMin, jMax, kMin, kMax;
  vtkIdType d01 = dimensions[0] * dimensions[1];
  iMin = iMax = jMin = jMax = kMin = kMax = 0;

  switch (dataDescription)
  {
    case VTK_STRUCTURED_SINGLE_POINT: // cellId can only be = 0
      numIds = 1;
      ptIds[0] = iMin + jMin * dimensions[0] + kMin * d01;
      break;

    case VTK_STRUCTURED_X_LINE:
      iMin = cellId;
      iMax = cellId + 1;
      numIds = 2;
      ptIds[0] = iMin + jMin * dimensions[0] + kMin * d01;
      ptIds[1] = iMax + jMin * dimensions[0] + kMin * d01;
      break;

    case VTK_STRUCTURED_Y_LINE:
      jMin = cellId;
      jMax = cellId + 1;
      numIds = 2;
      ptIds[0] = iMin + jMin * dimensions[0] + kMin * d01;
      ptIds[1] = iMin + jMax * dimensions[0] + kMin * d01;
      break;

    case VTK_STRUCTURED_Z_LINE:
      kMin = cellId;
      kMax = cellId + 1;
      numIds = 2;
      ptIds[0] = iMin + jMin * dimensions[0] + kMin * d01;
      ptIds[1] = iMin + jMin * dimensions[0] + kMax * d01;
      break;

    case VTK_STRUCTURED_XY_PLANE:
      iMin = cellId % (dimensions[0] - 1);
      iMax = iMin + 1;
      jMin = cellId / (dimensions[0] - 1);
      jMax = jMin + 1;
      numIds = 4;
      ptIds[0] = iMin + jMin * dimensions[0] + kMin * d01;
      ptIds[1] = iMax + jMin * dimensions[0] + kMin * d01;
      ptIds[2] = iMax + jMax * dimensions[0] + kMin * d01;
      ptIds[3] = iMin + jMax * dimensions[0] + kMin * d01;
      break;

    case VTK_STRUCTURED_YZ_PLANE:
      jMin = cellId % (dimensions[1] - 1);
      jMax = jMin + 1;
      kMin = cellId / (dimensions[1] - 1);
      kMax = kMin + 1;
      numIds = 4;
      ptIds[0] = iMin + jMin * dimensions[0] + kMin * d01;
      ptIds[1] = iMin + jMax * dimensions[0] + kMin * d01;
      ptIds[2] = iMin + jMax * dimensions[0] + kMax * d01;
      ptIds[3] = iMin + jMin * dimensions[0] + kMax * d01;
      break;

    case VTK_STRUCTURED_XZ_PLANE:
      iMin = cellId % (dimensions[0] - 1);
      iMax = iMin + 1;
      kMin = cellId / (dimensions[0] - 1);
      kMax = kMin + 1;
      numIds = 4;
      ptIds[0] = iMin + jMin * dimensions[0] + kMin * d01;
      ptIds[1] = iMax + jMin * dimensions[0] + kMin * d01;
      ptIds[2] = iMax + jMin * dimensions[0] + kMax * d01;
      ptIds[3] = iMin + jMin * dimensions[0] + kMax * d01;
      break;

    case VTK_STRUCTURED_XYZ_GRID:
      iMin = cellId % (dimensions[0] - 1);
      iMax = iMin + 1;
      jMin = (cellId / (dimensions[0] - 1)) % (dimensions[1] - 1);
      jMax = jMin + 1;
      kMin = cellId / ((dimensions[0] - 1) * (dimensions[1] - 1));
      kMax = kMin + 1;
      numIds = 8;
      ptIds[0] = iMin + jMin * dimensions[0] + kMin * d01;
      ptIds[1] = iMax + jMin * dimensions[0] + kMin * d01;
      ptIds[2] = iMax + jMax * dimensions[0] + kMin * d01;
      ptIds[3] = iMin + jMax * dimensions[0] + kMin * d01;
      ptIds[4] = iMin + jMin * dimensions[0] + kMax * d01;
      ptIds[5] = iMax + jMin * dimensions[0] + kMax * d01;
      ptIds[6] = iMax + jMax * dimensions[0] + kMax * d01;
      ptIds[7] = iMin + jMax * dimensions[0] + kMax * d01;
      break;

    default:
      return false;
  }

  for (int i = 0; i < numIds; i++)
  {
    if (!vtkStructuredData::IsPointVisible(ptIds[i], pointGhostArray))
    {
      return false;
    }
  }

  return true;
}

//------------------------------------------------------------------------------
void vtkStructuredData::GetCellNeighbors(
  vtkIdType cellId, vtkIdList* ptIds, vtkIdList* cellIds, int dim[3], int seedLoc[3])
{
  vtkIdType numPts = ptIds->GetNumberOfIds();
  cellIds->Reset();

  // Start by finding the "space" of the points in i-j-k space.
  // The points define a point, line, or plane in topological space,
  // which results in degrees of freedom in three, two or one direction.
  // The numbers of DOF determines which neighbors to select.

  // Start by finding a seed point
  vtkIdType id0 = seedLoc[0] + seedLoc[1] * dim[0] + seedLoc[2] * dim[0] * dim[1];

  // This defines the space around the seed
  int offset[8][3];
  offset[0][0] = -1;
  offset[0][1] = -1;
  offset[0][2] = -1;
  offset[1][0] = 0;
  offset[1][1] = -1;
  offset[1][2] = -1;
  offset[2][0] = -1;
  offset[2][1] = 0;
  offset[2][2] = -1;
  offset[3][0] = 0;
  offset[3][1] = 0;
  offset[3][2] = -1;
  offset[4][0] = -1;
  offset[4][1] = -1;
  offset[4][2] = 0;
  offset[5][0] = 0;
  offset[5][1] = -1;
  offset[5][2] = 0;
  offset[6][0] = -1;
  offset[6][1] = 0;
  offset[6][2] = 0;
  offset[7][0] = 0;
  offset[7][1] = 0;
  offset[7][2] = 0;

  // For the rest of the points, trim the seed region
  // This is essentially an intersection of edge neighbors.
  for (vtkIdType i = 0; i < numPts; i++)
  {
    //  Get the location of the point
    vtkIdType id = ptIds->GetId(i);
    if (id - 1 == id0)
    {
      offset[0][0] = -10;
      offset[2][0] = -10;
      offset[4][0] = -10;
      offset[6][0] = -10;
    }
    else if (id + 1 == id0)
    {
      offset[1][0] = -10;
      offset[3][0] = -10;
      offset[5][0] = -10;
      offset[7][0] = -10;
    }
    else if (id - dim[0] == id0)
    {
      offset[0][1] = -10;
      offset[1][1] = -10;
      offset[4][1] = -10;
      offset[5][1] = -10;
    }
    else if (id + dim[0] == id0)
    {
      offset[2][1] = -10;
      offset[3][1] = -10;
      offset[6][1] = -10;
      offset[7][1] = -10;
    }
    else if (id - dim[0] * dim[1] == id0)
    {
      offset[0][2] = -10;
      offset[1][2] = -10;
      offset[2][2] = -10;
      offset[3][2] = -10;
    }
    else if (id + dim[0] * dim[1] == id0)
    {
      offset[4][2] = -10;
      offset[5][2] = -10;
      offset[6][2] = -10;
      offset[7][2] = -10;
    }
  }

  // Load the non-trimmed cells
  int cellDim[3];
  for (int i = 0; i < 3; i++)
  {
    cellDim[i] = std::max(dim[i] - 1, 1);
  }

  for (int j = 0; j < 8; j++)
  {
    int i = 0, cellLoc[3];
    for (; i < 3; i++)
    {
      if (offset[j][i] != -10)
      {
        cellLoc[i] = seedLoc[i] + offset[j][i];
        if (cellLoc[i] < 0 || cellLoc[i] >= cellDim[i])
        {
          break;
        }
      }
      else
      {
        break;
      }
    }
    if (i >= 3) // add cell
    {
      vtkIdType id = cellLoc[0] + cellLoc[1] * cellDim[0] + cellLoc[2] * cellDim[0] * cellDim[1];
      if (id != cellId)
      {
        cellIds->InsertNextId(id);
      }
    }
  }
}

//------------------------------------------------------------------------------
void vtkStructuredData::GetCellNeighbors(
  vtkIdType cellId, vtkIdList* ptIds, vtkIdList* cellIds, int dim[3])
{
  int j, seedLoc[3], ptLoc[3], cellLoc[3];
  vtkIdType cellDim[3];
  int offset[8][3];
  vtkIdType numPts = ptIds->GetNumberOfIds(), id, i;

  cellIds->Reset();

  // Start by finding the "space" of the points in i-j-k space.
  // The points define a point, line, or plane in topological space,
  // which results in degrees of freedom in three, two or one direction.
  // The numbers of DOF determines which neighbors to select.

  // Start by finding a seed point
  id = ptIds->GetId(0);
  seedLoc[0] = id % dim[0];
  seedLoc[1] = (id / dim[0]) % dim[1];
  seedLoc[2] = id / (dim[0] * dim[1]);

  // This defines the space around the seed
  offset[0][0] = -1;
  offset[0][1] = -1;
  offset[0][2] = -1;
  offset[1][0] = 0;
  offset[1][1] = -1;
  offset[1][2] = -1;
  offset[2][0] = -1;
  offset[2][1] = 0;
  offset[2][2] = -1;
  offset[3][0] = 0;
  offset[3][1] = 0;
  offset[3][2] = -1;
  offset[4][0] = -1;
  offset[4][1] = -1;
  offset[4][2] = 0;
  offset[5][0] = 0;
  offset[5][1] = -1;
  offset[5][2] = 0;
  offset[6][0] = -1;
  offset[6][1] = 0;
  offset[6][2] = 0;
  offset[7][0] = 0;
  offset[7][1] = 0;
  offset[7][2] = 0;

  // For the rest of the points, trim the seed region
  // This is essentially an intersection of edge neighbors.
  for (i = 1; i < numPts; i++)
  {
    //  Get the location of the point
    id = ptIds->GetId(i);
    ptLoc[0] = id % dim[0];
    ptLoc[1] = (id / dim[0]) % dim[1];
    ptLoc[2] = id / (static_cast<vtkIdType>(dim[0]) * dim[1]);

    if ((ptLoc[0] - 1) == seedLoc[0])
    {
      offset[0][0] = -10;
      offset[2][0] = -10;
      offset[4][0] = -10;
      offset[6][0] = -10;
    }
    else if ((ptLoc[0] + 1) == seedLoc[0])
    {
      offset[1][0] = -10;
      offset[3][0] = -10;
      offset[5][0] = -10;
      offset[7][0] = -10;
    }
    else if ((ptLoc[1] - 1) == seedLoc[1])
    {
      offset[0][1] = -10;
      offset[1][1] = -10;
      offset[4][1] = -10;
      offset[5][1] = -10;
    }
    else if ((ptLoc[1] + 1) == seedLoc[1])
    {
      offset[2][1] = -10;
      offset[3][1] = -10;
      offset[6][1] = -10;
      offset[7][1] = -10;
    }
    else if ((ptLoc[2] - 1) == seedLoc[2])
    {
      offset[0][2] = -10;
      offset[1][2] = -10;
      offset[2][2] = -10;
      offset[3][2] = -10;
    }
    else if ((ptLoc[2] + 1) == seedLoc[2])
    {
      offset[4][2] = -10;
      offset[5][2] = -10;
      offset[6][2] = -10;
      offset[7][2] = -10;
    }
  }

  // Load the non-trimmed cells
  cellDim[0] = dim[0] - 1;
  cellDim[1] = dim[1] - 1;
  cellDim[2] = dim[2] - 1;

  for (i = 0; i < 3; i++)
  {
    cellDim[i] = std::max<vtkIdType>(cellDim[i], 1);
  }

  for (j = 0; j < 8; j++)
  {
    for (i = 0; i < 3; i++)
    {
      if (offset[j][i] != -10)
      {
        cellLoc[i] = seedLoc[i] + offset[j][i];
        if (cellLoc[i] < 0 || cellLoc[i] >= cellDim[i])
        {
          break;
        }
      }
      else
      {
        break;
      }
    }
    if (i >= 3) // add cell
    {
      id = cellLoc[0] + cellLoc[1] * cellDim[0] + cellLoc[2] * cellDim[0] * cellDim[1];
      if (id != cellId)
      {
        cellIds->InsertNextId(id);
      }
    }
  }
}
VTK_ABI_NAMESPACE_END