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

#include "vtkCellArray.h"
#include "vtkCellData.h"
#include "vtkExecutive.h"
#include "vtkFloatArray.h"
#include "vtkGenericCell.h"
#include "vtkIdList.h"
#include "vtkImageData.h"
#include "vtkIncrementalPointLocator.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkIntArray.h"
#include "vtkMergePoints.h"
#include "vtkObjectFactory.h"
#include "vtkPointData.h"
#include "vtkUnsignedCharArray.h"
#include "vtkUnstructuredGrid.h"

#include <cmath>
#include <vector>

VTK_ABI_NAMESPACE_BEGIN
vtkStandardNewMacro(vtkBoxClipDataSet);
vtkCxxSetObjectMacro(vtkBoxClipDataSet, Locator, vtkIncrementalPointLocator);
//------------------------------------------------------------------------------
vtkBoxClipDataSet::vtkBoxClipDataSet()
{
  this->Locator = nullptr;
  this->GenerateClipScalars = 0;

  this->GenerateClippedOutput = 0;
  // this->MergeTolerance = 0.01;

  this->SetNumberOfOutputPorts(2);

  this->Orientation = 1;

  this->PlaneNormal[0][0] = -1.0;
  this->PlaneNormal[0][1] = 0.0;
  this->PlaneNormal[0][2] = 0.0;

  this->PlaneNormal[1][0] = 1.0;
  this->PlaneNormal[1][1] = 0.0;
  this->PlaneNormal[1][2] = 0.0;

  this->PlaneNormal[2][0] = 0.0;
  this->PlaneNormal[2][1] = -1.0;
  this->PlaneNormal[2][2] = 0.0;

  this->PlaneNormal[3][0] = 0.0;
  this->PlaneNormal[3][1] = 1.0;
  this->PlaneNormal[3][2] = 0.0;

  this->PlaneNormal[4][0] = 0.0;
  this->PlaneNormal[4][1] = 0.0;
  this->PlaneNormal[4][2] = -1.0;

  this->PlaneNormal[5][0] = 0.0;
  this->PlaneNormal[5][1] = 0.0;
  this->PlaneNormal[5][2] = 1.0;

  this->PlanePoint[0][0] = 0.0;
  this->PlanePoint[0][1] = 0.0;
  this->PlanePoint[0][2] = 0.0;

  this->PlanePoint[1][0] = 1.0;
  this->PlanePoint[1][1] = 0.0;
  this->PlanePoint[1][2] = 0.0;

  this->PlanePoint[2][0] = 0.0;
  this->PlanePoint[2][1] = 0.0;
  this->PlanePoint[2][2] = 0.0;

  this->PlanePoint[3][0] = 0.0;
  this->PlanePoint[3][1] = 1.0;
  this->PlanePoint[3][2] = 0.0;

  this->PlanePoint[4][0] = 0.0;
  this->PlanePoint[4][1] = 0.0;
  this->PlanePoint[4][2] = 0.0;

  this->PlanePoint[5][0] = 0.0;
  this->PlanePoint[5][1] = 0.0;
  this->PlanePoint[5][2] = 1.0;

  int i = 0;
  while (i < 3)
  {
    this->BoundBoxClip[i][0] = 0.0;
    this->BoundBoxClip[i][1] = 1.0;
    ++i;
  }

  // by default process active point scalars
  this->SetInputArrayToProcess(
    0, 0, 0, vtkDataObject::FIELD_ASSOCIATION_POINTS, vtkDataSetAttributes::SCALARS);
}

//------------------------------------------------------------------------------
vtkBoxClipDataSet::~vtkBoxClipDataSet()
{
  this->SetLocator(nullptr);
}

//------------------------------------------------------------------------------
// Do not say we have two outputs unless we are generating the clipped output.
int vtkBoxClipDataSet::GetNumberOfOutputs()
{
  if (this->GenerateClippedOutput)
  {
    return 2;
  }
  return 1;
}

//------------------------------------------------------------------------------
// Overload standard modified time function. If Clip functions is modified,
// then this object is modified as well.
vtkMTimeType vtkBoxClipDataSet::GetMTime()
{
  vtkMTimeType mTime = this->Superclass::GetMTime();
  vtkMTimeType time;

  if (this->Locator != nullptr)
  {
    time = this->Locator->GetMTime();
    mTime = (time > mTime ? time : mTime);
  }

  return mTime;
}

//------------------------------------------------------------------------------
vtkUnstructuredGrid* vtkBoxClipDataSet::GetClippedOutput()
{
  if (this->GetNumberOfOutputPorts() < 2)
  {
    return nullptr;
  }

  return vtkUnstructuredGrid::SafeDownCast(this->GetExecutive()->GetOutputData(1));
}

//------------------------------------------------------------------------------
//
// Clip by box
//
int vtkBoxClipDataSet::RequestData(vtkInformation* vtkNotUsed(request),
  vtkInformationVector** inputVector, vtkInformationVector* outputVector)
{
  // get the info objects
  vtkInformation* inInfo = inputVector[0]->GetInformationObject(0);
  vtkInformation* outInfo = outputVector->GetInformationObject(0);

  // get the input and output
  vtkDataSet* input = vtkDataSet::SafeDownCast(inInfo->Get(vtkDataObject::DATA_OBJECT()));
  vtkUnstructuredGrid* output =
    vtkUnstructuredGrid::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT()));

  vtkUnstructuredGrid* clippedOutput = this->GetClippedOutput();

  vtkIdType i;
  vtkIdType npts;
  const vtkIdType* pts;
  vtkIdType estimatedSize;
  vtkIdType newCellId;
  vtkIdType numPts = input->GetNumberOfPoints();
  vtkIdType numCells = input->GetNumberOfCells();
  vtkPointData* inPD = input->GetPointData();
  vtkPointData* outPD[2];
  vtkCellData* inCD = input->GetCellData();
  vtkCellData* outCD[2];
  vtkPoints* newPoints;
  vtkPoints* cellPts;
  vtkDebugMacro(<< "Clip by Box\n");
  vtkUnsignedCharArray* types[2];

  int j;
  int cellType = 0;
  int numOutputs = 1;

  // Initialize self; create output objects
  //
  if (numPts < 1)
  {
    vtkDebugMacro(<< "No data to clip");
    return 1;
  }

  // allocate the output and associated helper classes
  estimatedSize = numCells;
  estimatedSize = estimatedSize / 1024 * 1024; // multiple of 1024
  estimatedSize = std::max<vtkIdType>(estimatedSize, 1024);
  vtkCellArray* conn[2];
  conn[0] = vtkCellArray::New();
  conn[0]->AllocateEstimate(estimatedSize, 1);
  conn[0]->InitTraversal();
  types[0] = vtkUnsignedCharArray::New();
  types[0]->Allocate(estimatedSize, estimatedSize / 2);

  if (this->GenerateClippedOutput)
  {
    numOutputs = 2;
    conn[1] = vtkCellArray::New();
    conn[1]->AllocateEstimate(estimatedSize, 1);
    conn[1]->InitTraversal();
    types[1] = vtkUnsignedCharArray::New();
    types[1]->Allocate(estimatedSize, estimatedSize / 2);
  }

  newPoints = vtkPoints::New();
  newPoints->Allocate(numPts, numPts / 2);

  // locator used to merge potentially duplicate points
  if (this->Locator == nullptr)
  {
    this->CreateDefaultLocator();
  }
  this->Locator->InitPointInsertion(newPoints, input->GetBounds());

  outPD[0] = output->GetPointData();
  if (this->GenerateClippedOutput)
  {
    outPD[1] = output->GetPointData();
  }

  vtkDataArray* scalars = this->GetInputArrayToProcess(0, inputVector);
  if (!this->GenerateClipScalars && !scalars)
  {
    outPD[0]->CopyScalarsOff();
    if (this->GenerateClippedOutput)
    {
      outPD[1]->CopyScalarsOff();
    }
  }
  else
  {
    outPD[0]->CopyScalarsOn();
    if (this->GenerateClippedOutput)
    {
      outPD[1]->CopyScalarsOn();
    }
  }
  outPD[0]->InterpolateAllocate(inPD, estimatedSize, estimatedSize / 2);
  outCD[0] = output->GetCellData();
  outCD[0]->CopyAllocate(inCD, estimatedSize, estimatedSize / 2);
  if (this->GenerateClippedOutput)
  {
    outPD[1] = clippedOutput->GetPointData();
    outPD[1]->InterpolateAllocate(inPD, estimatedSize, estimatedSize / 2);
    outCD[1] = clippedOutput->GetCellData();
    outCD[1]->CopyAllocate(inCD, estimatedSize, estimatedSize / 2);
  }

  // Process all cells and clip each in turn

  vtkIdType updateTime = numCells / 20 + 1; // update roughly every 5%
  vtkGenericCell* cell = vtkGenericCell::New();
  vtkIdType cellId;

  bool abort = false;
  int num[2];
  int numNew[2];

  num[0] = num[1] = 0;
  numNew[0] = numNew[1] = 0;

  unsigned int orientation = this->GetOrientation(); // Test if there is a transformation

  // clock_t init_tmp = clock();
  for (cellId = 0; cellId < numCells && !abort; cellId++)
  {
    if (!(cellId % updateTime))
    {
      this->UpdateProgress(static_cast<float>(cellId) / numCells);
      abort = this->CheckAbort();
    }

    input->GetCell(cellId, cell);
    cellPts = cell->GetPoints();
    npts = cellPts->GetNumberOfPoints();

    if (this->GenerateClippedOutput)
    {
      if ((cell->GetCellDimension()) == 3)
      {
        if (orientation)
        {
          this->ClipHexahedronInOut(
            newPoints, cell, this->Locator, conn, inPD, outPD, inCD, cellId, outCD);
        }
        else
        {
          this->ClipBoxInOut(
            newPoints, cell, this->Locator, conn, inPD, outPD, inCD, cellId, outCD);
        }

        numNew[0] = conn[0]->GetNumberOfCells() - num[0];
        numNew[1] = conn[1]->GetNumberOfCells() - num[1];
        num[0] = conn[0]->GetNumberOfCells();
        num[1] = conn[1]->GetNumberOfCells();
      }
      else if ((cell->GetCellDimension()) == 2)
      {
        if (orientation)
        {
          this->ClipHexahedronInOut2D(
            newPoints, cell, this->Locator, conn, inPD, outPD, inCD, cellId, outCD);
        }
        else
        {
          this->ClipBoxInOut2D(
            newPoints, cell, this->Locator, conn, inPD, outPD, inCD, cellId, outCD);
        }
        numNew[0] = conn[0]->GetNumberOfCells() - num[0];
        numNew[1] = conn[1]->GetNumberOfCells() - num[1];
        num[0] = conn[0]->GetNumberOfCells();
        num[1] = conn[1]->GetNumberOfCells();
      }
      else if (cell->GetCellDimension() == 1)
      {
        if (orientation)
        {
          this->ClipHexahedronInOut1D(
            newPoints, cell, this->Locator, conn, inPD, outPD, inCD, cellId, outCD);
        }
        else
        {
          this->ClipBoxInOut1D(
            newPoints, cell, this->Locator, conn, inPD, outPD, inCD, cellId, outCD);
        }
        numNew[0] = conn[0]->GetNumberOfCells() - num[0];
        numNew[1] = conn[1]->GetNumberOfCells() - num[1];
        num[0] = conn[0]->GetNumberOfCells();
        num[1] = conn[1]->GetNumberOfCells();
      }
      else if (cell->GetCellDimension() == 0)
      {
        if (orientation)
        {
          this->ClipHexahedronInOut0D(cell, this->Locator, conn, inPD, outPD, inCD, cellId, outCD);
        }
        else
        {
          this->ClipBoxInOut0D(cell, this->Locator, conn, inPD, outPD, inCD, cellId, outCD);
        }
        numNew[0] = conn[0]->GetNumberOfCells() - num[0];
        numNew[1] = conn[1]->GetNumberOfCells() - num[1];
        num[0] = conn[0]->GetNumberOfCells();
        num[1] = conn[1]->GetNumberOfCells();
      }
      else
      {
        vtkErrorMacro(<< "Do not support cells of dimension " << cell->GetCellDimension());
      }
    }
    else
    {
      if ((cell->GetCellDimension()) == 3)
      {
        if (orientation)
        {
          this->ClipHexahedron(
            newPoints, cell, this->Locator, conn[0], inPD, outPD[0], inCD, cellId, outCD[0]);
        }
        else
        {
          this->ClipBox(
            newPoints, cell, this->Locator, conn[0], inPD, outPD[0], inCD, cellId, outCD[0]);
        }

        numNew[0] = conn[0]->GetNumberOfCells() - num[0];
        num[0] = conn[0]->GetNumberOfCells();
      }
      else if ((cell->GetCellDimension()) == 2)
      {
        if (orientation)
        {
          this->ClipHexahedron2D(
            newPoints, cell, this->Locator, conn[0], inPD, outPD[0], inCD, cellId, outCD[0]);
        }
        else
        {
          this->ClipBox2D(
            newPoints, cell, this->Locator, conn[0], inPD, outPD[0], inCD, cellId, outCD[0]);
        }
        numNew[0] = conn[0]->GetNumberOfCells() - num[0];
        num[0] = conn[0]->GetNumberOfCells();
      }
      else if (cell->GetCellDimension() == 1)
      {
        if (orientation)
        {
          this->ClipHexahedron1D(
            newPoints, cell, this->Locator, conn[0], inPD, outPD[0], inCD, cellId, outCD[0]);
        }
        else
        {
          this->ClipBox1D(
            newPoints, cell, this->Locator, conn[0], inPD, outPD[0], inCD, cellId, outCD[0]);
        }
        numNew[0] = conn[0]->GetNumberOfCells() - num[0];
        num[0] = conn[0]->GetNumberOfCells();
      }
      else if (cell->GetCellDimension() == 0)
      {
        if (orientation)
        {
          this->ClipHexahedron0D(
            cell, this->Locator, conn[0], inPD, outPD[0], inCD, cellId, outCD[0]);
        }
        else
        {
          this->ClipBox0D(cell, this->Locator, conn[0], inPD, outPD[0], inCD, cellId, outCD[0]);
        }
        numNew[0] = conn[0]->GetNumberOfCells() - num[0];
        num[0] = conn[0]->GetNumberOfCells();
      }
      else
      {
        vtkErrorMacro(<< "Do not support cells of dimension " << cell->GetCellDimension());
      }
    }

    for (i = 0; i < numOutputs; i++) // for both outputs
    {
      for (j = 0; j < numNew[i]; j++)
      {
        conn[i]->GetNextCell(npts, pts);

        // For each new cell added, got to set the type of the cell
        switch (cell->GetCellDimension())
        {
          case 0: // points are generated-------------------------------
            cellType = (npts > 1 ? VTK_POLY_VERTEX : VTK_VERTEX);
            break;

          case 1: // lines are generated----------------------------------
            cellType = (npts > 2 ? VTK_POLY_LINE : VTK_LINE);
            break;

          case 2: // polygons are generated------------------------------
                  // NOLINTNEXTLINE(readability-avoid-nested-conditional-operator)
            cellType = (npts == 3 ? VTK_TRIANGLE : (npts == 4 ? VTK_QUAD : VTK_POLYGON));
            break;

          case 3: // tetrahedra are generated------------------------------
            cellType = VTK_TETRA;
            break;
        } // switch

        newCellId = types[i]->InsertNextValue(cellType);
        outCD[i]->CopyData(inCD, cellId, newCellId);
      } // for each new cell
    }   // for both outputs
  }     // for each cell

  cell->Delete();

  output->SetPoints(newPoints);
  output->SetCells(types[0], conn[0]);

  conn[0]->Delete();
  types[0]->Delete();

  if (this->GenerateClippedOutput)
  {
    clippedOutput->SetPoints(newPoints);
    clippedOutput->SetCells(types[1], conn[1]);
    conn[1]->Delete();
    types[1]->Delete();
  }
  newPoints->Delete();
  this->Locator->Initialize(); // release any extra memory
  output->Squeeze();

  return 1;
}

//------------------------------------------------------------------------------
// Specify a spatial locator for merging points. By default,
// an instance of vtkMergePoints is used.
void vtkBoxClipDataSet::CreateDefaultLocator()
{
  if (this->Locator == nullptr)
  {
    this->Locator = vtkMergePoints::New();
    this->Locator->Register(this);
    this->Locator->Delete();
  }
}

//------------------------------------------------------------------------------
// Set the box for clipping
// for each plane, specify the normal and one vertex on the plane.
//
void vtkBoxClipDataSet::SetBoxClip(const double* n0, const double* o0, const double* n1,
  const double* o1, const double* n2, const double* o2, const double* n3, const double* o3,
  const double* n4, const double* o4, const double* n5, const double* o5)
{
  if ((this->Orientation == 1)

    && (this->PlaneNormal[0][0] == n0[0]) && (this->PlaneNormal[0][1] == n0[1]) &&
    (this->PlaneNormal[0][2] == n0[2])

    && (this->PlaneNormal[1][0] == n1[0]) && (this->PlaneNormal[1][1] == n1[1]) &&
    (this->PlaneNormal[1][2] == n1[2])

    && (this->PlaneNormal[2][0] == n2[0]) && (this->PlaneNormal[2][1] == n2[1]) &&
    (this->PlaneNormal[2][2] == n2[2])

    && (this->PlaneNormal[3][0] == n3[0]) && (this->PlaneNormal[3][1] == n3[1]) &&
    (this->PlaneNormal[3][2] == n3[2])

    && (this->PlaneNormal[4][0] == n4[0]) && (this->PlaneNormal[4][1] == n4[1]) &&
    (this->PlaneNormal[4][2] == n4[2])

    && (this->PlaneNormal[5][0] == n5[0]) && (this->PlaneNormal[5][1] == n5[1]) &&
    (this->PlaneNormal[5][2] == n5[2])

    && (this->PlanePoint[0][0] == o0[0]) && (this->PlanePoint[0][1] == o0[1]) &&
    (this->PlanePoint[0][2] == o0[2])

    && (this->PlanePoint[1][0] == o1[0]) && (this->PlanePoint[1][1] == o1[1]) &&
    (this->PlanePoint[1][2] == o1[2])

    && (this->PlanePoint[2][0] == o2[0]) && (this->PlanePoint[2][1] == o2[1]) &&
    (this->PlanePoint[2][2] == o2[2])

    && (this->PlanePoint[3][0] == o3[0]) && (this->PlanePoint[3][1] == o3[1]) &&
    (this->PlanePoint[3][2] == o3[2])

    && (this->PlanePoint[4][0] == o4[0]) && (this->PlanePoint[4][1] == o4[1]) &&
    (this->PlanePoint[4][2] == o4[2])

    && (this->PlanePoint[5][0] == o5[0]) && (this->PlanePoint[5][1] == o5[1]) &&
    (this->PlanePoint[5][2] == o5[2]))
  {
    return;
  }

  this->SetOrientation(1);

  this->PlaneNormal[0][0] = n0[0];
  this->PlaneNormal[0][1] = n0[1];
  this->PlaneNormal[0][2] = n0[2];

  this->PlaneNormal[1][0] = n1[0];
  this->PlaneNormal[1][1] = n1[1];
  this->PlaneNormal[1][2] = n1[2];

  this->PlaneNormal[2][0] = n2[0];
  this->PlaneNormal[2][1] = n2[1];
  this->PlaneNormal[2][2] = n2[2];

  this->PlaneNormal[3][0] = n3[0];
  this->PlaneNormal[3][1] = n3[1];
  this->PlaneNormal[3][2] = n3[2];

  this->PlaneNormal[4][0] = n4[0];
  this->PlaneNormal[4][1] = n4[1];
  this->PlaneNormal[4][2] = n4[2];

  this->PlaneNormal[5][0] = n5[0];
  this->PlaneNormal[5][1] = n5[1];
  this->PlaneNormal[5][2] = n5[2];

  this->PlanePoint[0][0] = o0[0];
  this->PlanePoint[0][1] = o0[1];
  this->PlanePoint[0][2] = o0[2];

  this->PlanePoint[1][0] = o1[0];
  this->PlanePoint[1][1] = o1[1];
  this->PlanePoint[1][2] = o1[2];

  this->PlanePoint[2][0] = o2[0];
  this->PlanePoint[2][1] = o2[1];
  this->PlanePoint[2][2] = o2[2];

  this->PlanePoint[3][0] = o3[0];
  this->PlanePoint[3][1] = o3[1];
  this->PlanePoint[3][2] = o3[2];

  this->PlanePoint[4][0] = o4[0];
  this->PlanePoint[4][1] = o4[1];
  this->PlanePoint[4][2] = o4[2];

  this->PlanePoint[5][0] = o5[0];
  this->PlanePoint[5][1] = o5[1];
  this->PlanePoint[5][2] = o5[2];

  this->Modified();
}

//------------------------------------------------------------------------------
// Specify the bounding box for clipping

void vtkBoxClipDataSet::SetBoxClip(
  double xmin, double xmax, double ymin, double ymax, double zmin, double zmax)
{
  if ((this->Orientation == 0) && (this->BoundBoxClip[0][0] == xmin) &&
    (this->BoundBoxClip[0][1] == xmax) && (this->BoundBoxClip[1][0] == ymin) &&
    (this->BoundBoxClip[1][1] == ymax) && (this->BoundBoxClip[2][0] == zmin) &&
    (this->BoundBoxClip[2][1] == zmax))
  {
    return;
  }

  this->SetOrientation(0);
  this->BoundBoxClip[0][0] = xmin;
  this->BoundBoxClip[0][1] = xmax;
  this->BoundBoxClip[1][0] = ymin;
  this->BoundBoxClip[1][1] = ymax;
  this->BoundBoxClip[2][0] = zmin;
  this->BoundBoxClip[2][1] = zmax;

  this->Modified();
}

//------------------------------------------------------------------------------
int vtkBoxClipDataSet::FillInputPortInformation(int, vtkInformation* info)
{
  info->Set(vtkAlgorithm::INPUT_REQUIRED_DATA_TYPE(), "vtkDataSet");
  return 1;
}

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

  // os << indent << "Merge Tolerance: " << this->MergeTolerance << "\n";
  os << indent << "Orientation: " << this->Orientation << "\n";

  if (this->Locator)
  {
    os << indent << "Locator: " << this->Locator << "\n";
  }
  else
  {
    os << indent << "Locator: (none)\n";
  }

  os << indent << "Generate Clipped Output: " << (this->GenerateClippedOutput ? "Yes\n" : "Off\n");
  os << indent << "Generate Clip Scalars: " << (this->GenerateClipScalars ? "On\n" : "Off\n");
}

//------------------------------------------------------------------------------
// InterpolateEdge: Interpolate the data in a vtkDataSetAttributes along a line
//
// This method works very much like vtkDataSetAttributes::InterpolateEdge
// except that rather than take the interpolation information from an
// input and copy it to an output, the values to interpolate are already
// placed in the output arrays.  This is necessary because vtkBoxClipDataSet
// will continually cut (and consequently interpolate) the input until it
// fits within the bounds.
void vtkBoxClipDataSet::InterpolateEdge(
  vtkDataSetAttributes* attributes, vtkIdType toId, vtkIdType fromId1, vtkIdType fromId2, double t)
{
  int numArrays = attributes->GetNumberOfArrays();
  for (int i = 0; i < numArrays; i++)
  {
    vtkAbstractArray* array = attributes->GetAbstractArray(i);

    // We are ignoring any special interpolate flags (such as nearest neighbor
    // interpolation).  That kind of interpolation is not linear and could be
    // incorrect when multiple cuts are performed on the same primitive (which
    // can happen).

    array->InterpolateTuple(toId, fromId1, array, fromId2, array, t);
  }
}

//------------------------------------------------------------------------------
// CellGrid: Subdivide cells in consistent tetrahedra.
// Case : Voxel(11) or Hexahedron(12).
//
// MinEdgF search the smallest vertex index in linear order of a face(4 vertices)
//
void vtkBoxClipDataSet::MinEdgeF(
  const unsigned int* id_v, const vtkIdType* cellIds, unsigned int* edgF)
{

  int i;
  unsigned int id;
  int ids;
  int min_f;

  ids = 0;
  id = id_v[0]; // Face index
  min_f = cellIds[id_v[0]];

  for (i = 1; i < 4; i++)
  {
    if (min_f > cellIds[id_v[i]])
    {
      min_f = cellIds[id_v[i]];
      id = id_v[i];
      ids = i;
    }
  }

  switch (ids)
  {
    case 0:
      if (id < id_v[2])
      {
        edgF[0] = id;
        edgF[1] = id_v[2];
      }
      else
      {
        edgF[0] = id_v[2];
        edgF[1] = id;
      }
      break;
    case 1:
      if (id < id_v[3])
      {
        edgF[0] = id;
        edgF[1] = id_v[3];
      }
      else
      {
        edgF[0] = id_v[3];
        edgF[1] = id;
      }
      break;
    case 2:
      if (id < id_v[0])
      {
        edgF[0] = id;
        edgF[1] = id_v[0];
      }
      else
      {
        edgF[0] = id_v[0];
        edgF[1] = id;
      }
      break;
    case 3:
      if (id < id_v[1])
      {
        edgF[0] = id;
        edgF[1] = id_v[1];
      }
      else
      {
        edgF[0] = id_v[1];
        edgF[1] = id;
      }
      break;
  }
}

//------------------------------------------------------------------------------
// CellGrid: Subdivide cells in consistent tetrahedra.
//
// Case : Voxel or Hexahedron:
//       if ( subdivide voxel in 6 tetrahedra)
//           voxel : 2 wedges (*)
//       else
//           voxel : 5 tetrahedra
//
// Case : Wedge (*)
//
// ------------------------------------------------
//
//(*) WedgeToTetra: subdivide one wedge in 3 tetrahedra
//
//        wedge : 1 tetrahedron + 1 pyramid = 3 tetrahedra.
//
void vtkBoxClipDataSet::WedgeToTetra(
  const vtkIdType* wedgeId, const vtkIdType* cellIds, vtkCellArray* newCellArray)
{
  int i;
  int id;
  vtkIdType xmin;
  vtkIdType tab[4];
  vtkIdType tabpyram[5];

  constexpr vtkIdType vwedge[6][4] = {
    { 0, 4, 3, 5 },
    { 1, 4, 3, 5 },
    { 2, 4, 3, 5 },
    { 3, 0, 1, 2 },
    { 4, 0, 1, 2 },
    { 5, 0, 1, 2 },
  };

  // the table 'vwedge' set 6 possibilities of the smallest index
  //
  //             v5
  //             /\       .
  //         v3 /..\ v4
  //           /   /
  //        v2/\  /
  //       v0/__\/v1
  //   if(v0 index ) is the smallest index:
  //   wedge is subdivided in:
  //         1 tetrahedron-> vwedge[0]: {v0,v4,v3,v5}
  //     and 1 pyramid    -> vert[0]  : {v1,v2,v5,v4,v0}
  //

  id = 0;
  xmin = cellIds[wedgeId[0]];
  for (i = 1; i < 6; i++)
  {
    if (xmin > cellIds[wedgeId[i]])
    {
      xmin = cellIds[wedgeId[i]]; // the smallest global index
      id = i;                     // local index
    }
  }
  for (i = 0; i < 4; i++)
  {
    tab[i] = wedgeId[vwedge[id][i]];
  }
  newCellArray->InsertNextCell(4, tab);

  // Pyramid :create 2 tetrahedra
  constexpr vtkIdType vert[6][5] = {
    { 1, 2, 5, 4, 0 },
    { 2, 0, 3, 5, 1 },
    { 3, 0, 1, 4, 2 },
    { 1, 2, 5, 4, 3 },
    { 2, 0, 3, 5, 4 },
    { 3, 0, 1, 4, 5 },
  };
  for (i = 0; i < 5; i++)
  {
    tabpyram[i] = wedgeId[vert[id][i]];
  }
  this->PyramidToTetra(tabpyram, cellIds, newCellArray);
}

//------------------------------------------------------------------------------
// CellGrid: Subdivide cells in consistent tetrahedra.
//
// PyramidToTetra :Subdivide the pyramid in consistent tetrahedra.
//        Pyramid : 2 tetrahedra.
//
void vtkBoxClipDataSet::PyramidToTetra(
  const vtkIdType* pyramId, const vtkIdType* cellIds, vtkCellArray* newCellArray)
{
  vtkIdType xmin;
  unsigned int i, j, idpy;
  vtkIdType tab[4];

  // the table 'vpy' set 3 possibilities of the smallest index
  // vertices{v0,v1,v2,v3,v4}. {v0,v1,v2,v3} is a square face of pyramid
  //
  //                v4
  //                ^
  //
  //
  //           v3 _ _ __ _  v2
  //           /         /
  //        v0/_ _ _ _ _/v1
  //
  //   if(v0 index ) is the smallest index:
  //   the pyramid is subdivided in:
  //         2 tetrahedra-> vpy[0]: {v0,v1,v2,v4}
  //                        vpy[1]: {v0,v2,v3,v4}
  //
  constexpr vtkIdType vpy[8][4] = {
    { 0, 1, 2, 4 },
    { 0, 2, 3, 4 },
    { 1, 2, 3, 4 },
    { 1, 3, 0, 4 },
    { 2, 3, 0, 4 },
    { 2, 0, 1, 4 },
    { 3, 0, 1, 4 },
    { 3, 1, 2, 4 },
  };

  xmin = cellIds[pyramId[0]];
  idpy = 0;
  for (i = 1; i < 4; i++)
  {
    if (xmin > cellIds[pyramId[i]])
    {
      xmin = cellIds[pyramId[i]]; // global index
      idpy = i;                   // local index
    }
  }
  for (j = 0; j < 4; j++)
  {
    tab[j] = pyramId[vpy[2 * idpy][j]];
  }
  newCellArray->InsertNextCell(4, tab);

  for (j = 0; j < 4; j++)
  {
    tab[j] = pyramId[vpy[2 * idpy + 1][j]];
  }
  newCellArray->InsertNextCell(4, tab);
}

//------------------------------------------------------------------------------
// Tetra Grid : Subdivide cells in  consistent tetrahedra.
//             For each cell, search the smallest global index.
//
//  Case Tetrahedron(10): Just insert this cell in the newCellArray
//
//  Case Voxel(11) or Hexahedron(12):
//       - for each face: looking for the diagonal edge with the smallest index
//       - 2 possibilities: subdivide a cell in 5 or 6 tetrahedra
//
//         (I)Case 6 tetrahedra:
//                  - 6 possibilities: subdivide de cell in 2 wedges:
//                   (1) diagonal edges (v0,v5),(v2,v7)
//                       vwedge[0]={0,5,4,2,7,6},
//                       vwedge[1]={0,1,5,2,3,7},
//                   (2) diagonal edges (v4,v7),(v0,v3)
//                       vwedge[2]={4,7,6,0,3,2},
//                       vwedge[3]={4,5,7,0,1,3},
//                   (3) diagonal edges (v0,v6),(v1,v7)
//        vwedge[4]= {1,7,5,0,6,4},
//                       vwedge[5]{1,3,7,0,2,6},
//                       subdivide de cell in 2 wedges:
//                   (4) diagonal edges (v1,v2),(v5,v6)
//                       vwedge[6]={4,5,6,0,1,2},
//                       vwedge[7]={6,5,7,2,1,3},
//                   (5) diagonal edges (v2,v4),(v3,v5)
//                       vwedge[8]={3,7,5,2,6,4},
//                       vwedge[9]={1,3,5,0,2,4},
//                   (6) diagonal edges (v1,v4),(v3,v6)
//                       vwedge[10]={0,1,4,2,3,6}
//                       vwedge[11]={1,5,4,3,7,6}
//         v6 _ _ __ _  v7
//           /|        /|           VOXEL
//        v4/_|_ _ _ _/ |           opposite vertex of v0 is v7 and vice-versa
//          | |     v5| |           diagonal edges Edg_f[i]
//          |v2 _ _ _ |_|  v3
//          |/        |/
//        v0/_ _ _ _ _|/v1
//
//
//     (II)Case 5 tetrahedra:
//       - search the smallest vertex vi
//       - verify if the opposite vertices of vi do not belong to any diagonal edges Edg_f
//       - 2 possibilities: create 5 tetraedra
//           - if vi is ( 0 or 3 or 5 or 6)
//             vtetra[]={v0,v5,v3,v6},{v0,v4,v5,v6},{v0,v1,v3,v5},{v5,v3,v6,v7},{v0,v3,v2,v6}};
//           - if vi is ( 1 or 2 or 4 or 7)
//             vtetra[]={v1,v2,v4,v7},{v0,v1,v2,v4},{v1,v4,v5,v7},{v1,v3,v2,v7},{v2,v6,v4,v7}};
//  Case Wedge (13):
//       the table 'vwedge' set 6 possibilities of the smallest index
//
//                 v5
//               /\            .
//           v3 /..\ v4
//             /   /
//          v2/\  /
//         v0/__\/v1
//
//        if(v0 index ) is the smallest index:
//         wedge is subdivided in:
//              1 tetrahedron-> vwedge[0]: {v0,v4,v3,v5}
//               and 1 pyramid-> vert[0] : {v1,v2,v5,v4,v0}
//
//  Case Pyramid (14):
//       the table 'vpy' set 3 possibilities of the smallest index
//       vertices{v0,v1,v2,v3,v4}. {v0,v1,v2,v3} is a square face of pyramid
//
//                v4  (opposite vertex of face with 4 vertices)
//                ^
//
//
//          v3 _ _ __ _  v2
//           /         /
//        v0/_ _ _ _ _/v1
//
//        if(v0 index ) is the smallest index:
//        the pyramid is subdivided in:
//        2 tetrahedra-> vpyram[0]: {v0,v1,v2,v4}
//                       vpyram[1]: {v0,v2,v3,v4}
//
//
//------------------------------------------------------------------------------
void vtkBoxClipDataSet::CellGrid(
  vtkIdType typeobj, vtkIdType npts, const vtkIdType* cellIds, vtkCellArray* newCellArray)
{
  vtkIdType tab[4];
  vtkIdType tabp[5];
  vtkIdType ptstriangle = 3;
  vtkIdType ptstetra = 4;
  vtkIdType xmin;
  vtkIdType idt;
  int i, j;
  unsigned int id = 0;
  unsigned int idpy = 0;

  unsigned int Edg_f[6][2]; // edge selected of each face
  unsigned int idv[4];

  unsigned int idopos;
  unsigned int numbertetra;

  constexpr vtkIdType triPassThrough[3] = { 0, 1, 2 };
  vtkIdType tri[3];
  vtkIdType line[2];

  switch (typeobj)
  {
    case VTK_VERTEX:
    case VTK_POLY_VERTEX:
      for (idt = 0; idt < npts; idt++)
      {
        newCellArray->InsertNextCell(1, &idt);
      }
      break;

    case VTK_LINE:
    case VTK_POLY_LINE:
      for (idt = 0; idt < npts - 1; idt++)
      {
        line[0] = idt;
        line[1] = idt + 1;
        newCellArray->InsertNextCell(2, line);
      }
      break;

    case VTK_TRIANGLE: // 5
    case VTK_QUADRATIC_TRIANGLE:
    case VTK_BIQUADRATIC_TRIANGLE:
      newCellArray->InsertNextCell(ptstriangle, triPassThrough);
      break;

    case VTK_TRIANGLE_STRIP: // 6
      for (idt = 0; idt < npts - 2; idt++)
      {
        if (idt % 2 == 0)
        {
          tri[0] = idt;
          tri[1] = idt + 1;
          tri[2] = idt + 2;
        }
        else
        {
          tri[0] = idt;
          tri[1] = idt + 2;
          tri[2] = idt + 1;
        }
        newCellArray->InsertNextCell(3, tri);
      }
      break;

    case VTK_POLYGON: // 7 (Convex case)
      tri[0] = 0;
      for (idt = 2; idt < npts; idt++)
      {
        tri[1] = idt - 1;
        tri[2] = idt;
        newCellArray->InsertNextCell(3, tri);
      }
      break;

    case VTK_PIXEL: // 8
    {
      constexpr vtkIdType vtrip[2][3] = {
        { 0, 1, 3 },
        { 0, 3, 2 },
      };
      newCellArray->InsertNextCell(3, vtrip[0]);
      newCellArray->InsertNextCell(3, vtrip[1]);
    }
    break;

    case VTK_QUAD: // 9
    case VTK_QUADRATIC_QUAD:
    case VTK_BIQUADRATIC_QUAD:
    case VTK_QUADRATIC_LINEAR_QUAD:
    {
      constexpr vtkIdType vtriq[2][3] = {
        { 0, 1, 2 },
        { 0, 2, 3 },
      };
      newCellArray->InsertNextCell(3, vtriq[0]);
      newCellArray->InsertNextCell(3, vtriq[1]);
    }
    break;

    case VTK_TETRA: // 10
    case VTK_QUADRATIC_TETRA:
    {
      constexpr vtkIdType tetra[4] = { 0, 1, 2, 3 };
      newCellArray->InsertNextCell(ptstetra, tetra);
    }
    break;

    case VTK_VOXEL: // 11
      // each face: search edge with smallest global index
      // face 0 (0,1,5,4)
      idv[0] = 0;
      idv[1] = 1;
      idv[2] = 5;
      idv[3] = 4;
      this->MinEdgeF(idv, cellIds, Edg_f[0]);

      // face 1 (0,1,3,2)
      idv[0] = 0;
      idv[1] = 1;
      idv[2] = 3;
      idv[3] = 2;
      this->MinEdgeF(idv, cellIds, Edg_f[1]);
      // face 2 (0,2,6,4)
      idv[0] = 0;
      idv[1] = 2;
      idv[2] = 6;
      idv[3] = 4;
      this->MinEdgeF(idv, cellIds, Edg_f[2]);
      // face 3 (4,5,7,6)
      idv[0] = 4;
      idv[1] = 5;
      idv[2] = 7;
      idv[3] = 6;
      this->MinEdgeF(idv, cellIds, Edg_f[3]);
      // face 4 (2,3,7,6)
      idv[0] = 2;
      idv[1] = 3;
      idv[2] = 7;
      idv[3] = 6;
      this->MinEdgeF(idv, cellIds, Edg_f[4]);
      // face 5 (1,3,7,5)
      idv[0] = 1;
      idv[1] = 3;
      idv[2] = 7;
      idv[3] = 5;
      this->MinEdgeF(idv, cellIds, Edg_f[5]);

      // search the smallest global index of voxel
      xmin = cellIds[0];
      id = 0;
      for (i = 1; i < 8; i++)
      {
        if (xmin > cellIds[i])
        {
          xmin = cellIds[i]; // the smallest global index
          id = i;            // local index
        }
      }
      // two cases:
      idopos = 7 - id;
      numbertetra = 5;
      for (i = 0; i < 6; i++)
      {
        j = 0;
        if (idopos == Edg_f[i][j])
        {
          numbertetra = 6;
          break;
        }
        j = 1;
        if (idopos == Edg_f[i][j])
        {
          numbertetra = 6;
          break;
        }
      }

      if (numbertetra == 5)
      {
        // case 1: create  5 tetraedra
        if ((id == 0) || (id == 3) || (id == 5) || (id == 6))
        {
          constexpr vtkIdType vtetra[5][4] = {
            { 0, 5, 3, 6 },
            { 0, 4, 5, 6 },
            { 0, 1, 3, 5 },
            { 5, 3, 6, 7 },
            { 0, 3, 2, 6 },
          };
          for (i = 0; i < 5; i++)
          {
            newCellArray->InsertNextCell(4, vtetra[i]);
          }
        }
        else
        {
          constexpr vtkIdType vtetra[5][4] = {
            { 1, 2, 4, 7 },
            { 0, 1, 2, 4 },
            { 1, 4, 5, 7 },
            { 1, 3, 2, 7 },
            { 2, 6, 4, 7 },
          };

          for (i = 0; i < 5; i++)
          {
            newCellArray->InsertNextCell(4, vtetra[i]);
          }
        }
      }
      else
      {
        // case 2: create 2 wedges-> 6 tetrahedra
        constexpr vtkIdType vwedge[12][6] = {
          { 0, 5, 4, 2, 7, 6 },
          { 0, 1, 5, 2, 3, 7 },
          { 4, 7, 6, 0, 3, 2 },
          { 4, 5, 7, 0, 1, 3 },
          { 1, 7, 5, 0, 6, 4 },
          { 1, 3, 7, 0, 2, 6 },
          { 4, 5, 6, 0, 1, 2 },
          { 6, 5, 7, 2, 1, 3 },
          { 3, 7, 5, 2, 6, 4 },
          { 1, 3, 5, 0, 2, 4 },
          { 0, 1, 4, 2, 3, 6 },
          { 1, 5, 4, 3, 7, 6 },
        };
        unsigned int edgeId = 10 * Edg_f[i][0] + Edg_f[i][1];
        switch (edgeId)
        {
          case 5:  // edge(v0,v5):10*0 + 5
          case 27: // edge(v2,v7):10*2 + 7
            this->WedgeToTetra(vwedge[0], cellIds, newCellArray);
            this->WedgeToTetra(vwedge[1], cellIds, newCellArray);
            break;
          case 3:  // edge(v0,v2)
          case 47: // edge(v4,v6)
            this->WedgeToTetra(vwedge[2], cellIds, newCellArray);
            this->WedgeToTetra(vwedge[3], cellIds, newCellArray);
            break;
          case 6:
          case 17:
            this->WedgeToTetra(vwedge[4], cellIds, newCellArray);
            this->WedgeToTetra(vwedge[5], cellIds, newCellArray);
            break;
          case 12:
          case 56:
            this->WedgeToTetra(vwedge[6], cellIds, newCellArray);
            this->WedgeToTetra(vwedge[7], cellIds, newCellArray);
            break;
          case 24:
          case 35:
            this->WedgeToTetra(vwedge[8], cellIds, newCellArray);
            this->WedgeToTetra(vwedge[9], cellIds, newCellArray);
            break;
          case 14:
          case 36:
            this->WedgeToTetra(vwedge[10], cellIds, newCellArray);
            this->WedgeToTetra(vwedge[11], cellIds, newCellArray);
            break;
        }
      }
      break;

    case VTK_HEXAHEDRON: // 12
    case VTK_QUADRATIC_HEXAHEDRON:
    case VTK_TRIQUADRATIC_HEXAHEDRON:
    case VTK_BIQUADRATIC_QUADRATIC_HEXAHEDRON:
    {
      // each face: search edge with smallest global index
      // face 0 (0,1,5,4)
      idv[0] = 0;
      idv[1] = 1;
      idv[2] = 5;
      idv[3] = 4;
      this->MinEdgeF(idv, cellIds, Edg_f[0]);

      // face 1 (0,1,2,3)
      idv[0] = 0;
      idv[1] = 1;
      idv[2] = 2;
      idv[3] = 3;
      this->MinEdgeF(idv, cellIds, Edg_f[1]);
      // face 2 (0,3,7,4)
      idv[0] = 0;
      idv[1] = 3;
      idv[2] = 7;
      idv[3] = 4;
      this->MinEdgeF(idv, cellIds, Edg_f[2]);
      // face 3 (4,5,6,7)
      idv[0] = 4;
      idv[1] = 5;
      idv[2] = 6;
      idv[3] = 7;
      this->MinEdgeF(idv, cellIds, Edg_f[3]);
      // face 4 (3,2,6,7)
      idv[0] = 3;
      idv[1] = 2;
      idv[2] = 6;
      idv[3] = 7;
      this->MinEdgeF(idv, cellIds, Edg_f[4]);
      // face 5 (1,2,6,5)
      idv[0] = 1;
      idv[1] = 2;
      idv[2] = 6;
      idv[3] = 5;
      this->MinEdgeF(idv, cellIds, Edg_f[5]);

      // search the smallest global index of voxel
      xmin = cellIds[0];
      id = 0;
      for (i = 1; i < 8; i++)
      {
        if (xmin > cellIds[i])
        {
          xmin = cellIds[i]; // the smallest global index
          id = i;            // local index
        }
      }

      // two cases:
      constexpr unsigned int tabopos[8] = { 6, 7, 4, 5, 2, 3, 0, 1 };
      idopos = tabopos[id];
      numbertetra = 5;
      for (i = 0; i < 6; i++)
      {
        j = 0;
        if (idopos == Edg_f[i][j])
        {
          numbertetra = 6;
          break;
        }
        j = 1;
        if (idopos == Edg_f[i][j])
        {
          numbertetra = 6;
          break;
        }
      }

      if (numbertetra == 5)
      {
        // case 1: create  5 tetraedra
        if ((id == 0) || (id == 2) || (id == 5) || (id == 7))
        {
          constexpr vtkIdType vtetra[5][4] = {
            { 0, 5, 2, 7 },
            { 0, 4, 5, 7 },
            { 0, 1, 2, 5 },
            { 5, 2, 7, 6 },
            { 0, 2, 3, 7 },
          };
          for (i = 0; i < 5; i++)
          {
            newCellArray->InsertNextCell(4, vtetra[i]);
          }
        }
        else
        {
          constexpr vtkIdType vtetra[5][4] = {
            { 1, 3, 4, 6 },
            { 0, 1, 3, 4 },
            { 1, 4, 5, 6 },
            { 1, 2, 3, 6 },
            { 3, 7, 4, 6 },
          };

          for (i = 0; i < 5; i++)
          {
            newCellArray->InsertNextCell(4, vtetra[i]);
          }
        }
      }
      else
      {
        // case 2: create 2 wedges-> 6 tetrahedra
        constexpr vtkIdType vwedge[12][6] = {
          { 0, 5, 4, 3, 6, 7 },
          { 0, 1, 5, 3, 2, 6 },
          { 4, 6, 7, 0, 2, 3 },
          { 4, 5, 6, 0, 1, 2 },
          { 1, 6, 5, 0, 7, 4 },
          { 1, 2, 6, 0, 3, 7 },
          { 4, 5, 7, 0, 1, 3 },
          { 7, 5, 6, 3, 1, 2 },
          { 2, 6, 5, 3, 7, 4 },
          { 1, 2, 5, 0, 3, 4 },
          { 0, 1, 4, 3, 2, 7 },
          { 1, 5, 4, 2, 6, 7 },
        };
        unsigned int edgeId = 10 * Edg_f[i][0] + Edg_f[i][1];

        switch (edgeId)
        {
          case 5:  // edge(v0,v5):10*0 + 5
          case 36: // edge(v3,v6):10*3 + 6
            this->WedgeToTetra(vwedge[0], cellIds, newCellArray);
            this->WedgeToTetra(vwedge[1], cellIds, newCellArray);
            break;
          case 2:  // edge(v0,v2)
          case 46: // edge(v4,v6)
            this->WedgeToTetra(vwedge[2], cellIds, newCellArray);
            this->WedgeToTetra(vwedge[3], cellIds, newCellArray);
            break;
          case 7:
          case 16:
            this->WedgeToTetra(vwedge[4], cellIds, newCellArray);
            this->WedgeToTetra(vwedge[5], cellIds, newCellArray);
            break;
          case 13:
          case 57:
            this->WedgeToTetra(vwedge[6], cellIds, newCellArray);
            this->WedgeToTetra(vwedge[7], cellIds, newCellArray);
            break;
          case 34:
          case 25:
            this->WedgeToTetra(vwedge[8], cellIds, newCellArray);
            this->WedgeToTetra(vwedge[9], cellIds, newCellArray);
            break;
          case 14:
          case 27:
            this->WedgeToTetra(vwedge[10], cellIds, newCellArray);
            this->WedgeToTetra(vwedge[11], cellIds, newCellArray);
            break;
        }
      }
    }
    break;

    case VTK_WEDGE:
    case VTK_QUADRATIC_WEDGE:
    case VTK_QUADRATIC_LINEAR_WEDGE:
    case VTK_BIQUADRATIC_QUADRATIC_WEDGE:
      if (npts == 6) // create 3 tetrahedra
      {
        // first tetrahedron
        constexpr vtkIdType vwedge[6][4] = {
          { 0, 4, 3, 5 },
          { 1, 4, 3, 5 },
          { 2, 4, 3, 5 },
          { 3, 0, 1, 2 },
          { 4, 0, 1, 2 },
          { 5, 0, 1, 2 },
        };
        xmin = cellIds[0];
        id = 0;
        for (i = 1; i < 6; i++)
        {
          if (xmin > cellIds[i])
          {
            xmin = cellIds[i]; // the smallest global index
            id = i;            // local index
          }
        }
        newCellArray->InsertNextCell(4, vwedge[id]);

        // Pyramid :create 2 tetrahedra

        constexpr vtkIdType vert[6][5] = {
          { 1, 2, 5, 4, 0 },
          { 2, 0, 3, 5, 1 },
          { 3, 0, 1, 4, 2 },
          { 1, 2, 5, 4, 3 },
          { 2, 0, 3, 5, 4 },
          { 3, 0, 1, 4, 5 },
        };
        constexpr vtkIdType vpy[8][4] = {
          { 0, 1, 2, 4 },
          { 0, 2, 3, 4 },
          { 1, 2, 3, 4 },
          { 1, 3, 0, 4 },
          { 2, 3, 0, 4 },
          { 2, 0, 1, 4 },
          { 3, 0, 1, 4 },
          { 3, 1, 2, 4 },
        };
        xmin = cellIds[vert[id][0]];
        tabp[0] = vert[id][0];
        idpy = 0;
        for (i = 1; i < 4; i++)
        {
          tabp[i] = vert[id][i];
          if (xmin > cellIds[vert[id][i]])
          {
            xmin = cellIds[vert[id][i]]; // global index
            idpy = i;                    // local index
          }
        }
        tabp[4] = vert[id][4];
        for (j = 0; j < 4; j++)
        {
          tab[j] = tabp[vpy[2 * idpy][j]];
        }
        newCellArray->InsertNextCell(4, tab);

        for (j = 0; j < 4; j++)
        {
          tab[j] = tabp[vpy[2 * idpy + 1][j]];
        }
        newCellArray->InsertNextCell(4, tab);
      }
      else
      {
        vtkErrorMacro(<< " This cell is not a wedge\n");
        return;
      }
      break;

    case VTK_PYRAMID: // Create 2 tetrahedra
    case VTK_QUADRATIC_PYRAMID:
    case VTK_TRIQUADRATIC_PYRAMID:
      if (npts == 5)
      {
        // note: the first element vpyram[][0] is the smallest index of pyramid
        constexpr vtkIdType vpyram[8][4] = {
          { 0, 1, 2, 4 },
          { 0, 2, 3, 4 },
          { 1, 2, 3, 4 },
          { 1, 3, 0, 4 },
          { 2, 3, 0, 4 },
          { 2, 0, 1, 4 },
          { 3, 0, 1, 4 },
          { 3, 1, 2, 4 },
        };
        xmin = cellIds[0];
        id = 0;
        for (i = 1; i < 4; i++)
        {
          if (xmin > cellIds[i])
          {
            xmin = cellIds[i]; // the smallest global index of square face
            id = i;            // local index
          }
        }
        newCellArray->InsertNextCell(4, vpyram[2 * id]);
        newCellArray->InsertNextCell(4, vpyram[2 * id + 1]);
      }
      else
      {
        vtkErrorMacro(<< " This cell is not a pyramid\n");
        return;
      }
      break;
  }
}

//------------------------------------------------------------------------------
// The new cell created in intersection between tetrahedron and plane
// are tetrahedron or wedges or pyramids.
//
// CreateTetra is used to subdivide wedges and pyramids in tetrahedron.  The
// proper vertex ordering for wedges and pyramids can be found in "The
// Visualization Toolkit."  In the third edition, they are in Figure 5-2 on page
// 115 in section 5.4 ("Cell Types") in the "Basic Data Representation" chapter.
//
void vtkBoxClipDataSet::CreateTetra(
  vtkIdType npts, const vtkIdType* cellIds, vtkCellArray* newCellArray)
{
  vtkIdType tabp[5];
  vtkIdType tab[3][4];

  unsigned int i, j;
  unsigned int id = 0;
  unsigned int idpy;

  vtkIdType xmin;

  if (npts == 6)
  {
    // VTK_WEDGE: Create 3 tetrahedra
    // first tetrahedron
    constexpr vtkIdType vwedge[6][4] = {
      { 0, 4, 3, 5 },
      { 1, 4, 3, 5 },
      { 2, 4, 3, 5 },
      { 3, 0, 1, 2 },
      { 4, 0, 1, 2 },
      { 5, 0, 1, 2 },
    };
    xmin = cellIds[0];
    id = 0;
    for (i = 1; i < 6; i++)
    {
      if (xmin > cellIds[i])
      {
        xmin = cellIds[i]; // the smallest global index
        id = i;            // local index
      }
    }
    for (j = 0; j < 4; j++)
    {
      tab[0][j] = cellIds[vwedge[id][j]];
    }
    newCellArray->InsertNextCell(4, tab[0]);

    // Pyramid: create 2 tetrahedra

    constexpr vtkIdType vert[6][5] = {
      { 1, 2, 5, 4, 0 },
      { 2, 0, 3, 5, 1 },
      { 3, 0, 1, 4, 2 },
      { 1, 2, 5, 4, 3 },
      { 2, 0, 3, 5, 4 },
      { 3, 0, 1, 4, 5 },
    };
    constexpr vtkIdType vpy[8][4] = {
      { 0, 1, 2, 4 },
      { 0, 2, 3, 4 },
      { 1, 2, 3, 4 },
      { 1, 3, 0, 4 },
      { 2, 3, 0, 4 },
      { 2, 0, 1, 4 },
      { 3, 0, 1, 4 },
      { 3, 1, 2, 4 },
    };
    xmin = cellIds[vert[id][0]];
    tabp[0] = vert[id][0];
    idpy = 0;
    for (i = 1; i < 4; i++)
    {
      tabp[i] = vert[id][i];
      if (xmin > cellIds[vert[id][i]])
      {
        xmin = cellIds[vert[id][i]]; // global index
        idpy = i;                    // local index
      }
    }
    tabp[4] = vert[id][4];
    for (j = 0; j < 4; j++)
    {
      tab[1][j] = cellIds[tabp[vpy[2 * idpy][j]]];
    }
    newCellArray->InsertNextCell(4, tab[1]);

    for (j = 0; j < 4; j++)
    {
      tab[2][j] = cellIds[tabp[vpy[2 * idpy + 1][j]]];
    }
    newCellArray->InsertNextCell(4, tab[2]);
  }
  else
  {
    // VTK_PYRAMID: Create 2 tetrahedra
    // The first element in each set is the smallest index of pyramid
    constexpr vtkIdType vpyram[8][4] = {
      { 0, 1, 2, 4 },
      { 0, 2, 3, 4 },
      { 1, 2, 3, 4 },
      { 1, 3, 0, 4 },
      { 2, 3, 0, 4 },
      { 2, 0, 1, 4 },
      { 3, 0, 1, 4 },
      { 3, 1, 2, 4 },
    };
    xmin = cellIds[0];
    id = 0;
    for (i = 1; i < 4; i++)
    {
      if (xmin > cellIds[i])
      {
        xmin = cellIds[i]; // the smallest global index of face with 4 vertices
        id = i;            // local index
      }
    }
    for (j = 0; j < 4; j++)
    {
      tab[0][j] = cellIds[vpyram[2 * id][j]];
    }
    newCellArray->InsertNextCell(4, tab[0]);

    for (j = 0; j < 4; j++)
    {
      tab[1][j] = cellIds[vpyram[2 * id + 1][j]];
    }
    newCellArray->InsertNextCell(4, tab[1]);
  }
}

//------------------------------------------------------------------------------
// Clip each cell of an unstructured grid.
//
//------------------------------------------------------------------------------
//(1) How decide when the cell is NOT outside
//
//    Explaining with an example in 2D.
//    Look at 9 region in the picture and the triangle represented there.
//    v0,v1,v2 are vertices of triangle T.
//
//              |         |
//        1     |   2     |     3
//     _ _ _ _ _ _ _ _ _  _ _ _ _ _ ymax
//              |         | v1
//        4     |   5     |/\   6
//              |         /  \             .
//     _ _ _ _ _ _ _ _ _ /| _ \_ _ _ymin
//              |     v2/_|_ _ \v0
//        7     |   8     |     9
//            xmin       xmax
//
// set test={1,1,1,1} (one test for each plane: test={xmin,xmax,ymin,ymax} )
//    for each vertex, if the test is true set 0 in the test table:
//    vO > xmin ?, v0 < xmax ?, v0 > ymin ?, v0 < ymax ?
//    In the example: test={0,1,1,0}
//    v1 > xmin ?, v1 < xmax ?, v1 > ymin ?, v1 < ymax ?
//    In the example: test={0,1,0,0}
//    v2 > xmin ?, v2 < xmax ?, v2 > ymin ?, v2 < ymax ?
//    In the Example: test={0,0,0,0} -> triangle is NOT outside.
//
//
//    In general, look at the possibilities of each region
//
//    (1,0,0,1)  | (0,0,0,1)  | (0,1,0,1)
//    - - - - - - - - - - - - - - - - - -
//    (1,0,0,0)  | (0,0,0,0)  | (0,1,0,0)
//    - - - - - - - - - - - - - - - - - -
//    (1,0,1,0)  | (0,0,1,0)  | (0,1,1,0)
//
//    In the case above we have:(v0,v1,v2)
//    (1,1,1,1)->(0,1,1,0)->(0,1,0,0)->(0,0,0,0)
//
//    If you have one vertex in region 5, this triangle is NOT outside
//    The triangle IS outside if (v0,v1,v2) are in region like
//    (1,2,3): (1,0,0,1)->(0,0,0,1)->(0,0,0,1)
//    (1,4,7): (1,0,0,1)->(1,0,0,0)->(1,0,0,0)
//    (7,8,9): (1,0,1,0)->(0,0,1,0)->(0,0,1,0)
//    (9,6,3): (0,1,1,0)->(0,1,0,0)->(0,1,0,0)
//
//    Note that if the triangle T is on the region 5, T is NOT outside
//    In fact, T is inside
//
//    You can extend this idea to 3D.
//
//    Note: xmin = this->BoundBoxClip[0][0], xmax=  this->BoundBoxClip[0][1],...
//
//------------------------------------------------------------------------------
// (2) Intersection between Tetrahedron and Plane:
//     Description:
//         vertices of tetrahedron {v0,v1,v2,v3}
//         edge e1 : (v0,v1),   edge e2 : (v1,v2)
//         edge e3 : (v2,v0),   edge e4 : (v0,v3)
//         edge e5 : (v1,v3),   edge e6 : (v2,v3)
//
//         intersecting points p0, pi, ...
//
//         Note: The algorithm search intersection
//               points with these edge order.
//
//          v3                        v3
//       e6/|\ e5                     |
//        / | \            e2         |
//     v2/_ |_ \ v1    v2 - - -v1     |e4
//       \  |  /                      |
//      e3\ | /e1                     |
//         \|/                        |
//          v0                        v0
//
//    (a) Intersecting  4 edges: see tab4[]
//                    -> create: 2 wedges
//                    -> 3 cases
//                    ------------------------------------------
//        case 1246:
//                   if (plane intersecting  edges {e1,e2,e4,e6})
//                   then  p0 belongs to e1, p1 belongs to e2,
//                         p2 belongs to e4, p3 belongs to e6.
//
//                                    v3
//                                     |            v3
//                                     |            /
//                v1    v2 - *- -v1    |           * p3
//              /            p1        * p2       /
//             *p0                     |       v2/
//           /                         |
//          v0                        v0
//          (e1)          (e2)       (e4)        (e6)
//
//                   So, there are two wedges:
//                   {p0,v1,p1,p2,v3,p3}  {p2,v0,p0,p3,v2,p1}
//
//                   Note: if e1 and e2 are intersected by plane,
//                         and the plane intersects 4 edges,
//                         the edge e5 could not be intersected
//                         (skew lines, do not create a planar face)
//                         neither the edge e3((e1,e2,e3) is a face)
//
//                    ------------------------------------------
//        case 2345:
//                   if (plane intersecting  edges {e2,e3,e4,e5})
//                   The two wedges are:
//                   {p3,v3,p2,p0,v2,p1}, {p1,v0,p2,p0,v1,p3}
//
//                    ------------------------------------------
//        case 1356:
//                   if (plane intersecting  edges {e1,e3,e5,e6})
//                   The two wedges are:
//                   {p0,v0,p1,p2,v3,p3}, {p0,v1,p2,p1,v2,p3}
//
//                    -----------------------------------------
//
//    (b) Intersecting  3 edges: tab3[]
//                    ->create: 1 tetrahedron + 1 wedge
//                    -> 4 cases
//                    ------------------------------------------
//    case 134:
//                   if (plane intersecting  edges {e1,e3,e4})
//                   then  p0 belongs to e1, p1 belongs to e3,
//                         p2 belongs to e4.
//
//                           v3
//                           |
//                           |
//                     v2   *p2   v1
//                       \   |   /
//                      p1*  |  *p0
//                          \|/
//                          v0
//
//
//                   So, there are:
//                       1 tetrahedron {v0,p0,p2,p1)
//                       and 1 wedge:  {p0,p2,p1,v1,v3,v2}: tab3[0]
//
//                   Note:(a)if e1 and e3 are intersected by plane,
//                            and the plane intersects 3 edges,
//                            the edges e5 could not be intersected,
//                            if so, e6 be intersected too and the
//                            plane intersect 4 edges.
//                        (b) tetrahedron vertices:
//                            Use the first three indices of tab3:
//                            {v0,p(0),p(2),p(1)}
//
//                    ------------------------------------------
//        case 125:
//                   if (plane intersecting  edges {e1,e2,e5})
//         There are :
//                      1 tetrahedron: {v1,p0,p2,p1}
//                      1 wedge      : {p0,p2,p1,v0,v3,v2},
//                    ------------------------------------------
//    case 236:
//                   if (plane intersecting  edges {e2,e3,e6})
//         There are :
//                      1 tetrahedron: {v2,p0,p1,p2}
//                      1 wedge      : {p0,p1,p2,v1,v0,v3},
//                    ------------------------------------------
//        case 456:
//                  if (plane intersecting  edges {e4,e5,e6})
//         There are :
//                      1 tetrahedron: {v3,p0,p1,p2}
//                      1 wedge      : {p0,p1,p2,v0,v1,v2},
//                    ------------------------------------------
//
//    (c) Intersecting  2 edges and 1 vertex: tab2[]
//                      -> create: 1 tetrahedron + 1 pyramid
//      -> 12 cases
//                    ------------------------------------------
//        case 12:        v3        indexes:{0,0,1,2,3},
//                        *                  0 1 2 3 4
//                    e6/ | \ e5
//                     /  |p1\         tetrahedron: indices(*,4,1,2)
//                  v2/_ _|_*_\ v1
//                    \   |   /
//                   e3\  |p0*
//                      \ | /e1
//                        v0
//
//                  if (plane intersecting  edges {e1,e2}
//                      and one vertex)
//           There are
//                      1 tetrahedron: {v1,v3,p0,p1}
//                      1 pyramid    : {v0,p0,p1,v2,v3},
//
//                  Note:  if e1 and e2 are intersected by plane,
//                            and the plane intersects 2 edges,
//                            then the vertex is v3.
//                    ------------------------------------------
//
//        case 13:                   indexes{2,1,0,1,3}
//                      Intersecting edges {e1,e3}
//                      Intersectinf vertes v3,
//
//                      1 tetrahedron: {v0,v3,p1,p0}
//                      1 pyramid    : {v2,p1,p0,v1,v3},
//                    ------------------------------------------
//        all cases see tab2[]
//                    ------------------------------------------
//    (d) Intersecting  1 edges and 2 vertices: tab1[]
//                      -> create: 2 tetrahedra
//                      -> 6 cases
//
//        - edge e1, vertices {v2,v3}(e6)
//
//                        v3
//                        *              tetrahedra: {p0,v2,v3,v1},{p0,v3,v2,v0}
//                    e6/ | \ e5
//                     /  |  \                  .
//                 v2 *_ _|_ _\ v1
//                    \   |   /
//                   e3\  |p0*
//                      \ | /e1
//                        v0
//         - other cases see tab1[]
//------------------------------------------------------------------------------
//
void vtkBoxClipDataSet::ClipBox(vtkPoints* newPoints, vtkGenericCell* cell,
  vtkIncrementalPointLocator* locator, vtkCellArray* tets, vtkPointData* inPD, vtkPointData* outPD,
  vtkCellData* inCD, vtkIdType cellId, vtkCellData* outCD)
{
  vtkIdType cellType = cell->GetCellType();
  vtkIdList* cellIds = cell->GetPointIds();
  vtkCellArray* arraytetra = vtkCellArray::New();
  vtkPoints* cellPts = cell->GetPoints();
  vtkIdType npts = cellPts->GetNumberOfPoints();
  std::vector<vtkIdType> cellptId(npts);
  vtkIdType iid[4];
  const vtkIdType* v_id = nullptr;
  vtkIdType *verts, v1, v2;
  vtkIdType ptId;
  vtkIdType tab_id[6];
  vtkIdType ptstetra = 4;

  vtkIdType i, j;
  unsigned int allInside;
  unsigned int idcellnew;
  unsigned int cutInd;

  /* Edges Tetrahedron */
  vtkIdType edges[6][2] = {
    { 0, 1 },
    { 1, 2 },
    { 2, 0 },
    { 0, 3 },
    { 1, 3 },
    { 2, 3 },
  };
  double value, deltaScalar;
  double t;
  double *p1, *p2;
  double x[3], v[3];
  double v_tetra[4][3];

  for (i = 0; i < npts; i++)
  {
    cellptId[i] = cellIds->GetId(i);
  }

  // Convert all volume cells to tetrahedra
  this->CellGrid(cellType, npts, cellptId.data(), arraytetra);
  unsigned int totalnewtetra = arraytetra->GetNumberOfCells();
  unsigned int idtetranew;

  for (idtetranew = 0; idtetranew < totalnewtetra; idtetranew++)
  {
    arraytetra->GetNextCell(ptstetra, v_id);

    for (allInside = 1, i = 0; i < 4; i++)
    {
      cellPts->GetPoint(v_id[i], v);

      if (!(((v[0] >= this->BoundBoxClip[0][0]) && (v[0] <= this->BoundBoxClip[0][1]) &&
            (v[1] >= this->BoundBoxClip[1][0]) && (v[1] <= this->BoundBoxClip[1][1]) &&
            (v[2] >= this->BoundBoxClip[2][0]) && (v[2] <= this->BoundBoxClip[2][1]))))
      {
        // outside,its type might change later (nearby intersection)
        allInside = 0;
      }
    } // for all points of the tetrahedron.

    // Test Outside: see(1)
    if (!allInside)
    {
      unsigned int test[6] = { 1, 1, 1, 1, 1, 1 };
      for (i = 0; i < 4; i++)
      {
        cellPts->GetPoint(v_id[i], v);

        if (v[0] > this->BoundBoxClip[0][0])
        {
          test[0] = 0;
        }
        if (v[0] < this->BoundBoxClip[0][1])
        {
          test[1] = 0;
        }
        if (v[1] > this->BoundBoxClip[1][0])
        {
          test[2] = 0;
        }
        if (v[1] < this->BoundBoxClip[1][1])
        {
          test[3] = 0;
        }
        if (v[2] > this->BoundBoxClip[2][0])
        {
          test[4] = 0;
        }
        if (v[2] < this->BoundBoxClip[2][1])
        {
          test[5] = 0;
        }

      } // for all points of the cell.

      if ((test[0] == 1) || (test[1] == 1) || (test[2] == 1) || (test[3] == 1) || (test[4] == 1) ||
        (test[5] == 1))
      {
        continue; // Tetrahedron is outside.
      }
    } // if not all inside.

    for (i = 0; i < 4; i++)
    {
      ptId = cellIds->GetId(v_id[i]);
      cellPts->GetPoint(v_id[i], v);
      // Currently all points are injected because of the possibility
      // of intersection point merging.
      if (locator->InsertUniquePoint(v, iid[i]))
      {
        outPD->CopyData(inPD, ptId, iid[i]);
      }

    } // for all points of the tetrahedron.

    if (allInside)
    {
      // Tetrahedron inside.
      vtkIdType newCellId = tets->InsertNextCell(4, iid);
      outCD->CopyData(inCD, cellId, newCellId);
      continue;
    }

    double *pedg1, *pedg2;

    // Tetrahedron Intersection Cases
    constexpr unsigned int tab4[6][6] = {
      { 1, 1, 0, 3, 3, 2 },
      { 2, 0, 0, 3, 2, 1 },
      { 3, 3, 2, 0, 2, 1 },
      { 1, 0, 2, 0, 1, 3 },
      { 0, 0, 1, 2, 3, 3 },
      { 0, 1, 2, 1, 2, 3 },
    };
    constexpr unsigned int tab3[4][6] = {
      { 0, 2, 1, 1, 3, 2 },
      { 0, 1, 2, 0, 2, 3 },
      { 0, 1, 2, 1, 0, 3 },
      { 0, 1, 2, 0, 1, 2 },
    };
    constexpr unsigned int tab2[12][5] = {
      { 0, 0, 1, 2, 3 },
      { 2, 1, 0, 1, 3 },
      { 1, 0, 1, 0, 3 },
      { 2, 0, 1, 3, 0 },
      { 3, 1, 0, 1, 0 },
      { 1, 0, 1, 2, 0 },
      { 3, 1, 0, 2, 1 },
      { 2, 1, 0, 0, 1 },
      { 0, 0, 1, 3, 1 },
      { 1, 0, 1, 3, 2 },
      { 3, 1, 0, 0, 2 },
      { 0, 0, 1, 1, 2 },
    };
    constexpr unsigned int tab1[12][3] = {
      { 2, 3, 1 },
      { 3, 2, 0 },
      { 3, 0, 1 },
      { 0, 3, 2 },
      { 1, 3, 0 },
      { 3, 1, 2 },
      { 2, 1, 0 },
      { 1, 2, 3 },
      { 2, 0, 3 },
      { 0, 2, 1 },
      { 0, 1, 3 },
      { 1, 0, 2 },
    };

    vtkCellArray* cellarray = vtkCellArray::New();
    vtkIdType newCellId = cellarray->InsertNextCell(4, iid);
    unsigned int planes;

    // Test Cell intersection with each plane of box
    for (planes = 0; planes < 6; planes++)
    {
      // The index of the dimension of the cut plane (x == 0, y == 1, z == 2).
      cutInd = planes / 2;

      // The plane is always parallel to unitary cube.
      value = this->BoundBoxClip[cutInd][planes % 2];

      unsigned int totalnewcells = cellarray->GetNumberOfCells();
      vtkCellArray* newcellArray = vtkCellArray::New();
      int edgeNum;

      for (idcellnew = 0; idcellnew < totalnewcells; idcellnew++)
      {
        unsigned int num_inter = 0;
        unsigned int edges_inter = 0;
        unsigned int i0, i1;
        vtkIdType p_id[4];
        cellarray->GetNextCell(npts, v_id);

        newPoints->GetPoint(v_id[0], v_tetra[0]); // coord (x,y,z)
        newPoints->GetPoint(v_id[1], v_tetra[1]);
        newPoints->GetPoint(v_id[2], v_tetra[2]);
        newPoints->GetPoint(v_id[3], v_tetra[3]);

        for (edgeNum = 0; edgeNum < 6; edgeNum++)
        {
          verts = edges[edgeNum];

          p1 = v_tetra[verts[0]];
          p2 = v_tetra[verts[1]];

          if ((p1[cutInd] < value && value < p2[cutInd]) ||
            (p2[cutInd] < value && value < p1[cutInd]))
          {
            deltaScalar = p2[cutInd] - p1[cutInd];

            if (deltaScalar > 0)
            {
              pedg1 = p1;
              pedg2 = p2;
              v1 = verts[0];
              v2 = verts[1];
            }
            else
            {
              pedg1 = p2;
              pedg2 = p1;
              v1 = verts[1];
              v2 = verts[0];
              deltaScalar = -deltaScalar;
            }

            // linear interpolation
            t = (deltaScalar == 0.0 ? 0.0 : (value - pedg1[cutInd]) / deltaScalar);

            for (j = 0; j < 3; j++)
            {
              x[j] = pedg1[j] + t * (pedg2[j] - pedg1[j]);
            }

            // Incorporate point into output and interpolate edge data as necessary
            edges_inter = edges_inter * 10 + (edgeNum + 1);
            if (locator->InsertUniquePoint(x, p_id[num_inter]))
            {
              this->InterpolateEdge(outPD, p_id[num_inter], v_id[v1], v_id[v2], t);
            }
            num_inter++;
          } // if edge intersects value
        }   // for all edges

        if (num_inter == 0)
        {
          unsigned int outside = 0;
          for (i = 0; i < 4; i++)
          {
            if (((v_tetra[i][cutInd] < value) && ((planes % 2) == 0)) ||
              ((v_tetra[i][cutInd] > value) && ((planes % 2) == 1)))

            // If only one vertex is outside, so the tetrahedron is outside
            // because there is not intersection.
            {
              outside = 1;
              break;
            }
          }
          if (outside == 0)
          {
            // else it is possible intersection if other plane
            newCellId = newcellArray->InsertNextCell(4, v_id);
          }
          continue;
        }
        switch (num_inter)
        {
          case 4: // We have two wedges
            switch (edges_inter)
            {
              case 1246:
                i0 = 0;
                break;
              case 2345:
                i0 = 2;
                break;
              case 1356:
                i0 = 4;
                break;
              default:
                vtkErrorMacro(<< "Intersection not found: Num_inter = " << num_inter
                              << " Edges_inter = " << edges_inter);
                continue;
            }

            if (((v_tetra[3][cutInd] < value) && ((planes % 2) == 0)) ||
              ((v_tetra[3][cutInd] > value) && ((planes % 2) == 1)))
            {

              // The v_tetra[3] is outside box, so the first wedge is outside

              tab_id[0] = p_id[tab4[i0 + 1][0]];
              // ps: v_tetra[3] is always in first wedge (see tab)

              tab_id[1] = v_id[tab4[i0 + 1][1]];
              tab_id[2] = p_id[tab4[i0 + 1][2]];
              tab_id[3] = p_id[tab4[i0 + 1][3]];
              tab_id[4] = v_id[tab4[i0 + 1][4]];
              tab_id[5] = p_id[tab4[i0 + 1][5]];
              this->CreateTetra(6, tab_id, newcellArray);
            }
            else
            {
              tab_id[0] = p_id[tab4[i0][0]];
              tab_id[1] = v_id[tab4[i0][1]];
              tab_id[2] = p_id[tab4[i0][2]];
              tab_id[3] = p_id[tab4[i0][3]];
              tab_id[4] = v_id[tab4[i0][4]];
              tab_id[5] = p_id[tab4[i0][5]];
              this->CreateTetra(6, tab_id, newcellArray);
            }
            break;
          case 3: // We have one tetrahedron and one wedge
            // i0 gets the vertex on the tetrahedron.
            switch (edges_inter)
            {
              case 134:
                i0 = 0;
                break;
              case 125:
                i0 = 1;
                break;
              case 236:
                i0 = 2;
                break;
              case 456:
                i0 = 3;
                break;
              default:
                vtkErrorMacro(<< "Intersection not found: Num_inter = " << num_inter
                              << " Edges_inter = " << edges_inter);
                continue;
            }

            if (((v_tetra[i0][cutInd] < value) && ((planes % 2) == 0)) ||
              ((v_tetra[i0][cutInd] > value) && ((planes % 2) == 1)))
            {

              // Isolate vertex is outside box, so the tetrahedron is outside
              tab_id[0] = p_id[tab3[i0][0]];
              tab_id[1] = p_id[tab3[i0][1]];
              tab_id[2] = p_id[tab3[i0][2]];
              tab_id[3] = v_id[tab3[i0][3]];
              tab_id[4] = v_id[tab3[i0][4]];
              tab_id[5] = v_id[tab3[i0][5]];
              this->CreateTetra(6, tab_id, newcellArray);
            }
            else
            {
              tab_id[0] = p_id[tab3[i0][0]];
              tab_id[1] = p_id[tab3[i0][1]];
              tab_id[2] = p_id[tab3[i0][2]];
              tab_id[3] = v_id[i0];
              newCellId = newcellArray->InsertNextCell(4, tab_id);
            }
            break;

          case 2:                // We have one tetrahedron and one pyramid
            switch (edges_inter) // i1 = vertex of the tetrahedron
            {
              case 12:
                i0 = 0;
                i1 = 1;
                break;
              case 13:
                i0 = 1;
                i1 = 0;
                break;
              case 23:
                i0 = 2;
                i1 = 2;
                break;
              case 25:
                i0 = 3;
                i1 = 1;
                break;
              case 26:
                i0 = 4;
                i1 = 2;
                break;
              case 56:
                i0 = 5;
                i1 = 3;
                break;
              case 34:
                i0 = 6;
                i1 = 0;
                break;
              case 46:
                i0 = 7;
                i1 = 3;
                break;
              case 36:
                i0 = 8;
                i1 = 2;
                break;
              case 14:
                i0 = 9;
                i1 = 0;
                break;
              case 15:
                i0 = 10;
                i1 = 1;
                break;
              case 45:
                i0 = 11;
                i1 = 3;
                break;
              default:
                vtkErrorMacro(<< "Intersection not found: Num_inter = " << num_inter
                              << " Edges_inter = " << edges_inter);
                continue;
            }
            if (((v_tetra[i1][cutInd] < value) && ((planes % 2) == 0)) ||
              ((v_tetra[i1][cutInd] > value) && ((planes % 2) == 1)))
            {
              // Isolate vertex is outside box, so the tetrahedron is outside
              tab_id[0] = v_id[tab2[i0][0]];
              tab_id[1] = p_id[tab2[i0][1]];
              tab_id[2] = p_id[tab2[i0][2]];
              tab_id[3] = v_id[tab2[i0][3]];
              tab_id[4] = v_id[tab2[i0][4]];
              this->CreateTetra(5, tab_id, newcellArray);
            }
            else
            {
              tab_id[0] = v_id[i1];
              tab_id[1] = v_id[tab2[i0][4]];
              tab_id[2] = p_id[tab2[i0][2]];
              tab_id[3] = p_id[tab2[i0][1]];
              newCellId = newcellArray->InsertNextCell(4, tab_id);
            }
            break;

          case 1: // We have two tetrahedron.
            if ((edges_inter > 6) || (edges_inter < 1))
            {
              vtkErrorMacro(<< "Intersection not found: Num_inter = " << num_inter
                            << " Edges_inter = " << edges_inter);
              continue;
            }
            if (((v_tetra[tab1[2 * edges_inter - 1][2]][cutInd] < value) && ((planes % 2) == 0)) ||
              ((v_tetra[tab1[2 * edges_inter - 1][2]][cutInd] > value) && ((planes % 2) == 1)))
            {
              // Isolate vertex is outside box, so the tetrahedron is outside
              tab_id[0] = p_id[0];
              tab_id[1] = v_id[tab1[2 * edges_inter - 2][0]];
              tab_id[2] = v_id[tab1[2 * edges_inter - 2][1]];
              tab_id[3] = v_id[tab1[2 * edges_inter - 2][2]];
              newCellId = newcellArray->InsertNextCell(4, tab_id);
            }
            else
            {
              tab_id[0] = p_id[0];
              tab_id[1] = v_id[tab1[2 * edges_inter - 1][0]];
              tab_id[2] = v_id[tab1[2 * edges_inter - 1][1]];
              tab_id[3] = v_id[tab1[2 * edges_inter - 1][2]];
              newCellId = newcellArray->InsertNextCell(4, tab_id);
            }
            break;
        }
      } // for all new cells

      cellarray->Delete();
      cellarray = newcellArray;
    } // for all planes

    unsigned int totalnewcells = cellarray->GetNumberOfCells();

    for (idcellnew = 0; idcellnew < totalnewcells; idcellnew++)
    {
      cellarray->GetNextCell(npts, v_id);
      newCellId = tets->InsertNextCell(npts, v_id);
      outCD->CopyData(inCD, cellId, newCellId);
    }
    cellarray->Delete();
  }
  arraytetra->Delete();
}

//------------------------------------------------------------------------------
// ClipHexahedron: Box is like hexahedron.
//
// The difference between ClipBox and ClipHexahedron is the outside test.
// The ClipHexahedron use plane equation to decide who is outside.
//
void vtkBoxClipDataSet::ClipHexahedron(vtkPoints* newPoints, vtkGenericCell* cell,
  vtkIncrementalPointLocator* locator, vtkCellArray* tets, vtkPointData* inPD, vtkPointData* outPD,
  vtkCellData* inCD, vtkIdType cellId, vtkCellData* outCD)
{
  vtkIdType idcellnew;
  vtkIdType cellType = cell->GetCellType();
  vtkIdList* cellIds = cell->GetPointIds();
  vtkCellArray* arraytetra = vtkCellArray::New();
  vtkPoints* cellPts = cell->GetPoints();
  vtkIdType npts = cellPts->GetNumberOfPoints();
  std::vector<vtkIdType> cellptId(npts);
  vtkIdType iid[4];
  const vtkIdType* v_id = nullptr;
  vtkIdType *verts, v1, v2;
  vtkIdType ptId;
  vtkIdType tab_id[6];
  vtkIdType ptstetra = 4;

  vtkIdType i, j;
  unsigned int allInside, k;
  unsigned int planes;

  /* Edges Tetrahedron */
  vtkIdType edges[6][2] = {
    { 0, 1 },
    { 1, 2 },
    { 2, 0 },
    { 0, 3 },
    { 1, 3 },
    { 2, 3 },
  };
  double deltaScalar;
  double t;
  double *p1, *p2;
  double v[3], x[3];
  double p[6];
  double v_tetra[4][3];

  for (i = 0; i < npts; i++)
  {
    cellptId[i] = cellIds->GetId(i);
  }

  this->CellGrid(
    cellType, npts, cellptId.data(), arraytetra); // Convert all volume cells to tetrahedra

  unsigned int totalnewtetra = arraytetra->GetNumberOfCells();
  unsigned int idtetranew;

  for (idtetranew = 0; idtetranew < totalnewtetra; idtetranew++)
  {
    arraytetra->GetNextCell(ptstetra, v_id);

    for (allInside = 1, i = 0; i < 4; i++)
    {
      cellPts->GetPoint(v_id[i], v);

      for (k = 0; k < 6; k++)
      {
        p[k] = this->PlaneNormal[k][0] * (v[0] - this->PlanePoint[k][0]) +
          this->PlaneNormal[k][1] * (v[1] - this->PlanePoint[k][1]) +
          this->PlaneNormal[k][2] * (v[2] - this->PlanePoint[k][2]);
      }

      if (!((p[0] <= 0) && (p[1] <= 0) && (p[2] <= 0) && (p[3] <= 0) && (p[4] <= 0) && (p[5] <= 0)))
      {
        allInside = 0;
      }
    } // for all points of the cell.

    // Test Outside
    unsigned int test[6] = { 1, 1, 1, 1, 1, 1 };
    for (i = 0; i < 4; i++)
    {
      cellPts->GetPoint(v_id[i], v);

      // Use plane equation
      for (k = 0; k < 6; k++)
      {
        p[k] = this->PlaneNormal[k][0] * (v[0] - this->PlanePoint[k][0]) +
          this->PlaneNormal[k][1] * (v[1] - this->PlanePoint[k][1]) +
          this->PlaneNormal[k][2] * (v[2] - this->PlanePoint[k][2]);
      }

      for (k = 0; k < 3; k++)
      {
        if (p[2 * k] < 0)
        {
          test[2 * k] = 0;
        }
        if (p[2 * k + 1] < 0)
        {
          test[2 * k + 1] = 0;
        }
      }

    } // for all points of the cell.

    if (!allInside &&
      ((test[0] == 1) || (test[1] == 1) || (test[2] == 1) || (test[3] == 1) || (test[4] == 1) ||
        (test[5] == 1)))
    {
      continue; // Tetrahedron is outside.
    }

    for (i = 0; i < 4; i++)
    {
      ptId = cellIds->GetId(v_id[i]);
      cellPts->GetPoint(v_id[i], v);

      // Currently all points are injected because of the possibility
      // of intersection point merging.
      if (locator->InsertUniquePoint(v, iid[i]))
      {
        outPD->CopyData(inPD, ptId, iid[i]);
      }
    } // for all points of the tetrahedron.

    if (allInside)
    {
      vtkIdType newCellId = tets->InsertNextCell(4, iid); // Tetrahedron inside.
      outCD->CopyData(inCD, cellId, newCellId);
      continue;
    }

    double *pedg1, *pedg2;

    // Tetrahedron Intersection Cases
    constexpr unsigned int tab4[6][6] = {
      { 1, 1, 0, 3, 3, 2 },
      { 2, 0, 0, 3, 2, 1 },
      { 3, 3, 2, 0, 2, 1 },
      { 1, 0, 2, 0, 1, 3 },
      { 0, 0, 1, 2, 3, 3 },
      { 0, 1, 2, 1, 2, 3 },
    };
    constexpr unsigned int tab3[4][6] = {
      { 0, 2, 1, 1, 3, 2 },
      { 0, 1, 2, 0, 2, 3 },
      { 0, 1, 2, 1, 0, 3 },
      { 0, 1, 2, 0, 1, 2 },
    };
    constexpr unsigned int tab2[12][5] = {
      { 0, 0, 1, 2, 3 },
      { 2, 1, 0, 1, 3 },
      { 1, 0, 1, 0, 3 },
      { 2, 0, 1, 3, 0 },
      { 3, 1, 0, 1, 0 },
      { 1, 0, 1, 2, 0 },
      { 3, 1, 0, 2, 1 },
      { 2, 1, 0, 0, 1 },
      { 0, 0, 1, 3, 1 },
      { 1, 0, 1, 3, 2 },
      { 3, 1, 0, 0, 2 },
      { 0, 0, 1, 1, 2 },
    };
    constexpr unsigned int tab1[12][3] = {
      { 2, 3, 1 },
      { 3, 2, 0 },
      { 3, 0, 1 },
      { 0, 3, 2 },
      { 1, 3, 0 },
      { 3, 1, 2 },
      { 2, 1, 0 },
      { 1, 2, 3 },
      { 2, 0, 3 },
      { 0, 2, 1 },
      { 0, 1, 3 },
      { 1, 0, 2 },
    };

    vtkCellArray* cellarray = vtkCellArray::New();
    vtkIdType newCellId = cellarray->InsertNextCell(4, iid);

    // Test Cell intersection with each plane of box
    for (planes = 0; planes < 6; planes++)
    {
      vtkIdType totalnewcells = cellarray->GetNumberOfCells();
      vtkCellArray* newcellArray = vtkCellArray::New();

      for (idcellnew = 0; idcellnew < totalnewcells; idcellnew++)
      {
        unsigned int i0, i1;
        unsigned int num_inter = 0;
        unsigned int edges_inter = 0;
        vtkIdType p_id[4];

        cellarray->GetNextCell(npts, v_id);

        newPoints->GetPoint(v_id[0], v_tetra[0]); // coord (x,y,z)
        newPoints->GetPoint(v_id[1], v_tetra[1]);
        newPoints->GetPoint(v_id[2], v_tetra[2]);
        newPoints->GetPoint(v_id[3], v_tetra[3]);

        p[0] = this->PlaneNormal[planes][0] * (v_tetra[0][0] - this->PlanePoint[planes][0]) +
          this->PlaneNormal[planes][1] * (v_tetra[0][1] - this->PlanePoint[planes][1]) +
          this->PlaneNormal[planes][2] * (v_tetra[0][2] - this->PlanePoint[planes][2]);
        p[1] = this->PlaneNormal[planes][0] * (v_tetra[1][0] - this->PlanePoint[planes][0]) +
          this->PlaneNormal[planes][1] * (v_tetra[1][1] - this->PlanePoint[planes][1]) +
          this->PlaneNormal[planes][2] * (v_tetra[1][2] - this->PlanePoint[planes][2]);
        p[2] = this->PlaneNormal[planes][0] * (v_tetra[2][0] - this->PlanePoint[planes][0]) +
          this->PlaneNormal[planes][1] * (v_tetra[2][1] - this->PlanePoint[planes][1]) +
          this->PlaneNormal[planes][2] * (v_tetra[2][2] - this->PlanePoint[planes][2]);
        p[3] = this->PlaneNormal[planes][0] * (v_tetra[3][0] - this->PlanePoint[planes][0]) +
          this->PlaneNormal[planes][1] * (v_tetra[3][1] - this->PlanePoint[planes][1]) +
          this->PlaneNormal[planes][2] * (v_tetra[3][2] - this->PlanePoint[planes][2]);

        for (int edgeNum = 0; edgeNum < 6; edgeNum++)
        {
          verts = edges[edgeNum];

          p1 = v_tetra[verts[0]];
          p2 = v_tetra[verts[1]];
          double s1 = p[verts[0]];
          double s2 = p[verts[1]];
          if ((s1 * s2) < 0)
          {
            deltaScalar = s2 - s1;

            if (deltaScalar > 0)
            {
              pedg1 = p1;
              pedg2 = p2;
              v1 = verts[0];
              v2 = verts[1];
            }
            else
            {
              pedg1 = p2;
              pedg2 = p1;
              v1 = verts[1];
              v2 = verts[0];
              deltaScalar = -deltaScalar;
              t = s1;
              s1 = s2;
              s2 = t;
            }

            // linear interpolation
            t = (deltaScalar == 0.0 ? 0.0 : (-s1) / deltaScalar);

            for (j = 0; j < 3; j++)
            {
              x[j] = pedg1[j] + t * (pedg2[j] - pedg1[j]);
            }

            // Incorporate point into output and interpolate edge data as necessary
            edges_inter = edges_inter * 10 + (edgeNum + 1);

            if (locator->InsertUniquePoint(x, p_id[num_inter]))
            {
              this->InterpolateEdge(outPD, p_id[num_inter], v_id[v1], v_id[v2], t);
            }

            num_inter++;
          } // if edge intersects value
        }   // for all edges
        if (num_inter == 0)
        {
          unsigned int outside = 0;
          for (i = 0; i < 4; i++)
          {
            if (p[i] > 0)
            {
              // If only one vertex is outside, so the tetrahedron is outside
              // because there is not intersection.
              // some vertex could be on plane, so you need to test all vertex

              outside = 1;
              break;
            }
          }
          if (outside == 0)
          {
            // else it is possible intersection if other plane

            newCellId = newcellArray->InsertNextCell(4, v_id);
          }
          continue;
        }
        switch (num_inter)
        {
          case 4: // We have two wedges
            switch (edges_inter)
            {
              case 1246:
                i0 = 0;
                break;
              case 2345:
                i0 = 2;
                break;
              case 1356:
                i0 = 4;
                break;
              default:
                vtkErrorMacro(<< "Intersection not found: Num_inter = " << num_inter
                              << " Edges_inter = " << edges_inter);
                continue;
            }

            if (p[3] > 0)
            {
              // The v_tetra[3] is outside box, so the first wedge is outside
              // ps: v_tetra[3] is always in first wedge (see tab)

              tab_id[0] = p_id[tab4[i0 + 1][0]];
              tab_id[1] = v_id[tab4[i0 + 1][1]];
              tab_id[2] = p_id[tab4[i0 + 1][2]];
              tab_id[3] = p_id[tab4[i0 + 1][3]];
              tab_id[4] = v_id[tab4[i0 + 1][4]];
              tab_id[5] = p_id[tab4[i0 + 1][5]];
              this->CreateTetra(6, tab_id, newcellArray);
            }
            else
            {
              tab_id[0] = p_id[tab4[i0][0]];
              tab_id[1] = v_id[tab4[i0][1]];
              tab_id[2] = p_id[tab4[i0][2]];
              tab_id[3] = p_id[tab4[i0][3]];
              tab_id[4] = v_id[tab4[i0][4]];
              tab_id[5] = p_id[tab4[i0][5]];
              this->CreateTetra(6, tab_id, newcellArray);
            }
            break;
          case 3: // We have one tetrahedron and one wedge
            switch (edges_inter)
            {
              case 134:
                i0 = 0;
                break;
              case 125:
                i0 = 1;
                break;
              case 236:
                i0 = 2;
                break;
              case 456:
                i0 = 3;
                break;
              default:
                vtkErrorMacro(<< "Intersection not found: Num_inter = " << num_inter
                              << " Edges_inter = " << edges_inter);
                continue;
            }

            if (p[i0] > 0)
            {
              // Isolate vertex is outside box, so the tetrahedron is outside

              tab_id[0] = p_id[tab3[i0][0]];
              tab_id[1] = p_id[tab3[i0][1]];
              tab_id[2] = p_id[tab3[i0][2]];
              tab_id[3] = v_id[tab3[i0][3]];
              tab_id[4] = v_id[tab3[i0][4]];
              tab_id[5] = v_id[tab3[i0][5]];
              this->CreateTetra(6, tab_id, newcellArray);
            }
            else
            {
              tab_id[0] = p_id[tab3[i0][0]];
              tab_id[1] = p_id[tab3[i0][1]];
              tab_id[2] = p_id[tab3[i0][2]];
              tab_id[3] = v_id[i0];
              newCellId = newcellArray->InsertNextCell(4, tab_id);
            }
            break;
          case 2:                // We have one tetrahedron and one pyramid
            switch (edges_inter) // i1 = vertex of the tetrahedron
            {
              case 12:
                i0 = 0;
                i1 = 1;
                break;
              case 13:
                i0 = 1;
                i1 = 0;
                break;
              case 23:
                i0 = 2;
                i1 = 2;
                break;
              case 25:
                i0 = 3;
                i1 = 1;
                break;
              case 26:
                i0 = 4;
                i1 = 2;
                break;
              case 56:
                i0 = 5;
                i1 = 3;
                break;
              case 34:
                i0 = 6;
                i1 = 0;
                break;
              case 46:
                i0 = 7;
                i1 = 3;
                break;
              case 36:
                i0 = 8;
                i1 = 2;
                break;
              case 14:
                i0 = 9;
                i1 = 0;
                break;
              case 15:
                i0 = 10;
                i1 = 1;
                break;
              case 45:
                i0 = 11;
                i1 = 3;
                break;
              default:
                vtkErrorMacro(<< "Intersection not found: Num_inter = " << num_inter
                              << " Edges_inter = " << edges_inter);
                continue;
            }
            if (p[i1] > 0)
            {
              // Isolate vertex is outside box, so the tetrahedron is outside

              tab_id[0] = v_id[tab2[i0][0]];
              tab_id[1] = p_id[tab2[i0][1]];
              tab_id[2] = p_id[tab2[i0][2]];
              tab_id[3] = v_id[tab2[i0][3]];
              tab_id[4] = v_id[tab2[i0][4]];
              this->CreateTetra(5, tab_id, newcellArray);
            }
            else
            {
              tab_id[0] = v_id[i1];
              tab_id[1] = v_id[tab2[i0][4]];
              tab_id[2] = p_id[tab2[i0][2]];
              tab_id[3] = p_id[tab2[i0][1]];
              newCellId = newcellArray->InsertNextCell(4, tab_id);
            }
            break;
          case 1: // We have two tetrahedron.
            if ((edges_inter > 6) || (edges_inter < 1))
            {
              vtkErrorMacro(<< "Intersection not found: Num_inter = " << num_inter
                            << " Edges_inter = " << edges_inter);
              continue;
            }
            if (p[tab1[2 * edges_inter - 1][2]] > 0)
            {
              // Isolate vertex is outside box, so the tetrahedron is outside
              tab_id[0] = p_id[0];
              tab_id[1] = v_id[tab1[2 * edges_inter - 2][0]];
              tab_id[2] = v_id[tab1[2 * edges_inter - 2][1]];
              tab_id[3] = v_id[tab1[2 * edges_inter - 2][2]];
              newCellId = newcellArray->InsertNextCell(4, tab_id);
            }
            else
            {
              tab_id[0] = p_id[0];
              tab_id[1] = v_id[tab1[2 * edges_inter - 1][0]];
              tab_id[2] = v_id[tab1[2 * edges_inter - 1][1]];
              tab_id[3] = v_id[tab1[2 * edges_inter - 1][2]];
              newCellId = newcellArray->InsertNextCell(4, tab_id);
            }
            break;
        }
      } // for all new cells
      cellarray->Delete();
      cellarray = newcellArray;
    } // for all planes

    vtkIdType totalnewcells = cellarray->GetNumberOfCells();

    for (idcellnew = 0; idcellnew < totalnewcells; idcellnew++)
    {
      cellarray->GetNextCell(npts, v_id);
      newCellId = tets->InsertNextCell(npts, v_id);
      outCD->CopyData(inCD, cellId, newCellId);
    }
    cellarray->Delete();
  }
  arraytetra->Delete();
}
//------------------------------------------------------------------------------
// ClipBoxInOut
//
// The difference between ClipBox and ClipBoxInOut is the outputs.
// The ClipBoxInOut generate both outputs: inside and outside the clip box.
//
void vtkBoxClipDataSet::ClipBoxInOut(vtkPoints* newPoints, vtkGenericCell* cell,
  vtkIncrementalPointLocator* locator, vtkCellArray** tets, vtkPointData* inPD,
  vtkPointData** outPD, vtkCellData* inCD, vtkIdType cellId, vtkCellData** outCD)
{
  vtkIdType cellType = cell->GetCellType();
  vtkIdList* cellIds = cell->GetPointIds();
  vtkCellArray* arraytetra = vtkCellArray::New();
  vtkPoints* cellPts = cell->GetPoints();
  vtkIdType npts = cellPts->GetNumberOfPoints();
  std::vector<vtkIdType> cellptId(npts);
  vtkIdType iid[4];
  const vtkIdType* v_id = nullptr;
  vtkIdType *verts, v1, v2;
  vtkIdType ptId;
  vtkIdType ptIdout[4];
  vtkIdType tab_id[6];
  vtkIdType ptstetra = 4;

  int i, j;
  int allInside;
  int cutInd;

  unsigned int planes;
  unsigned int idcellnew;
  unsigned int idtetranew;

  /* Edges Tetrahedron */
  vtkIdType edges[6][2] = {
    { 0, 1 },
    { 1, 2 },
    { 2, 0 },
    { 0, 3 },
    { 1, 3 },
    { 2, 3 },
  };
  double value, deltaScalar;
  double t;
  double v[3], x[3];
  double v_tetra[4][3];
  double *p1, *p2;

  for (i = 0; i < npts; i++)
  {
    cellptId[i] = cellIds->GetId(i);
  }

  // Convert all volume cells to tetrahedra
  this->CellGrid(cellType, npts, cellptId.data(), arraytetra);
  unsigned int totalnewtetra = arraytetra->GetNumberOfCells();

  for (idtetranew = 0; idtetranew < totalnewtetra; idtetranew++)
  {
    arraytetra->GetNextCell(ptstetra, v_id);

    for (allInside = 1, i = 0; i < 4; i++)
    {
      cellPts->GetPoint(v_id[i], v);

      if (!(((v[0] >= this->BoundBoxClip[0][0]) && (v[0] <= this->BoundBoxClip[0][1]) &&
            (v[1] >= this->BoundBoxClip[1][0]) && (v[1] <= this->BoundBoxClip[1][1]) &&
            (v[2] >= this->BoundBoxClip[2][0]) && (v[2] <= this->BoundBoxClip[2][1]))))
      {
        // outside,its type might change later (nearby intersection)
        allInside = 0;
      }
    } // for all points of the cell.

    // Test Outside: see(1)
    if (!allInside)
    {
      unsigned int test[6] = { 1, 1, 1, 1, 1, 1 };
      for (i = 0; i < 4; i++)
      {
        ptIdout[i] = cellIds->GetId(v_id[i]);
        cellPts->GetPoint(v_id[i], v_tetra[i]);

        if (v_tetra[i][0] > this->BoundBoxClip[0][0])
        {
          test[0] = 0;
        }
        if (v_tetra[i][0] < this->BoundBoxClip[0][1])
        {
          test[1] = 0;
        }
        if (v_tetra[i][1] > this->BoundBoxClip[1][0])
        {
          test[2] = 0;
        }
        if (v_tetra[i][1] < this->BoundBoxClip[1][1])
        {
          test[3] = 0;
        }
        if (v_tetra[i][2] > this->BoundBoxClip[2][0])
        {
          test[4] = 0;
        }
        if (v_tetra[i][2] < this->BoundBoxClip[2][1])
        {
          test[5] = 0;
        }
      } // for all points of the cell.

      if ((test[0] == 1) || (test[1] == 1) || (test[2] == 1) || (test[3] == 1) || (test[4] == 1) ||
        (test[5] == 1))
      {
        for (i = 0; i < 4; i++)
        {
          if (locator->InsertUniquePoint(v_tetra[i], iid[i]))
          {
            outPD[0]->CopyData(inPD, ptIdout[i], iid[i]);
            outPD[1]->CopyData(inPD, ptIdout[i], iid[i]);
          }
        }
        int newCellId = tets[1]->InsertNextCell(4, iid);
        outCD[1]->CopyData(inCD, cellId, newCellId);
        continue; // Tetrahedron is outside.
      }
    } // if not allinside.

    for (i = 0; i < 4; i++)
    {
      ptId = cellIds->GetId(v_id[i]);
      cellPts->GetPoint(v_id[i], v);

      // Currently all points are injected because of the possibility
      // of intersection point merging.
      if (locator->InsertUniquePoint(v, iid[i]))
      {
        outPD[0]->CopyData(inPD, ptId, iid[i]);
        outPD[1]->CopyData(inPD, ptId, iid[i]);
      }

    } // for all points of the tetrahedron.

    if (allInside)
    {
      // Tetrahedron inside.
      int newCellId = tets[0]->InsertNextCell(4, iid);
      outCD[0]->CopyData(inCD, cellId, newCellId);
      continue;
    }

    double *pedg1, *pedg2;

    // Tetrahedron Intersection Cases
    constexpr unsigned int tab4[6][6] = {
      { 1, 1, 0, 3, 3, 2 },
      { 2, 0, 0, 3, 2, 1 },
      { 3, 3, 2, 0, 2, 1 },
      { 1, 0, 2, 0, 1, 3 },
      { 0, 0, 1, 2, 3, 3 },
      { 0, 1, 2, 1, 2, 3 },
    };
    constexpr unsigned int tab3[4][6] = {
      { 0, 2, 1, 1, 3, 2 },
      { 0, 1, 2, 0, 2, 3 },
      { 0, 1, 2, 1, 0, 3 },
      { 0, 1, 2, 0, 1, 2 },
    };
    constexpr unsigned int tab2[12][5] = {
      { 0, 0, 1, 2, 3 },
      { 2, 1, 0, 1, 3 },
      { 1, 0, 1, 0, 3 },
      { 2, 0, 1, 3, 0 },
      { 3, 1, 0, 1, 0 },
      { 1, 0, 1, 2, 0 },
      { 3, 1, 0, 2, 1 },
      { 2, 1, 0, 0, 1 },
      { 0, 0, 1, 3, 1 },
      { 1, 0, 1, 3, 2 },
      { 3, 1, 0, 0, 2 },
      { 0, 0, 1, 1, 2 },
    };
    constexpr unsigned int tab1[12][3] = {
      { 2, 3, 1 },
      { 3, 2, 0 },
      { 3, 0, 1 },
      { 0, 3, 2 },
      { 1, 3, 0 },
      { 3, 1, 2 },
      { 2, 1, 0 },
      { 1, 2, 3 },
      { 2, 0, 3 },
      { 0, 2, 1 },
      { 0, 1, 3 },
      { 1, 0, 2 },
    };

    vtkCellArray* cellarray = vtkCellArray::New();
    vtkCellArray* cellarrayout = vtkCellArray::New();
    int newCellId = cellarray->InsertNextCell(4, iid);

    // Test Cell intersection with each plane of box
    for (planes = 0; planes < 6; planes++)
    {
      // The index of the dimension of the cut plane (x == 0, y == 1, z == 2).
      cutInd = planes / 2;

      value =
        this->BoundBoxClip[cutInd][planes % 2]; // The plane is always parallel to unitary cube.

      unsigned int totalnewcells = cellarray->GetNumberOfCells();
      vtkCellArray* newcellArray = vtkCellArray::New();

      for (idcellnew = 0; idcellnew < totalnewcells; idcellnew++)
      {
        unsigned int num_inter = 0;
        unsigned int edges_inter = 0;
        unsigned int i0, i1;
        vtkIdType p_id[4];
        cellarray->GetNextCell(npts, v_id);

        newPoints->GetPoint(v_id[0], v_tetra[0]); // coord (x,y,z)
        newPoints->GetPoint(v_id[1], v_tetra[1]);
        newPoints->GetPoint(v_id[2], v_tetra[2]);
        newPoints->GetPoint(v_id[3], v_tetra[3]);
        for (int edgeNum = 0; edgeNum < 6; edgeNum++)
        {
          verts = edges[edgeNum];

          p1 = v_tetra[verts[0]];
          p2 = v_tetra[verts[1]];

          if ((p1[cutInd] < value && value < p2[cutInd]) ||
            (p2[cutInd] < value && value < p1[cutInd]))
          {
            deltaScalar = p2[cutInd] - p1[cutInd];

            if (deltaScalar > 0)
            {
              pedg1 = p1;
              pedg2 = p2;
              v1 = verts[0];
              v2 = verts[1];
            }
            else
            {
              pedg1 = p2;
              pedg2 = p1;
              v1 = verts[1];
              v2 = verts[0];
              deltaScalar = -deltaScalar;
            }

            // linear interpolation
            t = (deltaScalar == 0.0 ? 0.0 : (value - pedg1[cutInd]) / deltaScalar);

            for (j = 0; j < 3; j++)
            {
              x[j] = pedg1[j] + t * (pedg2[j] - pedg1[j]);
            }

            // Incorporate point into output and interpolate edge data as necessary
            edges_inter = edges_inter * 10 + (edgeNum + 1);
            if (locator->InsertUniquePoint(x, p_id[num_inter]))
            {
              this->InterpolateEdge(outPD[0], p_id[num_inter], v_id[v1], v_id[v2], t);
              this->InterpolateEdge(outPD[1], p_id[num_inter], v_id[v1], v_id[v2], t);
            }

            num_inter++;
          } // if edge intersects value
        }   // for all edges

        if (num_inter == 0)
        {
          unsigned int outside = 0;
          for (i = 0; i < 4; i++)
          {
            if (((v_tetra[i][cutInd] < value) && ((planes % 2) == 0)) ||
              ((v_tetra[i][cutInd] > value) && ((planes % 2) == 1)))
            {
              // If only one vertex is outside, so the tetrahedron is outside
              // because there is not intersection.
              outside = 1;
              break;
            }
          }
          if (outside == 0)
          {
            // else it is possible intersection if other plane
            newCellId = newcellArray->InsertNextCell(4, v_id);
          }
          else
          {
            newCellId = tets[1]->InsertNextCell(4, v_id);
            outCD[1]->CopyData(inCD, cellId, newCellId);
          }
          continue;
        }
        switch (num_inter)
        {
          case 4: // We have two wedges
            switch (edges_inter)
            {
              case 1246:
                i0 = 0;
                break;
              case 2345:
                i0 = 2;
                break;
              case 1356:
                i0 = 4;
                break;
              default:
                vtkErrorMacro(<< "Intersection not found: Num_inter = " << num_inter
                              << " Edges_inter = " << edges_inter);
                continue;
            }
            if (((v_tetra[3][cutInd] < value) && ((planes % 2) == 0)) ||
              ((v_tetra[3][cutInd] > value) && ((planes % 2) == 1)))
            {
              // The v_tetra[3] is outside box, so
              // the first wedge is outside
              // ps: v_tetra[3] is always in first wedge (see tab)

              tab_id[0] = p_id[tab4[i0 + 1][0]]; // Inside
              tab_id[1] = v_id[tab4[i0 + 1][1]];
              tab_id[2] = p_id[tab4[i0 + 1][2]];
              tab_id[3] = p_id[tab4[i0 + 1][3]];
              tab_id[4] = v_id[tab4[i0 + 1][4]];
              tab_id[5] = p_id[tab4[i0 + 1][5]];
              this->CreateTetra(6, tab_id, newcellArray);

              tab_id[0] = p_id[tab4[i0][0]]; // Outside
              tab_id[1] = v_id[tab4[i0][1]];
              tab_id[2] = p_id[tab4[i0][2]];
              tab_id[3] = p_id[tab4[i0][3]];
              tab_id[4] = v_id[tab4[i0][4]];
              tab_id[5] = p_id[tab4[i0][5]];
              this->CreateTetra(6, tab_id, cellarrayout);
            }
            else
            {
              tab_id[0] = p_id[tab4[i0][0]]; // Inside
              tab_id[1] = v_id[tab4[i0][1]];
              tab_id[2] = p_id[tab4[i0][2]];
              tab_id[3] = p_id[tab4[i0][3]];
              tab_id[4] = v_id[tab4[i0][4]];
              tab_id[5] = p_id[tab4[i0][5]];
              this->CreateTetra(6, tab_id, newcellArray);

              tab_id[0] = p_id[tab4[i0 + 1][0]]; // Outside
              tab_id[1] = v_id[tab4[i0 + 1][1]];
              tab_id[2] = p_id[tab4[i0 + 1][2]];
              tab_id[3] = p_id[tab4[i0 + 1][3]];
              tab_id[4] = v_id[tab4[i0 + 1][4]];
              tab_id[5] = p_id[tab4[i0 + 1][5]];
              this->CreateTetra(6, tab_id, cellarrayout);
            }
            break;
          case 3: // We have one tetrahedron and one wedge
            switch (edges_inter)
            {
              case 134:
                i0 = 0;
                break;
              case 125:
                i0 = 1;
                break;
              case 236:
                i0 = 2;
                break;
              case 456:
                i0 = 3;
                break;
              default:
                vtkErrorMacro(<< "Intersection not found: Num_inter = " << num_inter
                              << " Edges_inter = " << edges_inter);
                continue;
            }
            if (((v_tetra[i0][cutInd] < value) && ((planes % 2) == 0)) ||
              ((v_tetra[i0][cutInd] > value) && ((planes % 2) == 1)))
            {
              // Isolate vertex is outside box, so/ the tetrahedron is outside
              tab_id[0] = p_id[tab3[i0][0]]; // Inside
              tab_id[1] = p_id[tab3[i0][1]];
              tab_id[2] = p_id[tab3[i0][2]];
              tab_id[3] = v_id[tab3[i0][3]];
              tab_id[4] = v_id[tab3[i0][4]];
              tab_id[5] = v_id[tab3[i0][5]];
              this->CreateTetra(6, tab_id, newcellArray);

              tab_id[0] = p_id[tab3[i0][0]]; // Outside
              tab_id[1] = p_id[tab3[i0][1]];
              tab_id[2] = p_id[tab3[i0][2]];
              tab_id[3] = v_id[i0];
              newCellId = cellarrayout->InsertNextCell(4, tab_id);
            }
            else
            {
              tab_id[0] = p_id[tab3[i0][0]];
              tab_id[1] = p_id[tab3[i0][1]];
              tab_id[2] = p_id[tab3[i0][2]]; // Inside
              tab_id[3] = v_id[i0];
              newCellId = newcellArray->InsertNextCell(4, tab_id);

              tab_id[0] = p_id[tab3[i0][0]]; // Outside
              tab_id[1] = p_id[tab3[i0][1]];
              tab_id[2] = p_id[tab3[i0][2]];
              tab_id[3] = v_id[tab3[i0][3]];
              tab_id[4] = v_id[tab3[i0][4]];
              tab_id[5] = v_id[tab3[i0][5]];
              this->CreateTetra(6, tab_id, cellarrayout);
            }
            break;
          case 2:                // We have one tetrahedron and one pyramid
            switch (edges_inter) // i1 = vertex of the tetrahedron
            {
              case 12:
                i0 = 0;
                i1 = 1;
                break;
              case 13:
                i0 = 1;
                i1 = 0;
                break;
              case 23:
                i0 = 2;
                i1 = 2;
                break;
              case 25:
                i0 = 3;
                i1 = 1;
                break;
              case 26:
                i0 = 4;
                i1 = 2;
                break;
              case 56:
                i0 = 5;
                i1 = 3;
                break;
              case 34:
                i0 = 6;
                i1 = 0;
                break;
              case 46:
                i0 = 7;
                i1 = 3;
                break;
              case 36:
                i0 = 8;
                i1 = 2;
                break;
              case 14:
                i0 = 9;
                i1 = 0;
                break;
              case 15:
                i0 = 10;
                i1 = 1;
                break;
              case 45:
                i0 = 11;
                i1 = 3;
                break;
              default:
                vtkErrorMacro(<< "Intersection not found: Num_inter = " << num_inter
                              << " Edges_inter = " << edges_inter);
                continue;
            }
            if (((v_tetra[i1][cutInd] < value) && ((planes % 2) == 0)) ||
              ((v_tetra[i1][cutInd] > value) && ((planes % 2) == 1)))
            {
              // Isolate vertex is outside box, so the tetrahedron is outside
              tab_id[0] = v_id[tab2[i0][0]]; // Inside
              tab_id[1] = p_id[tab2[i0][1]];
              tab_id[2] = p_id[tab2[i0][2]];
              tab_id[3] = v_id[tab2[i0][3]];
              tab_id[4] = v_id[tab2[i0][4]];
              this->CreateTetra(5, tab_id, newcellArray);

              tab_id[0] = v_id[i1]; // Outside
              tab_id[1] = v_id[tab2[i0][4]];
              tab_id[2] = p_id[tab2[i0][2]];
              tab_id[3] = p_id[tab2[i0][1]];
              newCellId = cellarrayout->InsertNextCell(4, tab_id);
            }
            else
            {
              tab_id[0] = v_id[i1]; // Inside
              tab_id[1] = v_id[tab2[i0][4]];
              tab_id[2] = p_id[tab2[i0][2]];
              tab_id[3] = p_id[tab2[i0][1]];
              newCellId = newcellArray->InsertNextCell(4, tab_id);

              tab_id[0] = v_id[tab2[i0][0]]; // Outside
              tab_id[1] = p_id[tab2[i0][1]];
              tab_id[2] = p_id[tab2[i0][2]];
              tab_id[3] = v_id[tab2[i0][3]];
              tab_id[4] = v_id[tab2[i0][4]];
              this->CreateTetra(5, tab_id, cellarrayout);
            }
            break;
          case 1: // We have two tetrahedron.
            if ((edges_inter > 6) || (edges_inter < 1))
            {
              vtkErrorMacro(<< "Intersection not found: Num_inter = " << num_inter
                            << " Edges_inter = " << edges_inter);
              continue;
            }
            if (((v_tetra[tab1[2 * edges_inter - 1][2]][cutInd] < value) && ((planes % 2) == 0)) ||
              ((v_tetra[tab1[2 * edges_inter - 1][2]][cutInd] > value) && ((planes % 2) == 1)))
            {
              // Isolate vertex is outside box, so the tetrahedron is outside

              tab_id[0] = p_id[0]; // Inside
              tab_id[1] = v_id[tab1[2 * edges_inter - 2][0]];
              tab_id[2] = v_id[tab1[2 * edges_inter - 2][1]];
              tab_id[3] = v_id[tab1[2 * edges_inter - 2][2]];
              newCellId = newcellArray->InsertNextCell(4, tab_id);

              tab_id[0] = p_id[0]; // Outside
              tab_id[1] = v_id[tab1[2 * edges_inter - 1][0]];
              tab_id[2] = v_id[tab1[2 * edges_inter - 1][1]];
              tab_id[3] = v_id[tab1[2 * edges_inter - 1][2]];
              newCellId = cellarrayout->InsertNextCell(4, tab_id);
            }
            else
            {
              tab_id[0] = p_id[0]; // Inside
              tab_id[1] = v_id[tab1[2 * edges_inter - 1][0]];
              tab_id[2] = v_id[tab1[2 * edges_inter - 1][1]];
              tab_id[3] = v_id[tab1[2 * edges_inter - 1][2]];
              newCellId = newcellArray->InsertNextCell(4, tab_id);

              tab_id[0] = p_id[0]; // Outside
              tab_id[1] = v_id[tab1[2 * edges_inter - 2][0]];
              tab_id[2] = v_id[tab1[2 * edges_inter - 2][1]];
              tab_id[3] = v_id[tab1[2 * edges_inter - 2][2]];
              newCellId = cellarrayout->InsertNextCell(4, tab_id);
            }
            break;
        }
      } // for all new cells
      cellarray->Delete();
      cellarray = newcellArray;
    } // for all planes

    unsigned int totalnewcells = cellarray->GetNumberOfCells(); // Inside
    for (idcellnew = 0; idcellnew < totalnewcells; idcellnew++)
    {
      cellarray->GetNextCell(npts, v_id);
      newCellId = tets[0]->InsertNextCell(npts, v_id);
      outCD[0]->CopyData(inCD, cellId, newCellId);
    }
    cellarray->Delete();

    totalnewcells = cellarrayout->GetNumberOfCells(); // Outside

    for (idcellnew = 0; idcellnew < totalnewcells; idcellnew++)
    {
      cellarrayout->GetNextCell(npts, v_id);
      newCellId = tets[1]->InsertNextCell(npts, v_id);
      outCD[1]->CopyData(inCD, cellId, newCellId);
    }
    cellarrayout->Delete();
  }
  arraytetra->Delete();
}

//------------------------------------------------------------------------------
// ClipHexahedronInOut
//
// The difference between ClipHexahedron and ClipHexahedronInOut is the outputs.
// The ClipHexahedronInOut generate both outputs: inside and outside the clip hexahedron.
//
void vtkBoxClipDataSet::ClipHexahedronInOut(vtkPoints* newPoints, vtkGenericCell* cell,
  vtkIncrementalPointLocator* locator, vtkCellArray** tets, vtkPointData* inPD,
  vtkPointData** outPD, vtkCellData* inCD, vtkIdType cellId, vtkCellData** outCD)
{
  vtkIdType cellType = cell->GetCellType();
  vtkIdList* cellIds = cell->GetPointIds();
  vtkCellArray* arraytetra = vtkCellArray::New();
  vtkPoints* cellPts = cell->GetPoints();
  vtkIdType npts = cellPts->GetNumberOfPoints();
  std::vector<vtkIdType> cellptId(npts);
  vtkIdType iid[4];
  const vtkIdType* v_id = nullptr;
  vtkIdType *verts, v1, v2;
  vtkIdType ptId;
  vtkIdType ptIdout[4];
  vtkIdType tab_id[6];
  vtkIdType ptstetra = 4;

  int i, j, k;
  int allInside;

  /* Edges Tetrahedron */
  vtkIdType edges[6][2] = {
    { 0, 1 },
    { 1, 2 },
    { 2, 0 },
    { 0, 3 },
    { 1, 3 },
    { 2, 3 },
  };
  double deltaScalar;
  double p[6], t;
  double v_tetra[4][3];
  double v[3], x[3];
  double *p1, *p2;

  unsigned int idtetranew;
  unsigned int idcellnew;
  unsigned int planes;

  for (i = 0; i < npts; i++)
  {
    cellptId[i] = cellIds->GetId(i);
  }

  this->CellGrid(
    cellType, npts, cellptId.data(), arraytetra); // Convert all volume cells to tetrahedra

  unsigned int totalnewtetra = arraytetra->GetNumberOfCells();
  for (idtetranew = 0; idtetranew < totalnewtetra; idtetranew++)
  {
    arraytetra->GetNextCell(ptstetra, v_id);

    for (allInside = 1, i = 0; i < 4; i++)
    {
      cellPts->GetPoint(v_id[i], v);

      for (k = 0; k < 6; k++)
      {
        p[k] = this->PlaneNormal[k][0] * (v[0] - this->PlanePoint[k][0]) +
          this->PlaneNormal[k][1] * (v[1] - this->PlanePoint[k][1]) +
          this->PlaneNormal[k][2] * (v[2] - this->PlanePoint[k][2]);
      }

      if (!((p[0] <= 0) && (p[1] <= 0) && (p[2] <= 0) && (p[3] <= 0) && (p[4] <= 0) && (p[5] <= 0)))
      {
        allInside = 0;
      }
    } // for all points of the tetrahedron.

    // Test Outside
    unsigned int test[6] = { 1, 1, 1, 1, 1, 1 };
    for (i = 0; i < 4; i++)
    {
      ptIdout[i] = cellIds->GetId(v_id[i]);
      cellPts->GetPoint(v_id[i], v_tetra[i]);

      // Use plane equation
      for (k = 0; k < 6; k++)
      {
        p[k] = this->PlaneNormal[k][0] * (v_tetra[i][0] - this->PlanePoint[k][0]) +
          this->PlaneNormal[k][1] * (v_tetra[i][1] - this->PlanePoint[k][1]) +
          this->PlaneNormal[k][2] * (v_tetra[i][2] - this->PlanePoint[k][2]);
      }

      for (k = 0; k < 3; k++)
      {
        if (p[2 * k] < 0)
        {
          test[2 * k] = 0;
        }
        if (p[2 * k + 1] < 0)
        {
          test[2 * k + 1] = 0;
        }
      }

    } // for all points of the cell.

    if (!allInside &&
      ((test[0] == 1) || (test[1] == 1) || (test[2] == 1) || (test[3] == 1) || (test[4] == 1) ||
        (test[5] == 1)))
    {
      for (i = 0; i < 4; i++)
      {
        if (locator->InsertUniquePoint(v_tetra[i], iid[i]))
        {
          outPD[0]->CopyData(inPD, ptIdout[i], iid[i]);
          outPD[1]->CopyData(inPD, ptIdout[i], iid[i]);
        }
      }
      int newCellId = tets[1]->InsertNextCell(4, iid);
      outCD[1]->CopyData(inCD, cellId, newCellId);
      continue; // Tetrahedron is outside.
    }

    for (i = 0; i < 4; i++)
    {
      ptId = cellIds->GetId(v_id[i]);
      cellPts->GetPoint(v_id[i], v);

      // Currently all points are injected because of the possibility
      // of intersection point merging.

      if (locator->InsertUniquePoint(v, iid[i]))
      {
        outPD[0]->CopyData(inPD, ptId, iid[i]);
        outPD[1]->CopyData(inPD, ptId, iid[i]);
      }

    } // for all points of the tetrahedron.

    if (allInside)
    {
      int newCellId = tets[0]->InsertNextCell(4, iid); // Tetrahedron inside.
      outCD[0]->CopyData(inCD, cellId, newCellId);
      continue;
    }

    // float *pc1  , *pc2;
    double *pedg1, *pedg2;

    // Tetrahedron Intersection Cases
    constexpr unsigned int tab4[6][6] = {
      { 1, 1, 0, 3, 3, 2 },
      { 2, 0, 0, 3, 2, 1 },
      { 3, 3, 2, 0, 2, 1 },
      { 1, 0, 2, 0, 1, 3 },
      { 0, 0, 1, 2, 3, 3 },
      { 0, 1, 2, 1, 2, 3 },
    };
    constexpr unsigned int tab3[4][6] = {
      { 0, 2, 1, 1, 3, 2 },
      { 0, 1, 2, 0, 2, 3 },
      { 0, 1, 2, 1, 0, 3 },
      { 0, 1, 2, 0, 1, 2 },
    };
    constexpr unsigned int tab2[12][5] = {
      { 0, 0, 1, 2, 3 },
      { 2, 1, 0, 1, 3 },
      { 1, 0, 1, 0, 3 },
      { 2, 0, 1, 3, 0 },
      { 3, 1, 0, 1, 0 },
      { 1, 0, 1, 2, 0 },
      { 3, 1, 0, 2, 1 },
      { 2, 1, 0, 0, 1 },
      { 0, 0, 1, 3, 1 },
      { 1, 0, 1, 3, 2 },
      { 3, 1, 0, 0, 2 },
      { 0, 0, 1, 1, 2 },
    };
    constexpr unsigned int tab1[12][3] = {
      { 2, 3, 1 },
      { 3, 2, 0 },
      { 3, 0, 1 },
      { 0, 3, 2 },
      { 1, 3, 0 },
      { 3, 1, 2 },
      { 2, 1, 0 },
      { 1, 2, 3 },
      { 2, 0, 3 },
      { 0, 2, 1 },
      { 0, 1, 3 },
      { 1, 0, 2 },
    };

    vtkCellArray* cellarray = vtkCellArray::New();
    vtkCellArray* cellarrayout = vtkCellArray::New();
    int newCellId = cellarray->InsertNextCell(4, iid);

    // Test Cell intersection with each plane of box
    // FIXME: there is no difference in the for loop (planes has no influence!)
    for (planes = 0; planes < 6; planes++)
    {
      unsigned int totalnewcells = cellarray->GetNumberOfCells();
      vtkCellArray* newcellArray = vtkCellArray::New();

      for (idcellnew = 0; idcellnew < totalnewcells; idcellnew++)
      {
        unsigned int i0, i1;
        unsigned int num_inter = 0;
        unsigned int edges_inter = 0;
        vtkIdType p_id[4];

        cellarray->GetNextCell(npts, v_id);

        newPoints->GetPoint(v_id[0], v_tetra[0]); // coord (x,y,z)
        newPoints->GetPoint(v_id[1], v_tetra[1]);
        newPoints->GetPoint(v_id[2], v_tetra[2]);
        newPoints->GetPoint(v_id[3], v_tetra[3]);

        p[0] = this->PlaneNormal[planes][0] * (v_tetra[0][0] - this->PlanePoint[planes][0]) +
          this->PlaneNormal[planes][1] * (v_tetra[0][1] - this->PlanePoint[planes][1]) +
          this->PlaneNormal[planes][2] * (v_tetra[0][2] - this->PlanePoint[planes][2]);
        p[1] = this->PlaneNormal[planes][0] * (v_tetra[1][0] - this->PlanePoint[planes][0]) +
          this->PlaneNormal[planes][1] * (v_tetra[1][1] - this->PlanePoint[planes][1]) +
          this->PlaneNormal[planes][2] * (v_tetra[1][2] - this->PlanePoint[planes][2]);
        p[2] = this->PlaneNormal[planes][0] * (v_tetra[2][0] - this->PlanePoint[planes][0]) +
          this->PlaneNormal[planes][1] * (v_tetra[2][1] - this->PlanePoint[planes][1]) +
          this->PlaneNormal[planes][2] * (v_tetra[2][2] - this->PlanePoint[planes][2]);
        p[3] = this->PlaneNormal[planes][0] * (v_tetra[3][0] - this->PlanePoint[planes][0]) +
          this->PlaneNormal[planes][1] * (v_tetra[3][1] - this->PlanePoint[planes][1]) +
          this->PlaneNormal[planes][2] * (v_tetra[3][2] - this->PlanePoint[planes][2]);

        for (int edgeNum = 0; edgeNum < 6; edgeNum++)
        {
          verts = edges[edgeNum];

          p1 = v_tetra[verts[0]];
          p2 = v_tetra[verts[1]];
          double s1 = p[verts[0]];
          double s2 = p[verts[1]];
          if ((s1 * s2) < 0)
          {
            deltaScalar = s2 - s1;

            if (deltaScalar > 0)
            {
              pedg1 = p1;
              pedg2 = p2;
              v1 = verts[0];
              v2 = verts[1];
            }
            else
            {
              pedg1 = p2;
              pedg2 = p1;
              v1 = verts[1];
              v2 = verts[0];
              deltaScalar = -deltaScalar;
              t = s1;
              s1 = s2;
              s2 = t;
            }

            // linear interpolation
            t = (deltaScalar == 0.0 ? 0.0 : (-s1) / deltaScalar);

            for (j = 0; j < 3; j++)
            {
              x[j] = pedg1[j] + t * (pedg2[j] - pedg1[j]);
            }

            // Incorporate point into output and interpolate edge data as necessary
            edges_inter = edges_inter * 10 + (edgeNum + 1);

            if (locator->InsertUniquePoint(x, p_id[num_inter]))
            {
              this->InterpolateEdge(outPD[0], p_id[num_inter], v_id[v1], v_id[v2], t);
              this->InterpolateEdge(outPD[1], p_id[num_inter], v_id[v1], v_id[v2], t);
            }

            num_inter++;
          } // if edge intersects value
        }   // for all edges

        if (num_inter == 0)
        {
          unsigned int outside = 0;
          for (i = 0; i < 4; i++)
          {
            if (p[i] > 0)
            {              // If only one vertex is outside, so the tetrahedron is outside
              outside = 1; // because there is not intersection.
              break;       // some vertex could be on plane, so you need to test all vertex
            }
          }
          if (outside == 0)
          {
            // else it is possible intersection if other plane
            newCellId = newcellArray->InsertNextCell(4, v_id);
          }
          else
          {
            newCellId = tets[1]->InsertNextCell(4, v_id);
            outCD[1]->CopyData(inCD, cellId, newCellId);
          }
          continue;
        }

        switch (num_inter)
        {
          case 4: // We have two wedges
            switch (edges_inter)
            {
              case 1246:
                i0 = 0;
                break;
              case 2345:
                i0 = 2;
                break;
              case 1356:
                i0 = 4;
                break;
              default:
                vtkErrorMacro(<< "Intersection not found: Num_inter = " << num_inter
                              << " Edges_inter = " << edges_inter);
                continue;
            }
            if (p[3] > 0)
            {
              // The v_tetra[3] is outside box, so the first wedge is outside
              // ps: v_tetra[3] is always in first wedge (see tab)

              tab_id[0] = p_id[tab4[i0 + 1][0]]; // Inside
              tab_id[1] = v_id[tab4[i0 + 1][1]];
              tab_id[2] = p_id[tab4[i0 + 1][2]];
              tab_id[3] = p_id[tab4[i0 + 1][3]];
              tab_id[4] = v_id[tab4[i0 + 1][4]];
              tab_id[5] = p_id[tab4[i0 + 1][5]];
              this->CreateTetra(6, tab_id, newcellArray);

              tab_id[0] = p_id[tab4[i0][0]]; // Outside
              tab_id[1] = v_id[tab4[i0][1]];
              tab_id[2] = p_id[tab4[i0][2]];
              tab_id[3] = p_id[tab4[i0][3]];
              tab_id[4] = v_id[tab4[i0][4]];
              tab_id[5] = p_id[tab4[i0][5]];
              this->CreateTetra(6, tab_id, cellarrayout);
            }
            else
            {
              tab_id[0] = p_id[tab4[i0][0]]; // Inside
              tab_id[1] = v_id[tab4[i0][1]];
              tab_id[2] = p_id[tab4[i0][2]];
              tab_id[3] = p_id[tab4[i0][3]];
              tab_id[4] = v_id[tab4[i0][4]];
              tab_id[5] = p_id[tab4[i0][5]];
              this->CreateTetra(6, tab_id, newcellArray);

              tab_id[0] = p_id[tab4[i0 + 1][0]]; // Outside
              tab_id[1] = v_id[tab4[i0 + 1][1]];
              tab_id[2] = p_id[tab4[i0 + 1][2]];
              tab_id[3] = p_id[tab4[i0 + 1][3]];
              tab_id[4] = v_id[tab4[i0 + 1][4]];
              tab_id[5] = p_id[tab4[i0 + 1][5]];
              this->CreateTetra(6, tab_id, cellarrayout);
            }

            break;
          case 3: // We have one tetrahedron and one wedge
            switch (edges_inter)
            {
              case 134:
                i0 = 0;
                break;
              case 125:
                i0 = 1;
                break;
              case 236:
                i0 = 2;
                break;
              case 456:
                i0 = 3;
                break;
              default:
                vtkErrorMacro(<< "Intersection not found: Num_inter = " << num_inter
                              << " Edges_inter = " << edges_inter);
                continue;
            }
            if (p[i0] > 0)
            {
              // Isolate vertex is outside box, so the tetrahedron is outside
              tab_id[0] = p_id[tab3[i0][0]]; // Inside
              tab_id[1] = p_id[tab3[i0][1]];
              tab_id[2] = p_id[tab3[i0][2]];
              tab_id[3] = v_id[tab3[i0][3]];
              tab_id[4] = v_id[tab3[i0][4]];
              tab_id[5] = v_id[tab3[i0][5]];
              this->CreateTetra(6, tab_id, newcellArray);

              tab_id[0] = p_id[tab3[i0][0]]; // Outside
              tab_id[1] = p_id[tab3[i0][1]];
              tab_id[2] = p_id[tab3[i0][2]];
              tab_id[3] = v_id[i0];
              newCellId = cellarrayout->InsertNextCell(4, tab_id);
            }
            else
            {
              tab_id[0] = p_id[tab3[i0][0]]; // Inside
              tab_id[1] = p_id[tab3[i0][1]];
              tab_id[2] = p_id[tab3[i0][2]];
              tab_id[3] = v_id[i0];
              newCellId = newcellArray->InsertNextCell(4, tab_id);

              tab_id[0] = p_id[tab3[i0][0]]; // Outside
              tab_id[1] = p_id[tab3[i0][1]];
              tab_id[2] = p_id[tab3[i0][2]];
              tab_id[3] = v_id[tab3[i0][3]];
              tab_id[4] = v_id[tab3[i0][4]];
              tab_id[5] = v_id[tab3[i0][5]];
              this->CreateTetra(6, tab_id, cellarrayout);
            }
            break;
          case 2:                // We have one tetrahedron and one pyramid
            switch (edges_inter) // i1 = vertex of the tetrahedron
            {
              case 12:
                i0 = 0;
                i1 = 1;
                break;
              case 13:
                i0 = 1;
                i1 = 0;
                break;
              case 23:
                i0 = 2;
                i1 = 2;
                break;
              case 25:
                i0 = 3;
                i1 = 1;
                break;
              case 26:
                i0 = 4;
                i1 = 2;
                break;
              case 56:
                i0 = 5;
                i1 = 3;
                break;
              case 34:
                i0 = 6;
                i1 = 0;
                break;
              case 46:
                i0 = 7;
                i1 = 3;
                break;
              case 36:
                i0 = 8;
                i1 = 2;
                break;
              case 14:
                i0 = 9;
                i1 = 0;
                break;
              case 15:
                i0 = 10;
                i1 = 1;
                break;
              case 45:
                i0 = 11;
                i1 = 3;
                break;
              default:
                vtkErrorMacro(<< "Intersection not found: Num_inter = " << num_inter
                              << " Edges_inter = %" << edges_inter);
                continue;
            }

            if (p[i1] > 0)
            {
              // Isolate vertex is outside box, so the tetrahedron is outside
              tab_id[0] = v_id[tab2[i0][0]]; // Inside
              tab_id[1] = p_id[tab2[i0][1]];
              tab_id[2] = p_id[tab2[i0][2]];
              tab_id[3] = v_id[tab2[i0][3]];
              tab_id[4] = v_id[tab2[i0][4]];
              this->CreateTetra(5, tab_id, newcellArray);

              tab_id[0] = v_id[i1]; // Outside
              tab_id[1] = v_id[tab2[i0][4]];
              tab_id[2] = p_id[tab2[i0][2]];
              tab_id[3] = p_id[tab2[i0][1]];
              newCellId = cellarrayout->InsertNextCell(4, tab_id);
            }
            else
            {
              tab_id[0] = v_id[i1]; // Inside
              tab_id[1] = v_id[tab2[i0][4]];
              tab_id[2] = p_id[tab2[i0][2]];
              tab_id[3] = p_id[tab2[i0][1]];
              newCellId = newcellArray->InsertNextCell(4, tab_id);

              tab_id[0] = v_id[tab2[i0][0]]; // Outside
              tab_id[1] = p_id[tab2[i0][1]];
              tab_id[2] = p_id[tab2[i0][2]];
              tab_id[3] = v_id[tab2[i0][3]];
              tab_id[4] = v_id[tab2[i0][4]];
              this->CreateTetra(5, tab_id, cellarrayout);
            }
            break;
          case 1: // We have two tetrahedron.
            if ((edges_inter > 6) || (edges_inter < 1))
            {
              vtkErrorMacro(<< "Intersection not found: Num_inter = " << num_inter
                            << " Edges_inter = " << edges_inter);
              continue;
            }
            if (p[tab1[2 * edges_inter - 1][2]] > 0)
            {
              // Isolate vertex is outside box, so the tetrahedron is outside
              tab_id[0] = p_id[0]; // Inside
              tab_id[1] = v_id[tab1[2 * edges_inter - 2][0]];
              tab_id[2] = v_id[tab1[2 * edges_inter - 2][1]];
              tab_id[3] = v_id[tab1[2 * edges_inter - 2][2]];
              newCellId = newcellArray->InsertNextCell(4, tab_id);

              tab_id[0] = p_id[0]; // Outside
              tab_id[1] = v_id[tab1[2 * edges_inter - 1][0]];
              tab_id[2] = v_id[tab1[2 * edges_inter - 1][1]];
              tab_id[3] = v_id[tab1[2 * edges_inter - 1][2]];
              newCellId = cellarrayout->InsertNextCell(4, tab_id);
            }
            else
            {
              tab_id[0] = p_id[0]; // Inside
              tab_id[1] = v_id[tab1[2 * edges_inter - 1][0]];
              tab_id[2] = v_id[tab1[2 * edges_inter - 1][1]];
              tab_id[3] = v_id[tab1[2 * edges_inter - 1][2]];
              newCellId = newcellArray->InsertNextCell(4, tab_id);

              tab_id[0] = p_id[0]; // Outside
              tab_id[1] = v_id[tab1[2 * edges_inter - 2][0]];
              tab_id[2] = v_id[tab1[2 * edges_inter - 2][1]];
              tab_id[3] = v_id[tab1[2 * edges_inter - 2][2]];
              newCellId = cellarrayout->InsertNextCell(4, tab_id);
            }
            break;
        }
      } // for all new cells
      cellarray->Delete();
      cellarray = newcellArray;
    } // for all planes

    unsigned int totalnewcells = cellarray->GetNumberOfCells(); // Inside

    for (idcellnew = 0; idcellnew < totalnewcells; idcellnew++)
    {
      cellarray->GetNextCell(npts, v_id);
      newCellId = tets[0]->InsertNextCell(npts, v_id);
      outCD[0]->CopyData(inCD, cellId, newCellId);
    }
    cellarray->Delete();

    totalnewcells = cellarrayout->GetNumberOfCells(); // Outside

    for (idcellnew = 0; idcellnew < totalnewcells; idcellnew++)
    {
      cellarrayout->GetNextCell(npts, v_id);
      newCellId = tets[1]->InsertNextCell(npts, v_id);
      outCD[1]->CopyData(inCD, cellId, newCellId);
    }
    cellarrayout->Delete();
  }
  arraytetra->Delete();
}

//-------------------------------------------------------

//-------------------------------------------------------
void vtkBoxClipDataSet::ClipBox2D(vtkPoints* newPoints, vtkGenericCell* cell,
  vtkIncrementalPointLocator* locator, vtkCellArray* tets, vtkPointData* inPD, vtkPointData* outPD,
  vtkCellData* inCD, vtkIdType cellId, vtkCellData* outCD)
{
  vtkIdType cellType = cell->GetCellType();
  vtkIdList* cellIds = cell->GetPointIds();
  vtkCellArray* arraytriangle = vtkCellArray::New();
  vtkPoints* cellPts = cell->GetPoints();
  vtkIdType npts = cellPts->GetNumberOfPoints();
  std::vector<vtkIdType> cellptId(npts);
  vtkIdType iid[3];
  const vtkIdType* v_id = nullptr;
  vtkIdType *verts, v1, v2;
  vtkIdType ptId;
  vtkIdType tab_id[3];
  vtkIdType ptstriangle = 3;

  int i, j;
  unsigned int allInside;
  unsigned int planes;
  unsigned int cutInd;

  /* Edges Triangle*/
  vtkIdType edges[3][2] = {
    { 0, 1 },
    { 1, 2 },
    { 2, 0 },
  };
  double value, deltaScalar;
  double t, *p1, *p2;
  double v[3], x[3];
  double v_triangle[3][3];

  // To avoid false positive warning.
  tab_id[0] = 0;
  tab_id[1] = 0;
  tab_id[2] = 0;

  for (i = 0; i < npts; i++)
  {
    cellptId[i] = cellIds->GetId(i);
  }

  // Convert all 2d cells to triangle
  this->CellGrid(cellType, npts, cellptId.data(), arraytriangle);

  unsigned int totalnewtriangle = arraytriangle->GetNumberOfCells();
  unsigned int idtrianglenew;

  for (idtrianglenew = 0; idtrianglenew < totalnewtriangle; idtrianglenew++)
  {
    arraytriangle->GetNextCell(ptstriangle, v_id);

    for (allInside = 1, i = 0; i < 3; i++)
    {
      cellPts->GetPoint(v_id[i], v);

      if (!(((v[0] >= this->BoundBoxClip[0][0]) && (v[0] <= this->BoundBoxClip[0][1]) &&
            (v[1] >= this->BoundBoxClip[1][0]) && (v[1] <= this->BoundBoxClip[1][1]) &&
            (v[2] >= this->BoundBoxClip[2][0]) && (v[2] <= this->BoundBoxClip[2][1]))))
      {
        // outside,its type might change later (nearby intersection)
        allInside = 0;
      }
    }

    // Test Outside:
    if (!allInside)
    {
      unsigned int test[6] = { 1, 1, 1, 1, 1, 1 };
      for (i = 0; i < 3; i++)
      {
        cellPts->GetPoint(v_id[i], v);

        if (v[0] >= this->BoundBoxClip[0][0])
        {
          test[0] = 0;
        }
        if (v[0] <= this->BoundBoxClip[0][1])
        {
          test[1] = 0;
        }
        if (v[1] >= this->BoundBoxClip[1][0])
        {
          test[2] = 0;
        }
        if (v[1] <= this->BoundBoxClip[1][1])
        {
          test[3] = 0;
        }
        if (v[2] >= this->BoundBoxClip[2][0])
        {
          test[4] = 0;
        }
        if (v[2] <= this->BoundBoxClip[2][1])
        {
          test[5] = 0;
        }
      } // for all points of the triangle.

      if ((test[0] == 1) || (test[1] == 1) || (test[2] == 1) || (test[3] == 1) || (test[4] == 1) ||
        (test[5] == 1))
      {
        continue; // Triangle is outside.
      }
    }

    for (i = 0; i < 3; i++)
    {
      // Currently all points are injected because of the possibility
      // of intersection point merging.
      ptId = cellIds->GetId(v_id[i]);
      cellPts->GetPoint(v_id[i], v);
      if (locator->InsertUniquePoint(v, iid[i]))
      {
        outPD->CopyData(inPD, ptId, iid[i]);
      }

    } // for all points of the triangle.

    if (allInside)
    {
      // Triangle inside.
      int newCellId = tets->InsertNextCell(3, iid);
      outCD->CopyData(inCD, cellId, newCellId);
      continue;
    }

    // float *pc1  , *pc2;
    double *pedg1, *pedg2;

    // Triangle intersection cases

    unsigned int tab2[3][4] = {
      { 1, 2, 1, 0 },
      { 2, 0, 0, 1 },
      { 0, 1, 0, 1 },
    };
    unsigned int tab1[3][2] = {
      { 2, 1 },
      { 0, 2 },
      { 1, 0 },
    };

    vtkCellArray* cellarray = vtkCellArray::New();
    int newCellId = cellarray->InsertNextCell(3, iid);
    unsigned int idcellnew;

    // Test triangle intersection with each plane of box
    for (planes = 0; planes < 6; planes++)
    {
      // The index of the dimension of the cut plane (x == 0, y == 1, z == 2).
      cutInd = planes / 2;

      // The plane is always parallel to unitary cube.
      value = this->BoundBoxClip[cutInd][planes % 2];

      unsigned int totalnewcells = cellarray->GetNumberOfCells();
      vtkCellArray* newcellArray = vtkCellArray::New();

      for (idcellnew = 0; idcellnew < totalnewcells; idcellnew++)
      {
        unsigned int num_inter = 0;
        unsigned int edges_inter = 0;
        unsigned int i0;
        vtkIdType p_id[3];
        cellarray->GetNextCell(npts, v_id);

        newPoints->GetPoint(v_id[0], v_triangle[0]); // coord (x,y,z)
        newPoints->GetPoint(v_id[1], v_triangle[1]);
        newPoints->GetPoint(v_id[2], v_triangle[2]);
        for (int edgeNum = 0; edgeNum < 3; edgeNum++)
        {
          verts = edges[edgeNum];

          p1 = v_triangle[verts[0]];
          p2 = v_triangle[verts[1]];

          if ((p1[cutInd] < value && value < p2[cutInd]) ||
            (p2[cutInd] < value && value < p1[cutInd]))
          {
            deltaScalar = p2[cutInd] - p1[cutInd];

            if (deltaScalar > 0)
            {
              pedg1 = p1;
              pedg2 = p2;
              v1 = verts[0];
              v2 = verts[1];
            }
            else
            {
              pedg1 = p2;
              pedg2 = p1;
              v1 = verts[1];
              v2 = verts[0];
              deltaScalar = -deltaScalar;
            }

            // linear interpolation
            t = (deltaScalar == 0.0 ? 0.0 : (value - pedg1[cutInd]) / deltaScalar);

            for (j = 0; j < 3; j++)
            {
              x[j] = pedg1[j] + t * (pedg2[j] - pedg1[j]);
            }

            // Incorporate point into output and interpolate edge data as necessary
            edges_inter = edges_inter * 10 + (edgeNum + 1);
            if (locator->InsertUniquePoint(x, p_id[num_inter]))
            {
              this->InterpolateEdge(outPD, p_id[num_inter], v_id[v1], v_id[v2], t);
            }

            num_inter++;
          } // if edge intersects value
        }   // for all edges

        if (num_inter == 0)
        {
          unsigned int outside = 0;
          for (i = 0; i < 3; i++)
          {
            if (((v_triangle[i][cutInd] < value) && ((planes % 2) == 0)) ||
              // If only one vertex is outside, so the triangle is outside
              ((v_triangle[i][cutInd] > value) &&
                ((planes % 2) == 1))) // because there is not intersection.
            {
              outside = 1;
              break;
            }
          }
          if (outside == 0)
          {
            // else it is possible intersection if other plane
            newCellId = newcellArray->InsertNextCell(3, v_id);
          }
          continue;
        }
        switch (num_inter)
        {
          case 2: // We have one quad and one triangle
            switch (edges_inter)
            {
              case 12:
                i0 = 1;
                break;
              case 23:
                i0 = 2;
                break;
              case 13:
                i0 = 0;
                break;
              default:
                vtkErrorMacro(<< "Intersection not found: Num_inter = " << num_inter
                              << " Edges_inter = " << edges_inter);
                continue;
            }
            if (((v_triangle[i0][cutInd] < value) && ((planes % 2) == 0)) ||
              ((v_triangle[i0][cutInd] > value) && ((planes % 2) == 1)))
            {
              // The v_triangle[i0] is outside box, so
              // The Quad is inside: two triangles: (v0,v1,p0) and (p0,p1,v1)
              tab_id[0] = v_id[tab2[i0][0]];
              tab_id[1] = v_id[tab2[i0][1]];
              tab_id[2] = p_id[tab2[i0][2]];
              newCellId = newcellArray->InsertNextCell(3, tab_id);
              tab_id[0] = p_id[tab2[i0][2]];
              tab_id[1] = p_id[tab2[i0][3]];
              tab_id[2] = v_id[tab2[i0][0]];
              newCellId = newcellArray->InsertNextCell(3, tab_id);
            }
            else
            {
              // The Triangle is inside: (v0,p0,p1)
              // The correct winding of the new triangle depends on where the
              // plane intersected the original triangle.
              switch (edges_inter)
              {
                case 12:
                case 23:
                  tab_id[0] = v_id[i0];
                  tab_id[1] = p_id[1];
                  tab_id[2] = p_id[0];
                  break;
                case 13:
                  tab_id[0] = v_id[i0];
                  tab_id[1] = p_id[0];
                  tab_id[2] = p_id[1];
                  break;
              }
              newCellId = newcellArray->InsertNextCell(3, tab_id);
            }
            break;

          case 1: // We have two triangles
            switch (edges_inter)
            {
              case 1:
                i0 = 0;
                break;
              case 2:
                i0 = 1;
                break;
              case 3:
                i0 = 2;
                break;
              default:
                vtkErrorMacro(<< "Intersection not found: Num_inter = " << num_inter
                              << " Edges_inter = " << edges_inter);
                continue;
            }
            if (((v_triangle[i0][cutInd] < value) && ((planes % 2) == 0)) ||
              ((v_triangle[i0][cutInd] > value) && ((planes % 2) == 1)))
            {
              // Test one of the vertices  vertex i0 is outside
              tab_id[0] = v_id[tab1[i0][1]];
              tab_id[1] = v_id[tab1[i0][0]];
              tab_id[2] = p_id[0];
              newCellId = newcellArray->InsertNextCell(3, tab_id);
            }
            else
            {
              tab_id[0] = v_id[tab1[i0][0]];
              tab_id[1] = v_id[i0];
              tab_id[2] = p_id[0];
              newCellId = newcellArray->InsertNextCell(3, tab_id);
            }
            break;
        }
      } // for all new cells
      cellarray->Delete();
      cellarray = newcellArray;
    } // for all planes

    unsigned int totalnewcells = cellarray->GetNumberOfCells();

    for (idcellnew = 0; idcellnew < totalnewcells; idcellnew++)
    {
      cellarray->GetNextCell(npts, v_id);
      newCellId = tets->InsertNextCell(npts, v_id);
      outCD->CopyData(inCD, cellId, newCellId);
    }
    cellarray->Delete();
  }
  arraytriangle->Delete();
}
//-------------------------------------------------------

void vtkBoxClipDataSet::ClipBoxInOut2D(vtkPoints* newPoints, vtkGenericCell* cell,
  vtkIncrementalPointLocator* locator, vtkCellArray** tets, vtkPointData* inPD,
  vtkPointData** outPD, vtkCellData* inCD, vtkIdType cellId, vtkCellData** outCD)
{
  vtkIdType cellType = cell->GetCellType();
  vtkIdList* cellIds = cell->GetPointIds();
  vtkCellArray* arraytriangle = vtkCellArray::New();
  vtkPoints* cellPts = cell->GetPoints();
  vtkIdType npts = cellPts->GetNumberOfPoints();
  std::vector<vtkIdType> cellptId(npts);
  vtkIdType iid[3];
  const vtkIdType* v_id = nullptr;
  vtkIdType *verts, v1, v2;
  vtkIdType ptId;
  vtkIdType ptIdout[4];
  vtkIdType tab_id[3];
  vtkIdType ptstriangle = 3;

  int i, j;
  unsigned int allInside;
  unsigned int cutInd;
  unsigned int planes;
  unsigned int idcellnew;

  /* Edges Triangle */
  vtkIdType edges[3][2] = {
    { 0, 1 },
    { 1, 2 },
    { 2, 0 },
  };
  double value, deltaScalar;
  double t, *p1, *p2;
  double v[3], x[3];
  double v_triangle[3][3];

  // To avoid false positive warning.
  tab_id[0] = 0;
  tab_id[1] = 0;
  tab_id[2] = 0;

  for (i = 0; i < npts; i++)
  {
    cellptId[i] = cellIds->GetId(i);
  }

  // Convert all 2D cells to triangle
  this->CellGrid(cellType, npts, cellptId.data(), arraytriangle);
  unsigned int totalnewtriangle = arraytriangle->GetNumberOfCells();
  unsigned int idtrianglenew;

  for (idtrianglenew = 0; idtrianglenew < totalnewtriangle; idtrianglenew++)
  {
    arraytriangle->GetNextCell(ptstriangle, v_id);

    for (allInside = 1, i = 0; i < 3; i++)
    {
      cellPts->GetPoint(v_id[i], v);

      if (!(((v[0] >= this->BoundBoxClip[0][0]) && (v[0] <= this->BoundBoxClip[0][1]) &&
            (v[1] >= this->BoundBoxClip[1][0]) && (v[1] <= this->BoundBoxClip[1][1]) &&
            (v[2] >= this->BoundBoxClip[2][0]) && (v[2] <= this->BoundBoxClip[2][1]))))
      {
        // outside,its type might change later (nearby intersection)
        allInside = 0;
      }
    } // for all points of the cell.

    // Test Outside: see(1)
    if (!allInside)
    {
      unsigned int test[6] = { 1, 1, 1, 1, 1, 1 };
      for (i = 0; i < 3; i++)
      {
        ptIdout[i] = cellIds->GetId(v_id[i]);
        cellPts->GetPoint(v_id[i], v_triangle[i]);

        if (v_triangle[i][0] >= this->BoundBoxClip[0][0])
        {
          test[0] = 0;
        }
        if (v_triangle[i][0] <= this->BoundBoxClip[0][1])
        {
          test[1] = 0;
        }
        if (v_triangle[i][1] >= this->BoundBoxClip[1][0])
        {
          test[2] = 0;
        }
        if (v_triangle[i][1] <= this->BoundBoxClip[1][1])
        {
          test[3] = 0;
        }
        if (v_triangle[i][2] >= this->BoundBoxClip[2][0])
        {
          test[4] = 0;
        }
        if (v_triangle[i][2] <= this->BoundBoxClip[2][1])
        {
          test[5] = 0;
        }

      } // for all points of the cell.

      if ((test[0] == 1) || (test[1] == 1) || (test[2] == 1) || (test[3] == 1) || (test[4] == 1) ||
        (test[5] == 1))
      {
        for (i = 0; i < 3; i++)
        {
          if (locator->InsertUniquePoint(v_triangle[i], iid[i]))
          {
            outPD[0]->CopyData(inPD, ptIdout[i], iid[i]);
            outPD[1]->CopyData(inPD, ptIdout[i], iid[i]);
          }
        }

        int newCellId = tets[1]->InsertNextCell(3, iid);
        outCD[1]->CopyData(inCD, cellId, newCellId);
        continue; // Triangle is outside.
      }
    } // if not allInside.

    for (i = 0; i < 3; i++)
    {
      ptId = cellIds->GetId(v_id[i]);
      cellPts->GetPoint(v_id[i], v);

      // Currently all points are injected because of the possibility
      // of intersection point merging.
      if (locator->InsertUniquePoint(v, iid[i]))
      {
        outPD[0]->CopyData(inPD, ptId, iid[i]);
        outPD[1]->CopyData(inPD, ptId, iid[i]);
      }
    } // for all points of the triangle.

    if (allInside)
    {
      // Triangle inside.
      int newCellId = tets[0]->InsertNextCell(3, iid);
      outCD[0]->CopyData(inCD, cellId, newCellId);
      continue;
    }

    // float *pc1  , *pc2;
    double *pedg1, *pedg2;

    // Triangle intersection cases

    unsigned int tab2[3][4] = {
      { 1, 2, 1, 0 },
      { 2, 0, 0, 1 },
      { 0, 1, 0, 1 },
    };
    unsigned int tab1[3][2] = {
      { 2, 1 },
      { 0, 2 },
      { 1, 0 },
    };

    vtkCellArray* cellarray = vtkCellArray::New();
    vtkCellArray* cellarrayout = vtkCellArray::New();
    int newCellId = cellarray->InsertNextCell(3, iid);

    // Test Cell intersection with each plane of box
    for (planes = 0; planes < 6; planes++)
    {
      // The index of the dimension of the cut plane (x == 0, y == 1, z == 2).
      cutInd = planes / 2;

      // The plane is always parallel to unitary cube.
      value = this->BoundBoxClip[cutInd][planes % 2];

      unsigned int totalnewcells = cellarray->GetNumberOfCells();
      vtkCellArray* newcellArray = vtkCellArray::New();

      for (idcellnew = 0; idcellnew < totalnewcells; idcellnew++)
      {
        unsigned int num_inter = 0;
        unsigned int edges_inter = 0;
        unsigned int i0;
        vtkIdType p_id[3];
        cellarray->GetNextCell(npts, v_id);

        newPoints->GetPoint(v_id[0], v_triangle[0]); // coord (x,y,z)
        newPoints->GetPoint(v_id[1], v_triangle[1]);
        newPoints->GetPoint(v_id[2], v_triangle[2]);
        for (int edgeNum = 0; edgeNum < 3; edgeNum++)
        {
          verts = edges[edgeNum];

          p1 = v_triangle[verts[0]];
          p2 = v_triangle[verts[1]];

          if ((p1[cutInd] < value && value < p2[cutInd]) ||
            (p2[cutInd] < value && value < p1[cutInd]))
          {
            deltaScalar = p2[cutInd] - p1[cutInd];

            if (deltaScalar > 0)
            {
              pedg1 = p1;
              pedg2 = p2;
              v1 = verts[0];
              v2 = verts[1];
            }
            else
            {
              pedg1 = p2;
              pedg2 = p1;
              v1 = verts[1];
              v2 = verts[0];
              deltaScalar = -deltaScalar;
            }

            // linear interpolation
            t = (deltaScalar == 0.0 ? 0.0 : (value - pedg1[cutInd]) / deltaScalar);

            for (j = 0; j < 3; j++)
            {
              x[j] = pedg1[j] + t * (pedg2[j] - pedg1[j]);
            }

            // Incorporate point into output and interpolate edge data as necessary
            edges_inter = edges_inter * 10 + (edgeNum + 1);
            if (locator->InsertUniquePoint(x, p_id[num_inter]))
            {
              this->InterpolateEdge(outPD[0], p_id[num_inter], v_id[v1], v_id[v2], t);
              this->InterpolateEdge(outPD[1], p_id[num_inter], v_id[v1], v_id[v2], t);
            }

            num_inter++;
          } // if edge intersects value
        }   // for all edges

        if (num_inter == 0)
        {
          unsigned int outside = 0;
          for (i = 0; i < 3; i++)
          {
            if (((v_triangle[i][cutInd] < value) && ((planes % 2) == 0)) ||
              ((v_triangle[i][cutInd] > value) && ((planes % 2) == 1)))
            {
              // If only one vertex is outside, so the triangle is outside
              // because there is not intersection.
              outside = 1;
              break;
            }
          }
          if (outside == 0)
          {
            // else it is possible intersection if other plane
            newCellId = newcellArray->InsertNextCell(3, v_id);
          }
          else
          {
            newCellId = tets[1]->InsertNextCell(3, v_id);
            outCD[1]->CopyData(inCD, cellId, newCellId);
          }
          continue;
        }
        switch (num_inter)
        {
          case 2: // We have one quad and one triangle
            // i0 gets the index of the triangle point that lies alone on
            // one side of the plane.
            switch (edges_inter)
            {
              case 12:
                i0 = 1;
                break;
              case 23:
                i0 = 2;
                break;
              case 13:
                i0 = 0;
                break;
              default:
                vtkErrorMacro(<< "Intersection not found: Num_inter = " << num_inter
                              << " Edges_inter = " << edges_inter);
                continue;
            }
            if (((v_triangle[i0][cutInd] < value) && ((planes % 2) == 0)) ||
              ((v_triangle[i0][cutInd] > value) && ((planes % 2) == 1)))
            {
              // The v_triangle[i0] is outside box, so

              tab_id[0] = v_id[tab2[i0][0]]; // Quad Inside
              tab_id[1] = v_id[tab2[i0][1]];
              tab_id[2] = p_id[tab2[i0][2]];
              newCellId = newcellArray->InsertNextCell(3, tab_id);
              tab_id[0] = p_id[tab2[i0][2]];
              tab_id[1] = p_id[tab2[i0][3]];
              tab_id[2] = v_id[tab2[i0][0]];
              newCellId = newcellArray->InsertNextCell(3, tab_id);

              switch (edges_inter) // Triangle Outside
              {
                case 12:
                case 23:
                  tab_id[0] = v_id[i0];
                  tab_id[1] = p_id[1];
                  tab_id[2] = p_id[0];
                  break;
                case 13:
                  tab_id[0] = v_id[i0];
                  tab_id[1] = p_id[0];
                  tab_id[2] = p_id[1];
                  break;
              }
              newCellId = cellarrayout->InsertNextCell(3, tab_id);
            }
            else
            {
              // The Triangle is inside: (v0,p0,p1)
              switch (edges_inter)
              {
                case 12:
                case 23:
                  tab_id[0] = v_id[i0];
                  tab_id[1] = p_id[1];
                  tab_id[2] = p_id[0];
                  break;
                case 13:
                  tab_id[0] = v_id[i0];
                  tab_id[1] = p_id[0];
                  tab_id[2] = p_id[1];
                  break;
              }
              newCellId = newcellArray->InsertNextCell(3, tab_id);

              tab_id[0] = v_id[tab2[i0][0]]; // Quad is Outside
              tab_id[1] = v_id[tab2[i0][1]];
              tab_id[2] = p_id[tab2[i0][2]];
              newCellId = cellarrayout->InsertNextCell(3, tab_id);
              tab_id[0] = p_id[tab2[i0][2]];
              tab_id[1] = p_id[tab2[i0][3]];
              tab_id[2] = v_id[tab2[i0][0]];
              newCellId = cellarrayout->InsertNextCell(3, tab_id);
            }
            break;

          case 1: // We have two triangles
            switch (edges_inter)
            {
              case 1:
                i0 = 0;
                break;
              case 2:
                i0 = 1;
                break;
              case 3:
                i0 = 2;
                break;
              default:
                vtkErrorMacro(<< "Intersection not found: Num_inter = " << num_inter
                              << " Edges_inter = " << edges_inter);
                continue;
            }
            if (((v_triangle[i0][cutInd] < value) && ((planes % 2) == 0)) ||
              ((v_triangle[i0][cutInd] > value) && ((planes % 2) == 1)))
            {
              // Test one of the vertices vertex i0 is outside

              tab_id[0] = v_id[tab1[i0][1]]; // Inside
              tab_id[1] = v_id[tab1[i0][0]];
              tab_id[2] = p_id[0];
              newCellId = newcellArray->InsertNextCell(3, tab_id);

              tab_id[0] = v_id[tab1[i0][0]]; // Outside
              tab_id[1] = v_id[i0];
              tab_id[2] = p_id[0];
              newCellId = cellarrayout->InsertNextCell(3, tab_id);
            }
            else
            {
              tab_id[0] = v_id[tab1[i0][0]]; // Inside
              tab_id[1] = v_id[i0];
              tab_id[2] = p_id[0];
              newCellId = newcellArray->InsertNextCell(3, tab_id);

              tab_id[0] = v_id[tab1[i0][1]]; // Outside
              tab_id[1] = v_id[tab1[i0][0]];
              tab_id[2] = p_id[0];
              newCellId = cellarrayout->InsertNextCell(3, tab_id);
            }
            break;
        }
      } // for all new cells
      cellarray->Delete();
      cellarray = newcellArray;
    } // for all planes

    unsigned int totalnewcells = cellarray->GetNumberOfCells(); // Inside

    for (idcellnew = 0; idcellnew < totalnewcells; idcellnew++)
    {
      cellarray->GetNextCell(npts, v_id);
      newCellId = tets[0]->InsertNextCell(npts, v_id);
      outCD[0]->CopyData(inCD, cellId, newCellId);
    }
    cellarray->Delete();

    totalnewcells = cellarrayout->GetNumberOfCells(); // Outside

    for (idcellnew = 0; idcellnew < totalnewcells; idcellnew++)
    {
      cellarrayout->GetNextCell(npts, v_id);
      newCellId = tets[1]->InsertNextCell(npts, v_id);
      outCD[1]->CopyData(inCD, cellId, newCellId);
    }
    cellarrayout->Delete();
  }
  arraytriangle->Delete();
}

//-------------------------------------------------------

void vtkBoxClipDataSet::ClipHexahedron2D(vtkPoints* newPoints, vtkGenericCell* cell,
  vtkIncrementalPointLocator* locator, vtkCellArray* tets, vtkPointData* inPD, vtkPointData* outPD,
  vtkCellData* inCD, vtkIdType cellId, vtkCellData* outCD)
{
  vtkIdType cellType = cell->GetCellType();
  vtkIdList* cellIds = cell->GetPointIds();
  vtkCellArray* arraytriangle = vtkCellArray::New();
  vtkPoints* cellPts = cell->GetPoints();
  vtkIdType npts = cellPts->GetNumberOfPoints();
  std::vector<vtkIdType> cellptId(npts);
  vtkIdType iid[3];
  const vtkIdType* v_id = nullptr;
  vtkIdType *verts, v1, v2;
  vtkIdType ptId;
  vtkIdType tab_id[3];
  vtkIdType ptstriangle = 3;

  int i, j, k;
  unsigned int allInside;
  unsigned int idtrianglenew;
  unsigned int idcellnew;
  unsigned int planes;

  vtkIdType edges[3][2] = {
    { 0, 1 },
    { 1, 2 },
    { 2, 0 },
  };
  double deltaScalar;
  double t, *p1, *p2;
  double v[3], x[3];
  double p[6];
  double v_triangle[3][3];

  // To avoid false positive warning.
  tab_id[0] = 0;
  tab_id[1] = 0;
  tab_id[2] = 0;

  for (i = 0; i < npts; i++)
  {
    cellptId[i] = cellIds->GetId(i);
  }

  this->CellGrid(
    cellType, npts, cellptId.data(), arraytriangle); // Convert all volume cells to triangle

  unsigned int totalnewtriangle = arraytriangle->GetNumberOfCells();
  for (idtrianglenew = 0; idtrianglenew < totalnewtriangle; idtrianglenew++)
  {
    arraytriangle->GetNextCell(ptstriangle, v_id);

    for (allInside = 1, i = 0; i < 3; i++)
    {
      cellPts->GetPoint(v_id[i], v);

      for (k = 0; k < 6; k++)
      {
        p[k] = this->PlaneNormal[k][0] * (v[0] - this->PlanePoint[k][0]) +
          this->PlaneNormal[k][1] * (v[1] - this->PlanePoint[k][1]) +
          this->PlaneNormal[k][2] * (v[2] - this->PlanePoint[k][2]);
      }

      if (!((p[0] <= 0) && (p[1] <= 0) && (p[2] <= 0) && (p[3] <= 0) && (p[4] <= 0) && (p[5] <= 0)))
      {
        allInside = 0;
      }
    } // for all points of the triangle.

    // Test Outside
    unsigned int test[6] = { 1, 1, 1, 1, 1, 1 };
    for (i = 0; i < 3; i++)
    {
      cellPts->GetPoint(v_id[i], v);

      // Use plane equation
      for (k = 0; k < 6; k++)
      {
        p[k] = this->PlaneNormal[k][0] * (v[0] - this->PlanePoint[k][0]) +
          this->PlaneNormal[k][1] * (v[1] - this->PlanePoint[k][1]) +
          this->PlaneNormal[k][2] * (v[2] - this->PlanePoint[k][2]);
      }

      for (k = 0; k < 3; k++)
      {
        if (p[2 * k] <= 0)
        {
          test[2 * k] = 0;
        }
        if (p[2 * k + 1] <= 0)
        {
          test[2 * k + 1] = 0;
        }
      }
    } // for all points of the cell.

    if (!allInside &&
      ((test[0] == 1) || (test[1] == 1) || (test[2] == 1) || (test[3] == 1) || (test[4] == 1) ||
        (test[5] == 1)))
    {
      continue; // Triangle is outside.
    }

    for (i = 0; i < 3; i++)
    {
      ptId = cellIds->GetId(v_id[i]);
      cellPts->GetPoint(v_id[i], v);

      // Currently all points are injected because of the possibility
      // of intersection point merging.

      if (locator->InsertUniquePoint(v, iid[i]))
      {
        outPD->CopyData(inPD, ptId, iid[i]);
      }
    } // for all points of the triangle.

    if (allInside)
    {
      int newCellId = tets->InsertNextCell(3, iid); // Triangle inside.
      outCD->CopyData(inCD, cellId, newCellId);
      continue;
    }

    // float *pc1  , *pc2;
    double *pedg1, *pedg2;

    unsigned int tab2[3][4] = {
      { 1, 2, 1, 0 },
      { 2, 0, 0, 1 },
      { 0, 1, 0, 1 },
    };
    unsigned int tab1[3][2] = {
      { 2, 1 },
      { 0, 2 },
      { 1, 0 },
    };

    vtkCellArray* cellarray = vtkCellArray::New();
    int newCellId = cellarray->InsertNextCell(3, iid);

    // Test Cell intersection with each plane of box
    for (planes = 0; planes < 6; planes++)
    {
      unsigned int totalnewcells = cellarray->GetNumberOfCells();
      vtkCellArray* newcellArray = vtkCellArray::New();

      for (idcellnew = 0; idcellnew < totalnewcells; idcellnew++)
      {
        unsigned int i0;
        unsigned int num_inter = 0;
        unsigned int edges_inter = 0;
        vtkIdType p_id[3];

        cellarray->GetNextCell(npts, v_id);

        newPoints->GetPoint(v_id[0], v_triangle[0]); // coord (x,y,z)
        newPoints->GetPoint(v_id[1], v_triangle[1]);
        newPoints->GetPoint(v_id[2], v_triangle[2]);

        p[0] = this->PlaneNormal[planes][0] * (v_triangle[0][0] - this->PlanePoint[planes][0]) +
          this->PlaneNormal[planes][1] * (v_triangle[0][1] - this->PlanePoint[planes][1]) +
          this->PlaneNormal[planes][2] * (v_triangle[0][2] - this->PlanePoint[planes][2]);
        p[1] = this->PlaneNormal[planes][0] * (v_triangle[1][0] - this->PlanePoint[planes][0]) +
          this->PlaneNormal[planes][1] * (v_triangle[1][1] - this->PlanePoint[planes][1]) +
          this->PlaneNormal[planes][2] * (v_triangle[1][2] - this->PlanePoint[planes][2]);
        p[2] = this->PlaneNormal[planes][0] * (v_triangle[2][0] - this->PlanePoint[planes][0]) +
          this->PlaneNormal[planes][1] * (v_triangle[2][1] - this->PlanePoint[planes][1]) +
          this->PlaneNormal[planes][2] * (v_triangle[2][2] - this->PlanePoint[planes][2]);

        for (int edgeNum = 0; edgeNum < 3; edgeNum++)
        {
          verts = edges[edgeNum];

          p1 = v_triangle[verts[0]];
          p2 = v_triangle[verts[1]];
          double s1 = p[verts[0]];
          double s2 = p[verts[1]];
          if ((s1 * s2) < 0)
          {
            deltaScalar = s2 - s1;

            if (deltaScalar > 0)
            {
              pedg1 = p1;
              pedg2 = p2;
              v1 = verts[0];
              v2 = verts[1];
            }
            else
            {
              pedg1 = p2;
              pedg2 = p1;
              v1 = verts[1];
              v2 = verts[0];
              deltaScalar = -deltaScalar;
              t = s1;
              s1 = s2;
              s2 = t;
            }

            // linear interpolation
            t = (deltaScalar == 0.0 ? 0.0 : (-s1) / deltaScalar);

            for (j = 0; j < 3; j++)
            {
              x[j] = pedg1[j] + t * (pedg2[j] - pedg1[j]);
            }

            // Incorporate point into output and interpolate edge data as necessary
            edges_inter = edges_inter * 10 + (edgeNum + 1);

            if (locator->InsertUniquePoint(x, p_id[num_inter]))
            {
              this->InterpolateEdge(outPD, p_id[num_inter], v_id[v1], v_id[v2], t);
            }

            num_inter++;
          } // if edge intersects value
        }   // for all edges

        if (num_inter == 0)
        {
          unsigned int outside = 0;
          for (i = 0; i < 3; i++)
          {
            if (p[i] > 0) // If only one vertex is outside, so the triangle is outside
            {
              outside = 1; // because there is not intersection.
              break;       // some vertex could be on plane, so you need to test all vertex
            }
          }
          if (outside == 0)
          {
            // else it is possible intersection if other plane
            newCellId = newcellArray->InsertNextCell(3, v_id);
          }
          continue;
        }
        switch (num_inter)
        {
          case 2: // We have one quad and one triangle
                  // i0 gets the index of the triangle point that lies alone on
                  // one side of the plane.
            switch (edges_inter)
            {
              case 12:
                i0 = 1;
                break;
              case 23:
                i0 = 2;
                break;
              case 13:
                i0 = 0;
                break;
              default:
                vtkErrorMacro(<< "Intersection not found: Num_inter = " << num_inter
                              << " Edges_inter = " << edges_inter);
                continue;
            }
            if (p[i0] > 0)
            {
              // The v_triangle[3] is outside box, so
              // the first wedge is outside

              // The v_triangle[3] is outside box, so the quad is outside
              // The Quad is inside: two triangles: (v0,v1,p0) and (p0,p1,v1)
              tab_id[0] = v_id[tab2[i0][0]];
              tab_id[1] = v_id[tab2[i0][1]];
              tab_id[2] = p_id[tab2[i0][2]];
              newCellId = newcellArray->InsertNextCell(3, tab_id);
              tab_id[0] = p_id[tab2[i0][2]];
              tab_id[1] = p_id[tab2[i0][3]];
              tab_id[2] = v_id[tab2[i0][0]];
              newCellId = newcellArray->InsertNextCell(3, tab_id);
            }
            else
            {
              // The Triangle is inside: (v0,p0,p1)
              // The correct winding of the new triangle depends on where the
              // plane intersected the original triangle.
              switch (edges_inter)
              {
                case 12:
                case 23:
                  tab_id[0] = v_id[i0];
                  tab_id[1] = p_id[1];
                  tab_id[2] = p_id[0];
                  break;
                case 13:
                  tab_id[0] = v_id[i0];
                  tab_id[1] = p_id[0];
                  tab_id[2] = p_id[1];
                  break;
              }
              newCellId = newcellArray->InsertNextCell(3, tab_id);
            }
            break;

          case 1: // We have two triangles
            switch (edges_inter)
            {
              case 1:
                i0 = 0;
                break;
              case 2:
                i0 = 1;
                break;
              case 3:
                i0 = 2;
                break;
              default:
                vtkErrorMacro(<< "Intersection not found: Num_inter = " << num_inter
                              << " Edges_inter = " << edges_inter);
                continue;
            }
            if (p[i0] > 0)
            {
              // Isolate vertex is outside box, so the triangle is outside
              tab_id[0] = v_id[tab1[i0][1]];
              tab_id[1] = v_id[tab1[i0][0]];
              tab_id[2] = p_id[0];
              newCellId = newcellArray->InsertNextCell(3, tab_id);
            }
            else
            {
              tab_id[0] = v_id[tab1[i0][0]];
              tab_id[1] = v_id[i0];
              tab_id[2] = p_id[0];
              newCellId = newcellArray->InsertNextCell(3, tab_id);
            }
            break;
        }
      } // for all new cells
      cellarray->Delete();
      cellarray = newcellArray;
    } // for all planes

    unsigned int totalnewcells = cellarray->GetNumberOfCells();

    for (idcellnew = 0; idcellnew < totalnewcells; idcellnew++)
    {
      cellarray->GetNextCell(npts, v_id);
      newCellId = tets->InsertNextCell(npts, v_id);
      outCD->CopyData(inCD, cellId, newCellId);
    }
    cellarray->Delete();
  }
  arraytriangle->Delete();
}

//-------------------------------------------------------

void vtkBoxClipDataSet::ClipHexahedronInOut2D(vtkPoints* newPoints, vtkGenericCell* cell,
  vtkIncrementalPointLocator* locator, vtkCellArray** tets, vtkPointData* inPD,
  vtkPointData** outPD, vtkCellData* inCD, vtkIdType cellId, vtkCellData** outCD)
{
  vtkIdType cellType = cell->GetCellType();
  vtkIdList* cellIds = cell->GetPointIds();
  vtkCellArray* arraytriangle = vtkCellArray::New();
  vtkPoints* cellPts = cell->GetPoints();
  vtkIdType npts = cellPts->GetNumberOfPoints();
  std::vector<vtkIdType> cellptId(npts);
  vtkIdType iid[3];
  const vtkIdType* v_id = nullptr;
  vtkIdType *verts, v1, v2;
  vtkIdType ptId;
  vtkIdType ptIdout[3];
  vtkIdType tab_id[3];
  vtkIdType ptstriangle = 3;

  int i, j, k;
  unsigned int allInside;
  unsigned int idtrianglenew;
  unsigned int idcellnew;
  unsigned int planes;

  /* Edges Triangle */
  vtkIdType edges[3][2] = {
    { 0, 1 },
    { 1, 2 },
    { 2, 0 },
  };
  double deltaScalar;
  double t, *p1, *p2;
  double v[3], x[3];
  double p[6];
  double v_triangle[3][3];

  // To avoid false positive warning.
  tab_id[0] = 0;
  tab_id[1] = 0;
  tab_id[2] = 0;

  for (i = 0; i < npts; i++)
  {
    cellptId[i] = cellIds->GetId(i);
  }

  // Convert all polygon cells to triangles
  this->CellGrid(cellType, npts, cellptId.data(), arraytriangle);

  unsigned int totalnewtriangle = arraytriangle->GetNumberOfCells();
  for (idtrianglenew = 0; idtrianglenew < totalnewtriangle; idtrianglenew++)
  {
    arraytriangle->GetNextCell(ptstriangle, v_id);

    for (allInside = 1, i = 0; i < 3; i++)
    {
      cellPts->GetPoint(v_id[i], v);

      for (k = 0; k < 6; k++)
      {
        p[k] = this->PlaneNormal[k][0] * (v[0] - this->PlanePoint[k][0]) +
          this->PlaneNormal[k][1] * (v[1] - this->PlanePoint[k][1]) +
          this->PlaneNormal[k][2] * (v[2] - this->PlanePoint[k][2]);
      }

      if (!((p[0] <= 0) && (p[1] <= 0) && (p[2] <= 0) && (p[3] <= 0) && (p[4] <= 0) && (p[5] <= 0)))
      {
        allInside = 0;
      }
    } // for all points of the trianglehedron.

    // Test Outside
    unsigned int test[6] = { 1, 1, 1, 1, 1, 1 };
    for (i = 0; i < 3; i++)
    {
      ptIdout[i] = cellIds->GetId(v_id[i]);
      cellPts->GetPoint(v_id[i], v_triangle[i]);

      // Use plane equation
      for (k = 0; k < 6; k++)
      {
        p[k] = this->PlaneNormal[k][0] * (v_triangle[i][0] - this->PlanePoint[k][0]) +
          this->PlaneNormal[k][1] * (v_triangle[i][1] - this->PlanePoint[k][1]) +
          this->PlaneNormal[k][2] * (v_triangle[i][2] - this->PlanePoint[k][2]);
      }

      for (k = 0; k < 3; k++)
      {
        if (p[2 * k] <= 0)
        {
          test[2 * k] = 0;
        }
        if (p[2 * k + 1] <= 0)
        {
          test[2 * k + 1] = 0;
        }
      }
    } // for all points of the cell.

    if (!allInside &&
      ((test[0] == 1) || (test[1] == 1) || (test[2] == 1) || (test[3] == 1) || (test[4] == 1) ||
        (test[5] == 1)))
    {
      for (i = 0; i < 3; i++)
      {
        if (locator->InsertUniquePoint(v_triangle[i], iid[i]))
        {
          outPD[0]->CopyData(inPD, ptIdout[i], iid[i]);
          outPD[1]->CopyData(inPD, ptIdout[i], iid[i]);
        }
      }

      int newCellId = tets[1]->InsertNextCell(3, iid);
      outCD[1]->CopyData(inCD, cellId, newCellId);
      continue; // Triangle is outside.
    }

    for (i = 0; i < 3; i++)
    {
      ptId = cellIds->GetId(v_id[i]);
      cellPts->GetPoint(v_id[i], v);

      // Currently all points are injected because of the possibility
      // of intersection point merging.
      if (locator->InsertUniquePoint(v, iid[i]))
      {
        outPD[0]->CopyData(inPD, ptId, iid[i]);
        outPD[1]->CopyData(inPD, ptId, iid[i]);
      }

    } // for all points of the trianglehedron.

    if (allInside)
    {
      int newCellId = tets[0]->InsertNextCell(3, iid); // Tetrahedron inside.
      outCD[0]->CopyData(inCD, cellId, newCellId);
      continue;
    }

    double *pedg1, *pedg2;

    unsigned int tab2[3][4] = {
      { 1, 2, 1, 0 },
      { 2, 0, 0, 1 },
      { 0, 1, 0, 1 },
    };
    unsigned int tab1[3][2] = {
      { 2, 1 },
      { 0, 2 },
      { 1, 0 },
    };

    vtkCellArray* cellarray = vtkCellArray::New();
    vtkCellArray* cellarrayout = vtkCellArray::New();
    int newCellId = cellarray->InsertNextCell(3, iid);

    // Test Cell intersection with each plane of box
    for (planes = 0; planes < 6; planes++)
    {
      unsigned int totalnewcells = cellarray->GetNumberOfCells();
      vtkCellArray* newcellArray = vtkCellArray::New();

      for (idcellnew = 0; idcellnew < totalnewcells; idcellnew++)
      {
        unsigned int i0;
        unsigned int num_inter = 0;
        unsigned int edges_inter = 0;
        vtkIdType p_id[3];

        cellarray->GetNextCell(npts, v_id);

        newPoints->GetPoint(v_id[0], v_triangle[0]); // coord (x,y,z)
        newPoints->GetPoint(v_id[1], v_triangle[1]);
        newPoints->GetPoint(v_id[2], v_triangle[2]);

        p[0] = this->PlaneNormal[planes][0] * (v_triangle[0][0] - this->PlanePoint[planes][0]) +
          this->PlaneNormal[planes][1] * (v_triangle[0][1] - this->PlanePoint[planes][1]) +
          this->PlaneNormal[planes][2] * (v_triangle[0][2] - this->PlanePoint[planes][2]);
        p[1] = this->PlaneNormal[planes][0] * (v_triangle[1][0] - this->PlanePoint[planes][0]) +
          this->PlaneNormal[planes][1] * (v_triangle[1][1] - this->PlanePoint[planes][1]) +
          this->PlaneNormal[planes][2] * (v_triangle[1][2] - this->PlanePoint[planes][2]);
        p[2] = this->PlaneNormal[planes][0] * (v_triangle[2][0] - this->PlanePoint[planes][0]) +
          this->PlaneNormal[planes][1] * (v_triangle[2][1] - this->PlanePoint[planes][1]) +
          this->PlaneNormal[planes][2] * (v_triangle[2][2] - this->PlanePoint[planes][2]);

        for (int edgeNum = 0; edgeNum < 3; edgeNum++)
        {
          verts = edges[edgeNum];

          p1 = v_triangle[verts[0]];
          p2 = v_triangle[verts[1]];
          double s1 = p[verts[0]];
          double s2 = p[verts[1]];
          if ((s1 * s2) < 0)
          {
            deltaScalar = s2 - s1;

            if (deltaScalar > 0)
            {
              pedg1 = p1;
              pedg2 = p2;
              v1 = verts[0];
              v2 = verts[1];
            }
            else
            {
              pedg1 = p2;
              pedg2 = p1;
              v1 = verts[1];
              v2 = verts[0];
              deltaScalar = -deltaScalar;
              t = s1;
              s1 = s2;
              s2 = t;
            }

            // linear interpolation
            t = (deltaScalar == 0.0 ? 0.0 : (-s1) / deltaScalar);

            for (j = 0; j < 3; j++)
            {
              x[j] = pedg1[j] + t * (pedg2[j] - pedg1[j]);
            }

            // Incorporate point into output and interpolate edge data as necessary
            edges_inter = edges_inter * 10 + (edgeNum + 1);

            if (locator->InsertUniquePoint(x, p_id[num_inter]))
            {
              this->InterpolateEdge(outPD[0], p_id[num_inter], v_id[v1], v_id[v2], t);
              this->InterpolateEdge(outPD[1], p_id[num_inter], v_id[v1], v_id[v2], t);
            }

            num_inter++;
          } // if edge intersects value
        }   // for all edges

        if (num_inter == 0)
        {
          unsigned int outside = 0;
          for (i = 0; i < 3; i++)
          {
            if (p[i] > 0) // If only one vertex is outside, so the trianglehedron is outside
            {
              outside = 1; // because there is not intersection.
              break;       // some vertex could be on plane, so you need to test all vertex
            }
          }
          if (outside == 0)
          {
            // else it is possible intersection if other plane
            newCellId = newcellArray->InsertNextCell(3, v_id);
          }
          else
          {
            newCellId = tets[1]->InsertNextCell(3, v_id);
            outCD[1]->CopyData(inCD, cellId, newCellId);
          }
          continue;
        }
        switch (num_inter)
        {
          case 2: // We have one quad and one triangle
            switch (edges_inter)
            {
              case 12:
                i0 = 1;
                break;
              case 23:
                i0 = 2;
                break;
              case 13:
                i0 = 0;
                break;
              default:
                vtkErrorMacro(<< "Intersection not found: Num_inter = " << num_inter
                              << " Edges_inter = " << edges_inter);
                continue;
            }

            if (p[i0] > 0)
            {
              // The v_triangle[3] is outside box, so
              // the first wedge is outside

              // The Quad is inside: two triangles: (v0,v1,p0) and (p0,p1,v1)
              tab_id[0] = v_id[tab2[i0][0]];
              tab_id[1] = v_id[tab2[i0][1]];
              tab_id[2] = p_id[tab2[i0][2]];
              newCellId = newcellArray->InsertNextCell(3, tab_id);
              tab_id[0] = p_id[tab2[i0][2]];
              tab_id[1] = p_id[tab2[i0][3]];
              tab_id[2] = v_id[tab2[i0][0]];
              newCellId = newcellArray->InsertNextCell(3, tab_id);

              switch (edges_inter) // Triangle Outside
              {
                case 12:
                case 23:
                  tab_id[0] = v_id[i0];
                  tab_id[1] = p_id[1];
                  tab_id[2] = p_id[0];
                  break;
                case 13:
                  tab_id[0] = v_id[i0];
                  tab_id[1] = p_id[0];
                  tab_id[2] = p_id[1];
                  break;
              }
              newCellId = cellarrayout->InsertNextCell(3, tab_id);
            }
            else
            {
              // The Triangle is inside: (v0,p0,p1)
              switch (edges_inter)
              {
                case 12:
                case 23:
                  tab_id[0] = v_id[i0];
                  tab_id[1] = p_id[1];
                  tab_id[2] = p_id[0];
                  break;
                case 13:
                  tab_id[0] = v_id[i0];
                  tab_id[1] = p_id[0];
                  tab_id[2] = p_id[1];
                  break;
              }
              newCellId = newcellArray->InsertNextCell(3, tab_id);

              // The Quad is outside: two triangles: (v0,v1,p0) and (p0,p1,v1)
              tab_id[0] = v_id[tab2[i0][0]];
              tab_id[1] = v_id[tab2[i0][1]];
              tab_id[2] = p_id[tab2[i0][2]];
              newCellId = cellarrayout->InsertNextCell(3, tab_id);
              tab_id[0] = p_id[tab2[i0][2]];
              tab_id[1] = p_id[tab2[i0][3]];
              tab_id[2] = v_id[tab2[i0][0]];
              newCellId = cellarrayout->InsertNextCell(3, tab_id);
            }
            break;

          case 1: // We have two triangles
            switch (edges_inter)
            {
              case 1:
                i0 = 0;
                break;
              case 2:
                i0 = 1;
                break;
              case 3:
                i0 = 2;
                break;
              default:
                vtkErrorMacro(<< "Intersection not found: Num_inter = " << num_inter
                              << " Edges_inter = " << edges_inter);
                continue;
            }
            if (p[i0] > 0)
            {
              // Isolate vertex is outside box, so
              // the triangle is outside
              tab_id[0] = v_id[tab1[i0][1]]; // Inside
              tab_id[1] = v_id[tab1[i0][0]];
              tab_id[2] = p_id[0];
              newCellId = newcellArray->InsertNextCell(3, tab_id);

              tab_id[0] = v_id[tab1[i0][0]]; // Outside
              tab_id[1] = v_id[i0];
              tab_id[2] = p_id[0];
              newCellId = cellarrayout->InsertNextCell(3, tab_id);
            }
            else
            {
              tab_id[0] = v_id[tab1[i0][0]]; // Inside
              tab_id[1] = v_id[i0];
              tab_id[2] = p_id[0];
              newCellId = newcellArray->InsertNextCell(3, tab_id);

              tab_id[0] = v_id[tab1[i0][1]]; // Outside
              tab_id[1] = v_id[tab1[i0][0]];
              tab_id[2] = p_id[0];
              newCellId = cellarrayout->InsertNextCell(3, tab_id);
            }
            break;
        }
      } // for all new cells
      cellarray->Delete();
      cellarray = newcellArray;
    } // for all planes

    unsigned int totalnewcells = cellarray->GetNumberOfCells(); // Inside

    for (idcellnew = 0; idcellnew < totalnewcells; idcellnew++)
    {
      cellarray->GetNextCell(npts, v_id);
      newCellId = tets[0]->InsertNextCell(npts, v_id);
      outCD[0]->CopyData(inCD, cellId, newCellId);
    }
    cellarray->Delete();

    totalnewcells = cellarrayout->GetNumberOfCells(); // Outside

    for (idcellnew = 0; idcellnew < totalnewcells; idcellnew++)
    {
      cellarrayout->GetNextCell(npts, v_id);
      newCellId = tets[1]->InsertNextCell(npts, v_id);
      outCD[1]->CopyData(inCD, cellId, newCellId);
    }
    cellarrayout->Delete();
  }
  arraytriangle->Delete();
}

//------------------------------------------------------------------------------

void vtkBoxClipDataSet::ClipBox1D(vtkPoints* newPoints, vtkGenericCell* cell,
  vtkIncrementalPointLocator* locator, vtkCellArray* lines, vtkPointData* inPD, vtkPointData* outPD,
  vtkCellData* inCD, vtkIdType cellId, vtkCellData* outCD)
{
  vtkIdType cellType = cell->GetCellType();
  vtkIdList* cellIds = cell->GetPointIds();
  vtkCellArray* arrayline = vtkCellArray::New();
  vtkPoints* cellPts = cell->GetPoints();
  vtkIdType npts = cellPts->GetNumberOfPoints();
  std::vector<vtkIdType> cellptId(npts);
  vtkIdType iid[2];
  const vtkIdType* v_id = nullptr;
  vtkIdType ptId;
  vtkIdType tab_id[2];
  vtkIdType ptsline = 2;

  int i, j;
  unsigned int allInside;
  unsigned int planes;
  unsigned int cutInd;

  double value;
  double t;
  double v[3], x[3];
  double v_line[2][3];

  for (i = 0; i < npts; i++)
  {
    cellptId[i] = cellIds->GetId(i);
  }

  // Convert all 1d cells to single line.
  this->CellGrid(cellType, npts, cellptId.data(), arrayline);

  unsigned int totalnewline = arrayline->GetNumberOfCells();
  for (unsigned int idlinenew = 0; idlinenew < totalnewline; idlinenew++)
  {
    arrayline->GetNextCell(ptsline, v_id);

    for (allInside = 1, i = 0; i < 2; i++)
    {
      cellPts->GetPoint(v_id[i], v);

      if (!((v[0] >= this->BoundBoxClip[0][0]) && (v[0] <= this->BoundBoxClip[0][1]) &&
            (v[1] >= this->BoundBoxClip[1][0]) && (v[1] <= this->BoundBoxClip[1][1]) &&
            (v[2] >= this->BoundBoxClip[2][0]) && (v[2] <= this->BoundBoxClip[2][1])))
      {
        // outside
        allInside = 0;
      }
    }

    // Test Outside:
    if (!allInside)
    {
      unsigned int test[6] = { 1, 1, 1, 1, 1, 1 };
      for (i = 0; i < 2; i++)
      {
        cellPts->GetPoint(v_id[i], v);

        if (v[0] >= this->BoundBoxClip[0][0])
        {
          test[0] = 0;
        }
        if (v[0] <= this->BoundBoxClip[0][1])
        {
          test[1] = 0;
        }
        if (v[1] >= this->BoundBoxClip[1][0])
        {
          test[2] = 0;
        }
        if (v[1] <= this->BoundBoxClip[1][1])
        {
          test[3] = 0;
        }
        if (v[2] >= this->BoundBoxClip[2][0])
        {
          test[4] = 0;
        }
        if (v[2] <= this->BoundBoxClip[2][1])
        {
          test[5] = 0;
        }
      } // for all points of the line.

      if ((test[0] == 1) || (test[1] == 1) || (test[2] == 1) || (test[3] == 1) || (test[4] == 1) ||
        (test[5] == 1))
      {
        continue; // Line is outside.
      }
    } // if not allInside

    for (i = 0; i < 2; i++)
    {
      // Currently all points are injected because of the possibility
      // of intersection point merging.
      ptId = cellIds->GetId(v_id[i]);
      cellPts->GetPoint(v_id[i], v);
      if (locator->InsertUniquePoint(v, iid[i]))
      {
        outPD->CopyData(inPD, ptId, iid[i]);
      }
    } // for all points of the triangle.

    if (allInside)
    {
      // Triangle inside.
      int newCellId = lines->InsertNextCell(2, iid);
      outCD->CopyData(inCD, cellId, newCellId);
      continue;
    }

    vtkCellArray* cellarray = vtkCellArray::New();
    int newCellId = cellarray->InsertNextCell(2, iid);
    unsigned int idcellnew;

    // Test triangle intersection with each plane of box
    for (planes = 0; planes < 6; planes++)
    {
      // The index of the dimension of the cut plane (x == 0, y == 1, z == 2).
      cutInd = planes / 2;

      // The plane is always parallel to unitary cube.
      value = this->BoundBoxClip[cutInd][planes % 2];

      unsigned int totalnewcells = cellarray->GetNumberOfCells();
      vtkCellArray* newcellArray = vtkCellArray::New();

      for (idcellnew = 0; idcellnew < totalnewcells; idcellnew++)
      {
        vtkIdType p_id;
        cellarray->GetNextCell(npts, v_id);

        newPoints->GetPoint(v_id[0], v_line[0]); // coord (x,y,z)
        newPoints->GetPoint(v_id[1], v_line[1]);

        // Check to see if line is inside plane.
        if (((planes % 2 == 0) && (v_line[0][cutInd] >= value) && (v_line[1][cutInd] >= value)) ||
          ((planes % 2 == 1) && (v_line[0][cutInd] <= value) && (v_line[1][cutInd] <= value)))
        {
          newCellId = newcellArray->InsertNextCell(2, v_id);
          continue;
        }

        // Check to see if line is outside plane.
        if (((planes % 2 == 0) && (v_line[0][cutInd] <= value) && (v_line[1][cutInd] <= value)) ||
          ((planes % 2 == 1) && (v_line[0][cutInd] >= value) && (v_line[1][cutInd] >= value)))
        {
          continue;
        }

        // If we are here, the plane intersects the line segment.
        t = (value - v_line[0][cutInd]) / (v_line[1][cutInd] - v_line[0][cutInd]);
        for (j = 0; j < 3; j++)
        {
          x[j] = (v_line[1][j] - v_line[0][j]) * t + v_line[0][j];
        }

        // Incorporate point into output and interpolate edge data as
        // necessary.
        if (locator->InsertUniquePoint(x, p_id))
        {
          this->InterpolateEdge(outPD, p_id, v_id[0], v_id[1], t);
        }

        // Add the clipped line to the output.
        if (((planes % 2 == 0) && (v_line[0][cutInd] >= value)) ||
          ((planes % 2 == 1) && (v_line[0][cutInd] <= value)))
        {
          // First point of line is inside.
          tab_id[0] = v_id[0];
          tab_id[1] = p_id;
          newcellArray->InsertNextCell(2, tab_id);
        }
        else
        {
          // Second point of line is inside.
          tab_id[0] = p_id;
          tab_id[1] = v_id[1];
          newcellArray->InsertNextCell(2, tab_id);
        }
      } // for all new cells
      cellarray->Delete();
      cellarray = newcellArray;
    } // for all planes

    unsigned int totalnewcells = cellarray->GetNumberOfCells();

    for (idcellnew = 0; idcellnew < totalnewcells; idcellnew++)
    {
      cellarray->GetNextCell(npts, v_id);
      newCellId = lines->InsertNextCell(npts, v_id);
      outCD->CopyData(inCD, cellId, newCellId);
    }
    cellarray->Delete();
  }
  arrayline->Delete();
}

//------------------------------------------------------------------------------

void vtkBoxClipDataSet::ClipBoxInOut1D(vtkPoints* newPoints, vtkGenericCell* cell,
  vtkIncrementalPointLocator* locator, vtkCellArray** lines, vtkPointData* inPD,
  vtkPointData** outPD, vtkCellData* inCD, vtkIdType cellId, vtkCellData** outCD)
{
  vtkIdType cellType = cell->GetCellType();
  vtkIdList* cellIds = cell->GetPointIds();
  vtkCellArray* arrayline = vtkCellArray::New();
  vtkPoints* cellPts = cell->GetPoints();
  vtkIdType npts = cellPts->GetNumberOfPoints();
  std::vector<vtkIdType> cellptId(npts);
  vtkIdType iid[2];
  const vtkIdType* v_id = nullptr;
  vtkIdType ptId;
  vtkIdType tab_id[2];
  vtkIdType ptsline = 2;

  int i, j;
  unsigned int allInside;
  unsigned int planes;
  unsigned int cutInd;

  double value;
  double t;
  double v[3], x[3];
  double v_line[2][3];

  for (i = 0; i < npts; i++)
  {
    cellptId[i] = cellIds->GetId(i);
  }

  // Convert all 1d cells to single line.
  this->CellGrid(cellType, npts, cellptId.data(), arrayline);

  unsigned int totalnewline = arrayline->GetNumberOfCells();
  for (unsigned int idlinenew = 0; idlinenew < totalnewline; idlinenew++)
  {
    arrayline->GetNextCell(ptsline, v_id);

    for (allInside = 1, i = 0; i < 2; i++)
    {
      ptId = cellIds->GetId(v_id[i]);
      cellPts->GetPoint(v_id[i], v);

      if (!((v[0] >= this->BoundBoxClip[0][0]) && (v[0] <= this->BoundBoxClip[0][1]) &&
            (v[1] >= this->BoundBoxClip[1][0]) && (v[1] <= this->BoundBoxClip[1][1]) &&
            (v[2] >= this->BoundBoxClip[2][0]) && (v[2] <= this->BoundBoxClip[2][1])))
      {
        // outside
        allInside = 0;
      }

      if (locator->InsertUniquePoint(v, iid[i]))
      {
        outPD[0]->CopyData(inPD, ptId, iid[i]);
        outPD[1]->CopyData(inPD, ptId, iid[i]);
      }
    }

    // Test Outside:
    if (!allInside)
    {
      unsigned int test[6] = { 1, 1, 1, 1, 1, 1 };
      for (i = 0; i < 2; i++)
      {
        cellPts->GetPoint(v_id[i], v);

        if (v[0] >= this->BoundBoxClip[0][0])
        {
          test[0] = 0;
        }
        if (v[0] <= this->BoundBoxClip[0][1])
        {
          test[1] = 0;
        }
        if (v[1] >= this->BoundBoxClip[1][0])
        {
          test[2] = 0;
        }
        if (v[1] <= this->BoundBoxClip[1][1])
        {
          test[3] = 0;
        }
        if (v[2] >= this->BoundBoxClip[2][0])
        {
          test[4] = 0;
        }
        if (v[2] <= this->BoundBoxClip[2][1])
        {
          test[5] = 0;
        }
      } // for all points of the line.

      if ((test[0] == 1) || (test[1] == 1) || (test[2] == 1) || (test[3] == 1) || (test[4] == 1) ||
        (test[5] == 1))
      {
        int newCellId = lines[1]->InsertNextCell(2, iid);
        outCD[1]->CopyData(inCD, cellId, newCellId);
        continue; // Line is outside.
      }
    } // if not allInside

    if (allInside)
    {
      // Triangle inside.
      int newCellId = lines[0]->InsertNextCell(2, iid);
      outCD[0]->CopyData(inCD, cellId, newCellId);
      continue;
    }

    vtkCellArray* cellarray = vtkCellArray::New();
    vtkCellArray* cellarrayout = vtkCellArray::New();
    int newCellId = cellarray->InsertNextCell(2, iid);
    unsigned int idcellnew;

    // Test triangle intersection with each plane of box
    for (planes = 0; planes < 6; planes++)
    {
      // The index of the dimension of the cut plane (x == 0, y == 1, z == 2).
      cutInd = planes / 2;

      // The plane is always parallel to unitary cube.
      value = this->BoundBoxClip[cutInd][planes % 2];

      unsigned int totalnewcells = cellarray->GetNumberOfCells();
      vtkCellArray* newcellArray = vtkCellArray::New();

      for (idcellnew = 0; idcellnew < totalnewcells; idcellnew++)
      {
        vtkIdType p_id;
        cellarray->GetNextCell(npts, v_id);

        newPoints->GetPoint(v_id[0], v_line[0]); // coord (x,y,z)
        newPoints->GetPoint(v_id[1], v_line[1]);

        // Check to see if line is inside plane.
        if (((planes % 2 == 0) && (v_line[0][cutInd] >= value) && (v_line[1][cutInd] >= value)) ||
          ((planes % 2 == 1) && (v_line[0][cutInd] <= value) && (v_line[1][cutInd] <= value)))
        {
          newCellId = newcellArray->InsertNextCell(2, v_id);
          continue;
        }

        // Check to see if line is outside plane.
        if (((planes % 2 == 0) && (v_line[0][cutInd] <= value) && (v_line[1][cutInd] <= value)) ||
          ((planes % 2 == 1) && (v_line[0][cutInd] >= value) && (v_line[1][cutInd] >= value)))
        {
          newCellId = lines[1]->InsertNextCell(2, v_id);
          outCD[1]->CopyData(inCD, cellId, newCellId);
          continue;
        }

        // If we are here, the plane intersects the line segment.
        t = (value - v_line[0][cutInd]) / (v_line[1][cutInd] - v_line[0][cutInd]);
        for (j = 0; j < 3; j++)
        {
          x[j] = (v_line[1][j] - v_line[0][j]) * t + v_line[0][j];
        }

        // Incorporate point into output and interpolate edge data as
        // necessary.
        if (locator->InsertUniquePoint(x, p_id))
        {
          this->InterpolateEdge(outPD[0], p_id, v_id[0], v_id[1], t);
          this->InterpolateEdge(outPD[1], p_id, v_id[0], v_id[1], t);
        }

        // Add the clipped line to the output.
        if (((planes % 2 == 0) && (v_line[0][cutInd] >= value)) ||
          ((planes % 2 == 1) && (v_line[0][cutInd] <= value)))
        {
          // First point of line is inside.
          tab_id[0] = v_id[0];
          tab_id[1] = p_id;
          newcellArray->InsertNextCell(2, tab_id);

          // Second point of line is outside.
          tab_id[0] = p_id;
          tab_id[1] = v_id[1];
          cellarrayout->InsertNextCell(2, tab_id);
        }
        else
        {
          // Second point of line is inside.
          tab_id[0] = p_id;
          tab_id[1] = v_id[1];
          newcellArray->InsertNextCell(2, tab_id);

          // First point of line is outside.
          tab_id[0] = v_id[0];
          tab_id[1] = p_id;
          cellarrayout->InsertNextCell(2, tab_id);
        }
      } // for all new cells
      cellarray->Delete();
      cellarray = newcellArray;
    } // for all planes

    unsigned int totalnewcells = cellarray->GetNumberOfCells(); // Inside
    for (idcellnew = 0; idcellnew < totalnewcells; idcellnew++)
    {
      cellarray->GetNextCell(npts, v_id);
      newCellId = lines[0]->InsertNextCell(npts, v_id);
      outCD[0]->CopyData(inCD, cellId, newCellId);
    }
    cellarray->Delete();

    totalnewcells = cellarrayout->GetNumberOfCells(); // Outside
    for (idcellnew = 0; idcellnew < totalnewcells; idcellnew++)
    {
      cellarrayout->GetNextCell(npts, v_id);
      newCellId = lines[1]->InsertNextCell(npts, v_id);
      outCD[1]->CopyData(inCD, cellId, newCellId);
    }
    cellarrayout->Delete();
  }
  arrayline->Delete();
}

//------------------------------------------------------------------------------

void vtkBoxClipDataSet::ClipHexahedron1D(vtkPoints* newPoints, vtkGenericCell* cell,
  vtkIncrementalPointLocator* locator, vtkCellArray* lines, vtkPointData* inPD, vtkPointData* outPD,
  vtkCellData* inCD, vtkIdType cellId, vtkCellData* outCD)
{
  vtkIdType cellType = cell->GetCellType();
  vtkIdList* cellIds = cell->GetPointIds();
  vtkCellArray* arrayline = vtkCellArray::New();
  vtkPoints* cellPts = cell->GetPoints();
  vtkIdType npts = cellPts->GetNumberOfPoints();
  std::vector<vtkIdType> cellptId(npts);
  vtkIdType iid[2];
  const vtkIdType* v_id = nullptr;
  vtkIdType ptId;
  vtkIdType tab_id[2];
  vtkIdType ptsline = 2;

  int i, j;
  unsigned int allInside;
  unsigned int planes;

  double values[2];
  double t;
  double v[3], x[3];
  double v_line[2][3];

  for (i = 0; i < npts; i++)
  {
    cellptId[i] = cellIds->GetId(i);
  }

  // Convert all 1d cells to single line.
  this->CellGrid(cellType, npts, cellptId.data(), arrayline);

  unsigned int totalnewline = arrayline->GetNumberOfCells();
  for (unsigned int idlinenew = 0; idlinenew < totalnewline; idlinenew++)
  {
    arrayline->GetNextCell(ptsline, v_id);

    for (allInside = 1, i = 0; i < 2; i++)
    {
      cellPts->GetPoint(v_id[i], v);

      for (int k = 0; k < 6; k++)
      {
        values[0] = (this->PlaneNormal[k][0] * (v[0] - this->PlanePoint[k][0]) +
          this->PlaneNormal[k][1] * (v[1] - this->PlanePoint[k][1]) +
          this->PlaneNormal[k][2] * (v[2] - this->PlanePoint[k][2]));
        if (values[0] > 0)
        {
          // outside
          allInside = 0;
        }
      }
    }

    // Test Outside:
    if (!allInside)
    {
      unsigned int test[6] = { 1, 1, 1, 1, 1, 1 };
      for (i = 0; i < 2; i++)
      {
        cellPts->GetPoint(v_id[i], v);

        // Use plane equation.
        for (int k = 0; k < 6; k++)
        {
          values[0] = (this->PlaneNormal[k][0] * (v[0] - this->PlanePoint[k][0]) +
            this->PlaneNormal[k][1] * (v[1] - this->PlanePoint[k][1]) +
            this->PlaneNormal[k][2] * (v[2] - this->PlanePoint[k][2]));
          if (values[0] <= 0)
          {
            test[k] = 0;
          }
        } // for all planes of the hexahedron.
      }   // for all points of the line.

      if ((test[0] == 1) || (test[1] == 1) || (test[2] == 1) || (test[3] == 1) || (test[4] == 1) ||
        (test[5] == 1))
      {
        continue; // Line is outside.
      }
    } // if not allInside

    for (i = 0; i < 2; i++)
    {
      // Currently all points are injected because of the possibility
      // of intersection point merging.
      ptId = cellIds->GetId(v_id[i]);
      cellPts->GetPoint(v_id[i], v);
      if (locator->InsertUniquePoint(v, iid[i]))
      {
        outPD->CopyData(inPD, ptId, iid[i]);
      }
    } // for all points of the triangle.

    if (allInside)
    {
      // Triangle inside.
      int newCellId = lines->InsertNextCell(2, iid);
      outCD->CopyData(inCD, cellId, newCellId);
      continue;
    }

    vtkCellArray* cellarray = vtkCellArray::New();
    int newCellId = cellarray->InsertNextCell(2, iid);
    unsigned int idcellnew;

    // Test triangle intersection with each plane of box
    for (planes = 0; planes < 6; planes++)
    {
      unsigned int totalnewcells = cellarray->GetNumberOfCells();
      vtkCellArray* newcellArray = vtkCellArray::New();

      for (idcellnew = 0; idcellnew < totalnewcells; idcellnew++)
      {
        vtkIdType p_id;
        cellarray->GetNextCell(npts, v_id);

        newPoints->GetPoint(v_id[0], v_line[0]); // coord (x,y,z)
        newPoints->GetPoint(v_id[1], v_line[1]);

        const double* planeNormal = this->PlaneNormal[planes];
        const double* planePoint = this->PlanePoint[planes];
        values[0] = (planeNormal[0] * (v_line[0][0] - planePoint[0]) +
          planeNormal[1] * (v_line[0][1] - planePoint[1]) +
          planeNormal[2] * (v_line[0][2] - planePoint[2]));
        values[1] = (planeNormal[0] * (v_line[1][0] - planePoint[0]) +
          planeNormal[1] * (v_line[1][1] - planePoint[1]) +
          planeNormal[2] * (v_line[1][2] - planePoint[2]));

        // Check to see if line is inside plane.
        if ((values[0] <= 0) && (values[1] <= 0))
        {
          newCellId = newcellArray->InsertNextCell(2, v_id);
          continue;
        }

        // Check to see if line is outside plane.
        if ((values[0] >= 0) && (values[1] >= 0))
        {
          continue;
        }

        // If we are here, the plane intersects the line segment.
        t = values[0] / (values[0] - values[1]);
        for (j = 0; j < 3; j++)
        {
          x[j] = (v_line[1][j] - v_line[0][j]) * t + v_line[0][j];
        }

        // Incorporate point into output and interpolate edge data as
        // necessary.
        if (locator->InsertUniquePoint(x, p_id))
        {
          this->InterpolateEdge(outPD, p_id, v_id[0], v_id[1], t);
        }

        // Add the clipped line to the output.
        if (values[0] <= 0)
        {
          // First point of line is inside.
          tab_id[0] = v_id[0];
          tab_id[1] = p_id;
          newcellArray->InsertNextCell(2, tab_id);
        }
        else
        {
          // Second point of line is inside.
          tab_id[0] = p_id;
          tab_id[1] = v_id[1];
          newcellArray->InsertNextCell(2, tab_id);
        }
      } // for all new cells
      cellarray->Delete();
      cellarray = newcellArray;
    } // for all planes

    unsigned int totalnewcells = cellarray->GetNumberOfCells();

    for (idcellnew = 0; idcellnew < totalnewcells; idcellnew++)
    {
      cellarray->GetNextCell(npts, v_id);
      newCellId = lines->InsertNextCell(npts, v_id);
      outCD->CopyData(inCD, cellId, newCellId);
    }
    cellarray->Delete();
  }
  arrayline->Delete();
}

//------------------------------------------------------------------------------

void vtkBoxClipDataSet::ClipHexahedronInOut1D(vtkPoints* newPoints, vtkGenericCell* cell,
  vtkIncrementalPointLocator* locator, vtkCellArray** lines, vtkPointData* inPD,
  vtkPointData** outPD, vtkCellData* inCD, vtkIdType cellId, vtkCellData** outCD)
{
  vtkIdType cellType = cell->GetCellType();
  vtkIdList* cellIds = cell->GetPointIds();
  vtkCellArray* arrayline = vtkCellArray::New();
  vtkPoints* cellPts = cell->GetPoints();
  vtkIdType npts = cellPts->GetNumberOfPoints();
  std::vector<vtkIdType> cellptId(npts);
  vtkIdType iid[2];
  const vtkIdType* v_id = nullptr;
  vtkIdType ptId;
  vtkIdType tab_id[2];
  vtkIdType ptsline = 2;

  int i, j;
  unsigned int allInside;
  unsigned int planes;

  double values[2];
  double t;
  double v[3], x[3];
  double v_line[2][3];

  for (i = 0; i < npts; i++)
  {
    cellptId[i] = cellIds->GetId(i);
  }

  // Convert all 1d cells to single line.
  this->CellGrid(cellType, npts, cellptId.data(), arrayline);

  unsigned int totalnewline = arrayline->GetNumberOfCells();
  for (unsigned int idlinenew = 0; idlinenew < totalnewline; idlinenew++)
  {
    arrayline->GetNextCell(ptsline, v_id);

    for (allInside = 1, i = 0; i < 2; i++)
    {
      ptId = cellIds->GetId(v_id[i]);
      cellPts->GetPoint(v_id[i], v);

      for (int k = 0; k < 6; k++)
      {
        values[0] = (this->PlaneNormal[k][0] * (v[0] - this->PlanePoint[k][0]) +
          this->PlaneNormal[k][1] * (v[1] - this->PlanePoint[k][1]) +
          this->PlaneNormal[k][2] * (v[2] - this->PlanePoint[k][2]));
        if (values[0] > 0)
        {
          // outside
          allInside = 0;
        }
      }

      if (locator->InsertUniquePoint(v, iid[i]))
      {
        outPD[0]->CopyData(inPD, ptId, iid[i]);
        outPD[1]->CopyData(inPD, ptId, iid[i]);
      }
    }

    // Test Outside:
    if (!allInside)
    {
      unsigned int test[6] = { 1, 1, 1, 1, 1, 1 };
      for (i = 0; i < 2; i++)
      {
        cellPts->GetPoint(v_id[i], v);

        // Use plane equation.
        for (int k = 0; k < 6; k++)
        {
          values[0] = (this->PlaneNormal[k][0] * (v[0] - this->PlanePoint[k][0]) +
            this->PlaneNormal[k][1] * (v[1] - this->PlanePoint[k][1]) +
            this->PlaneNormal[k][2] * (v[2] - this->PlanePoint[k][2]));
          if (values[0] <= 0)
          {
            test[k] = 0;
          }
        } // for all planes of the hexahedron.
      }   // for all points of the line.

      if ((test[0] == 1) || (test[1] == 1) || (test[2] == 1) || (test[3] == 1) || (test[4] == 1) ||
        (test[5] == 1))
      {
        int newCellId = lines[1]->InsertNextCell(2, iid);
        outCD[1]->CopyData(inCD, cellId, newCellId);
        continue; // Line is outside.
      }
    } // if not allInside

    if (allInside)
    {
      // Triangle inside.
      int newCellId = lines[0]->InsertNextCell(2, iid);
      outCD[0]->CopyData(inCD, cellId, newCellId);
      continue;
    }

    vtkCellArray* cellarray = vtkCellArray::New();
    vtkCellArray* cellarrayout = vtkCellArray::New();
    int newCellId = cellarray->InsertNextCell(2, iid);
    unsigned int idcellnew;

    // Test triangle intersection with each plane of box
    for (planes = 0; planes < 6; planes++)
    {
      unsigned int totalnewcells = cellarray->GetNumberOfCells();
      vtkCellArray* newcellArray = vtkCellArray::New();

      for (idcellnew = 0; idcellnew < totalnewcells; idcellnew++)
      {
        vtkIdType p_id;
        cellarray->GetNextCell(npts, v_id);

        newPoints->GetPoint(v_id[0], v_line[0]); // coord (x,y,z)
        newPoints->GetPoint(v_id[1], v_line[1]);

        const double* planeNormal = this->PlaneNormal[planes];
        const double* planePoint = this->PlanePoint[planes];
        values[0] = (planeNormal[0] * (v_line[0][0] - planePoint[0]) +
          planeNormal[1] * (v_line[0][1] - planePoint[1]) +
          planeNormal[2] * (v_line[0][2] - planePoint[2]));
        values[1] = (planeNormal[0] * (v_line[1][0] - planePoint[0]) +
          planeNormal[1] * (v_line[1][1] - planePoint[1]) +
          planeNormal[2] * (v_line[1][2] - planePoint[2]));

        // Check to see if line is inside plane.
        if ((values[0] <= 0) && (values[1] <= 0))
        {
          newCellId = newcellArray->InsertNextCell(2, v_id);
          continue;
        }

        // Check to see if line is outside plane.
        if ((values[0] >= 0) && (values[1] >= 0))
        {
          newCellId = lines[1]->InsertNextCell(2, v_id);
          outCD[1]->CopyData(inCD, cellId, newCellId);
          continue;
        }

        // If we are here, the plane intersects the line segment.
        t = values[0] / (values[0] - values[1]);
        for (j = 0; j < 3; j++)
        {
          x[j] = (v_line[1][j] - v_line[0][j]) * t + v_line[0][j];
        }

        // Incorporate point into output and interpolate edge data as
        // necessary.
        if (locator->InsertUniquePoint(x, p_id))
        {
          this->InterpolateEdge(outPD[0], p_id, v_id[0], v_id[1], t);
          this->InterpolateEdge(outPD[1], p_id, v_id[0], v_id[1], t);
        }

        // Add the clipped line to the output.
        if (values[0] <= 0)
        {
          // First point of line is inside.
          tab_id[0] = v_id[0];
          tab_id[1] = p_id;
          newcellArray->InsertNextCell(2, tab_id);

          // Second point of line is outside.
          tab_id[0] = p_id;
          tab_id[1] = v_id[1];
          cellarrayout->InsertNextCell(2, tab_id);
        }
        else
        {
          // Second point of line is inside.
          tab_id[0] = p_id;
          tab_id[1] = v_id[1];
          newcellArray->InsertNextCell(2, tab_id);

          // First point of line is outside.
          tab_id[0] = v_id[0];
          tab_id[1] = p_id;
          cellarrayout->InsertNextCell(2, tab_id);
        }
      } // for all new cells
      cellarray->Delete();
      cellarray = newcellArray;
    } // for all planes

    unsigned int totalnewcells = cellarray->GetNumberOfCells(); // Inside
    for (idcellnew = 0; idcellnew < totalnewcells; idcellnew++)
    {
      cellarray->GetNextCell(npts, v_id);
      newCellId = lines[0]->InsertNextCell(npts, v_id);
      outCD[0]->CopyData(inCD, cellId, newCellId);
    }
    cellarray->Delete();

    totalnewcells = cellarrayout->GetNumberOfCells(); // Outside
    for (idcellnew = 0; idcellnew < totalnewcells; idcellnew++)
    {
      cellarrayout->GetNextCell(npts, v_id);
      newCellId = lines[1]->InsertNextCell(npts, v_id);
      outCD[1]->CopyData(inCD, cellId, newCellId);
    }
    cellarrayout->Delete();
  }
  arrayline->Delete();
}

//------------------------------------------------------------------------------

void vtkBoxClipDataSet::ClipBox0D(vtkGenericCell* cell, vtkIncrementalPointLocator* locator,
  vtkCellArray* verts, vtkPointData* inPD, vtkPointData* outPD, vtkCellData* inCD, vtkIdType cellId,
  vtkCellData* outCD)
{
  vtkIdType cellType = cell->GetCellType();
  vtkIdList* cellIds = cell->GetPointIds();
  vtkCellArray* arrayvert = vtkCellArray::New();
  vtkPoints* cellPts = cell->GetPoints();
  vtkIdType npts = cellPts->GetNumberOfPoints();
  std::vector<vtkIdType> cellptId(npts);
  vtkIdType iid;
  const vtkIdType* v_id = nullptr;
  vtkIdType ptId;
  vtkIdType ptsvert = 1;

  int i;

  double v[3];

  for (i = 0; i < npts; i++)
  {
    cellptId[i] = cellIds->GetId(i);
  }

  // Convert all 0d cells to single vert.
  this->CellGrid(cellType, npts, cellptId.data(), arrayvert);

  unsigned int totalnewvert = arrayvert->GetNumberOfCells();
  for (unsigned int idlinenew = 0; idlinenew < totalnewvert; idlinenew++)
  {
    arrayvert->GetNextCell(ptsvert, v_id);

    // Clipping verts is easy.  Either it is inside the box or it isn't.
    cellPts->GetPoint(v_id[0], v);
    if ((v[0] >= this->BoundBoxClip[0][0]) && (v[0] <= this->BoundBoxClip[0][1]) &&
      (v[1] >= this->BoundBoxClip[1][0]) && (v[1] <= this->BoundBoxClip[1][1]) &&
      (v[2] >= this->BoundBoxClip[2][0]) && (v[2] <= this->BoundBoxClip[2][1]))
    {
      // Vert is inside.
      ptId = cellIds->GetId(v_id[0]);
      if (locator->InsertUniquePoint(v, iid))
      {
        outPD->CopyData(inPD, ptId, iid);
      }

      int newCellId = verts->InsertNextCell(1, &iid);
      outCD->CopyData(inCD, cellId, newCellId);
    }
  }
  arrayvert->Delete();
}

//------------------------------------------------------------------------------

void vtkBoxClipDataSet::ClipBoxInOut0D(vtkGenericCell* cell, vtkIncrementalPointLocator* locator,
  vtkCellArray** verts, vtkPointData* inPD, vtkPointData** outPD, vtkCellData* inCD,
  vtkIdType cellId, vtkCellData** outCD)
{
  vtkIdType cellType = cell->GetCellType();
  vtkIdList* cellIds = cell->GetPointIds();
  vtkCellArray* arrayvert = vtkCellArray::New();
  vtkPoints* cellPts = cell->GetPoints();
  vtkIdType npts = cellPts->GetNumberOfPoints();
  std::vector<vtkIdType> cellptId(npts);
  vtkIdType iid;
  const vtkIdType* v_id = nullptr;
  vtkIdType ptId;
  vtkIdType ptsvert = 1;

  int i;

  double v[3];

  for (i = 0; i < npts; i++)
  {
    cellptId[i] = cellIds->GetId(i);
  }

  // Convert all 0d cells to single vert.
  this->CellGrid(cellType, npts, cellptId.data(), arrayvert);

  unsigned int totalnewvert = arrayvert->GetNumberOfCells();
  for (unsigned int idlinenew = 0; idlinenew < totalnewvert; idlinenew++)
  {
    arrayvert->GetNextCell(ptsvert, v_id);

    // One way or another, we are adding the point.
    ptId = cellIds->GetId(v_id[0]);
    cellPts->GetPoint(v_id[0], v);

    if (locator->InsertUniquePoint(v, iid))
    {
      outPD[0]->CopyData(inPD, ptId, iid);
      outPD[1]->CopyData(inPD, ptId, iid);
    }

    // Clipping verts is easy.  Either it is inside the box or it isn't.
    if ((v[0] >= this->BoundBoxClip[0][0]) && (v[0] <= this->BoundBoxClip[0][1]) &&
      (v[1] >= this->BoundBoxClip[1][0]) && (v[1] <= this->BoundBoxClip[1][1]) &&
      (v[2] >= this->BoundBoxClip[2][0]) && (v[2] <= this->BoundBoxClip[2][1]))
    {
      // Vert is inside.
      int newCellId = verts[0]->InsertNextCell(1, &iid);
      outCD[0]->CopyData(inCD, cellId, newCellId);
    }
    else
    {
      // Vert is outside.
      int newCellId = verts[1]->InsertNextCell(1, &iid);
      outCD[1]->CopyData(inCD, cellId, newCellId);
    }
  }
  arrayvert->Delete();
}

//------------------------------------------------------------------------------

void vtkBoxClipDataSet::ClipHexahedron0D(vtkGenericCell* cell, vtkIncrementalPointLocator* locator,
  vtkCellArray* verts, vtkPointData* inPD, vtkPointData* outPD, vtkCellData* inCD, vtkIdType cellId,
  vtkCellData* outCD)
{
  vtkIdType cellType = cell->GetCellType();
  vtkIdList* cellIds = cell->GetPointIds();
  vtkCellArray* arrayvert = vtkCellArray::New();
  vtkPoints* cellPts = cell->GetPoints();
  vtkIdType npts = cellPts->GetNumberOfPoints();
  std::vector<vtkIdType> cellptId(npts);
  vtkIdType iid;
  const vtkIdType* v_id = nullptr;
  vtkIdType ptId;
  vtkIdType ptsvert = 1;

  int i;

  double v[3];

  for (i = 0; i < npts; i++)
  {
    cellptId[i] = cellIds->GetId(i);
  }

  // Convert all 0d cells to single vert.
  this->CellGrid(cellType, npts, cellptId.data(), arrayvert);

  unsigned int totalnewvert = arrayvert->GetNumberOfCells();
  for (unsigned int idlinenew = 0; idlinenew < totalnewvert; idlinenew++)
  {
    arrayvert->GetNextCell(ptsvert, v_id);

    // Clipping verts is easy.  Either it is inside the hexahedron or it isn't.
    cellPts->GetPoint(v_id[0], v);
    int inside = 1;

    for (int k = 0; k < 6; k++)
    {
      double value = (this->PlaneNormal[k][0] * (v[0] - this->PlanePoint[k][0]) +
        this->PlaneNormal[k][1] * (v[1] - this->PlanePoint[k][1]) +
        this->PlaneNormal[k][2] * (v[2] - this->PlanePoint[k][2]));
      if (value > 0)
      {
        inside = 0;
      }
    }

    if (inside)
    {
      // Vert is inside.
      ptId = cellIds->GetId(v_id[0]);
      if (locator->InsertUniquePoint(v, iid))
      {
        outPD->CopyData(inPD, ptId, iid);
      }

      int newCellId = verts->InsertNextCell(1, &iid);
      outCD->CopyData(inCD, cellId, newCellId);
    }
  }
  arrayvert->Delete();
}

//------------------------------------------------------------------------------

void vtkBoxClipDataSet::ClipHexahedronInOut0D(vtkGenericCell* cell,
  vtkIncrementalPointLocator* locator, vtkCellArray** verts, vtkPointData* inPD,
  vtkPointData** outPD, vtkCellData* inCD, vtkIdType cellId, vtkCellData** outCD)
{
  vtkIdType cellType = cell->GetCellType();
  vtkIdList* cellIds = cell->GetPointIds();
  vtkCellArray* arrayvert = vtkCellArray::New();
  vtkPoints* cellPts = cell->GetPoints();
  vtkIdType npts = cellPts->GetNumberOfPoints();
  std::vector<vtkIdType> cellptId(npts);
  vtkIdType iid;
  const vtkIdType* v_id = nullptr;
  vtkIdType ptId;
  vtkIdType ptsvert = 1;

  int i;

  double v[3];

  for (i = 0; i < npts; i++)
  {
    cellptId[i] = cellIds->GetId(i);
  }

  // Convert all 0d cells to single vert.
  this->CellGrid(cellType, npts, cellptId.data(), arrayvert);

  unsigned int totalnewvert = arrayvert->GetNumberOfCells();
  for (unsigned int idlinenew = 0; idlinenew < totalnewvert; idlinenew++)
  {
    arrayvert->GetNextCell(ptsvert, v_id);

    // One way or another, we are adding the point.
    ptId = cellIds->GetId(v_id[0]);
    cellPts->GetPoint(v_id[0], v);

    if (locator->InsertUniquePoint(v, iid))
    {
      outPD[0]->CopyData(inPD, ptId, iid);
      outPD[1]->CopyData(inPD, ptId, iid);
    }

    int inside = 1;
    for (int k = 0; k < 6; k++)
    {
      double value = (this->PlaneNormal[k][0] * (v[0] - this->PlanePoint[k][0]) +
        this->PlaneNormal[k][1] * (v[1] - this->PlanePoint[k][1]) +
        this->PlaneNormal[k][2] * (v[2] - this->PlanePoint[k][2]));
      if (value > 0)
      {
        inside = 0;
      }
    }

    // Clipping verts is easy.  Either it is inside the box or it isn't.
    if (inside)
    {
      // Vert is inside.
      int newCellId = verts[0]->InsertNextCell(1, &iid);
      outCD[0]->CopyData(inCD, cellId, newCellId);
    }
    else
    {
      // Vert is outside.
      int newCellId = verts[1]->InsertNextCell(1, &iid);
      outCD[1]->CopyData(inCD, cellId, newCellId);
    }
  }
  arrayvert->Delete();
}
VTK_ABI_NAMESPACE_END