// SPDX-FileCopyrightText: Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
// SPDX-License-Identifier: BSD-3-Clause

#include "vtkDiscreteMarchingCubes.h"

#include "vtkArrayDispatch.h"
#include "vtkCellArray.h"
#include "vtkCellData.h"
#include "vtkDataArrayRange.h"
#include "vtkDoubleArray.h"
#include "vtkFloatArray.h"
#include "vtkImageTransform.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkMarchingCubesTriangleCases.h"
#include "vtkMath.h"
#include "vtkMergePoints.h"
#include "vtkObjectFactory.h"
#include "vtkPointData.h"
#include "vtkPolyData.h"
#include "vtkStreamingDemandDrivenPipeline.h"
#include "vtkStructuredPoints.h"

VTK_ABI_NAMESPACE_BEGIN
vtkStandardNewMacro(vtkDiscreteMarchingCubes);

void vtkDiscreteMarchingCubes::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os, indent);
  os << indent << "ComputeAdjacentScalars: " << this->ComputeAdjacentScalars << endl;
}

// Description:
// Construct object with initial range (0,1) and single contour value
// of 0.0. ComputeNormals is off, ComputeGradients is off, ComputeScalars is on and
// ComputeAdjacentScalars is off.
vtkDiscreteMarchingCubes::vtkDiscreteMarchingCubes()
{
  this->ComputeNormals = 0;
  this->ComputeGradients = 0;
  this->ComputeScalars = 1;
  this->ComputeAdjacentScalars = 0;
}

vtkDiscreteMarchingCubes::~vtkDiscreteMarchingCubes() = default;

//
// Contouring filter specialized for volumes and "short int" data values.
//
template <int TupleSize>
struct vtkDiscreteMarchingCubesComputeGradientFunctor
{
  template <class TArray>
  void operator()(TArray* scalarArray, vtkDiscreteMarchingCubes* self, int dims[3],
    vtkIncrementalPointLocator* locator, vtkDataArray* newCellScalars,
    vtkDataArray* newPointScalars, vtkCellArray* newPolys, double* values, int numValues)
  {
    double s[8], value;
    int i, j, k, pts[8][3], xp, yp, zp, *x1, *x2;
    vtkIdType sliceSize, rowSize;
    static const int CASE_MASK[8] = { 1, 2, 4, 8, 16, 32, 64, 128 };
    vtkMarchingCubesTriangleCases *triCase, *triCases;
    int* edge;
    int contNum, ii, index, *vert;
    vtkIdType jOffset, kOffset, idx;
    vtkIdType ptIds[3];
    int extent[6];
    vtkTypeBool ComputeScalars = newCellScalars != nullptr;
    int ComputeAdjacentScalars = newPointScalars != nullptr;
    double t, x[3], min, max;
    static int edges[12][2] = { { 0, 1 }, { 1, 2 }, { 3, 2 }, { 0, 3 }, { 4, 5 }, { 5, 6 },
      { 7, 6 }, { 4, 7 }, { 0, 4 }, { 1, 5 }, { 3, 7 }, { 2, 6 } };

    vtkInformation* inInfo = self->GetExecutive()->GetInputInformation(0, 0);
    inInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), extent);

    triCases = vtkMarchingCubesTriangleCases::GetCases();
    auto scalars = vtk::DataArrayValueRange<TupleSize>(scalarArray).begin();

    //
    // Get min/max contour values
    //
    if (numValues < 1)
    {
      return;
    }
    for (min = max = values[0], i = 1; i < numValues; i++)
    {
      min = std::min(values[i], min);
      max = std::max(values[i], max);
    }
    //
    // Traverse all voxel cells, generating triangles
    // using marching cubes algorithm.
    //
    rowSize = dims[0];
    sliceSize = rowSize * dims[1];
    for (k = 0; k < (dims[2] - 1); k++)
    {
      self->UpdateProgress(static_cast<double>(k) / (dims[2] - 1));
      if (self->CheckAbort())
      {
        break;
      }
      kOffset = k * sliceSize;
      pts[0][2] = k + extent[4];
      zp = pts[0][2] + 1;
      for (j = 0; j < (dims[1] - 1); j++)
      {
        jOffset = j * rowSize;
        pts[0][1] = j + extent[2];
        yp = pts[0][1] + 1;
        for (i = 0; i < (dims[0] - 1); i++)
        {
          // get scalar values
          idx = i + jOffset + kOffset;
          s[0] = scalars[idx];
          s[1] = scalars[idx + 1];
          s[2] = scalars[idx + 1 + dims[0]];
          s[3] = scalars[idx + dims[0]];
          s[4] = scalars[idx + sliceSize];
          s[5] = scalars[idx + 1 + sliceSize];
          s[6] = scalars[idx + 1 + dims[0] + sliceSize];
          s[7] = scalars[idx + dims[0] + sliceSize];

          if ((s[0] < min && s[1] < min && s[2] < min && s[3] < min && s[4] < min && s[5] < min &&
                s[6] < min && s[7] < min) ||
            (s[0] > max && s[1] > max && s[2] > max && s[3] > max && s[4] > max && s[5] > max &&
              s[6] > max && s[7] > max))
          {
            continue; // no contours possible
          }

          // create voxel points
          pts[0][0] = i + extent[0];
          xp = pts[0][0] + 1;

          pts[1][0] = xp;
          pts[1][1] = pts[0][1];
          pts[1][2] = pts[0][2];

          pts[2][0] = xp;
          pts[2][1] = yp;
          pts[2][2] = pts[0][2];

          pts[3][0] = pts[0][0];
          pts[3][1] = yp;
          pts[3][2] = pts[0][2];

          pts[4][0] = pts[0][0];
          pts[4][1] = pts[0][1];
          pts[4][2] = zp;

          pts[5][0] = xp;
          pts[5][1] = pts[0][1];
          pts[5][2] = zp;

          pts[6][0] = xp;
          pts[6][1] = yp;
          pts[6][2] = zp;

          pts[7][0] = pts[0][0];
          pts[7][1] = yp;
          pts[7][2] = zp;

          for (contNum = 0; contNum < numValues; contNum++)
          {
            value = values[contNum];
            // Build the case table
            for (ii = 0, index = 0; ii < 8; ii++)
            {
              // for discrete marching cubes, we are looking for an
              // exact match of a scalar at a vertex to a value
              if (s[ii] == value)
              {
                index |= CASE_MASK[ii];
              }
            }
            if (index == 0 || index == 255) // no surface
            {
              continue;
            }

            triCase = triCases + index;
            edge = triCase->edges;

            for (; edge[0] > -1; edge += 3)
            {
              for (ii = 0; ii < 3; ii++) // insert triangle
              {
                vert = edges[edge[ii]];
                // for discrete marching cubes, the interpolation point
                // is always 0.5.
                t = 0.5;
                x1 = pts[vert[0]];
                x2 = pts[vert[1]];
                x[0] = x1[0] + t * (x2[0] - x1[0]);
                x[1] = x1[1] + t * (x2[1] - x1[1]);
                x[2] = x1[2] + t * (x2[2] - x1[2]);

                // add point
                if (locator->InsertUniquePoint(x, ptIds[ii]))
                {
                  if (ComputeAdjacentScalars)
                  {
                    // check which vert holds the neighbour value
                    if (s[vert[0]] == value)
                    {
                      newPointScalars->InsertTuple(ptIds[ii], &s[vert[1]]);
                    }
                    else
                    {
                      newPointScalars->InsertTuple(ptIds[ii], &s[vert[0]]);
                    }
                  }
                }
              }
              // check for degenerate triangle
              if (ptIds[0] != ptIds[1] && ptIds[0] != ptIds[2] && ptIds[1] != ptIds[2])
              {
                newPolys->InsertNextCell(3, ptIds);
                // Note that DiscreteMarchingCubes stores the scalar
                // data in the cells. It does not use the point data
                // since cells from different labeled segments may use
                // the same point.
                if (ComputeScalars)
                {
                  newCellScalars->InsertNextTuple(&value);
                }
              }
            } // for each triangle
          }   // for all contours
        }     // for i
      }       // for j
    }         // for k
  }
};

//
// Contouring filter specialized for volumes and "short int" data values.
//
int vtkDiscreteMarchingCubes::RequestData(vtkInformation* vtkNotUsed(request),
  vtkInformationVector** inputVector, vtkInformationVector* outputVector)
{
  vtkInformation* inInfo = inputVector[0]->GetInformationObject(0);
  vtkInformation* outInfo = outputVector->GetInformationObject(0);

  vtkPoints* newPts;
  vtkCellArray* newPolys;
  vtkFloatArray* newCellScalars;
  vtkFloatArray* newPointScalars;
  vtkImageData* input = vtkImageData::SafeDownCast(inInfo->Get(vtkDataObject::DATA_OBJECT()));
  vtkPointData* pd;
  vtkDataArray* inScalars;
  int dims[3], extent[6];
  vtkIdType estimatedSize;
  double bounds[6];
  vtkPolyData* output = vtkPolyData::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT()));
  vtkIdType numContours = this->ContourValues->GetNumberOfContours();
  double* values = this->ContourValues->GetValues();

  vtkDebugMacro(<< "Executing marching cubes");

  // initialize and check input
  pd = input->GetPointData();
  if (pd == nullptr)
  {
    vtkErrorMacro(<< "PointData is nullptr");
    return 1;
  }
  vtkInformationVector* inArrayVec = this->Information->Get(INPUT_ARRAYS_TO_PROCESS());
  if (inArrayVec)
  { // we have been passed an input array
    inScalars = this->GetInputArrayToProcess(0, inputVector);
  }
  else
  {
    inScalars = pd->GetScalars();
  }
  if (inScalars == nullptr)
  {
    vtkErrorMacro(<< "Scalars must be defined for contouring");
    return 1;
  }

  if (input->GetDataDimension() != 3)
  {
    vtkErrorMacro(<< "Cannot contour data of dimension != 3");
    return 1;
  }
  input->GetDimensions(dims);

  inInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), extent);

  // estimate the number of points from the volume dimensions
  estimatedSize = dims[0];
  estimatedSize *= dims[1]; // The "*=" ensures coercion to vtkIdType,
  estimatedSize *= dims[2]; // which might be wider than "int"
  estimatedSize = static_cast<vtkIdType>(pow(static_cast<double>(estimatedSize), .75));
  estimatedSize = estimatedSize / 1024 * 1024; // multiple of 1024
  estimatedSize = std::max<vtkIdType>(estimatedSize, 1024);
  vtkDebugMacro(<< "Estimated allocation size is " << estimatedSize);
  newPts = vtkPoints::New();
  newPts->Allocate(estimatedSize, estimatedSize / 2);

  // compute bounds for merging points
  for (int i = 0; i < 3; i++)
  {
    bounds[2 * i] = extent[2 * i];
    bounds[2 * i + 1] = extent[2 * i + 1];
  }
  if (this->Locator == nullptr)
  {
    this->CreateDefaultLocator();
  }
  this->Locator->InitPointInsertion(newPts, bounds, estimatedSize);

  newPolys = vtkCellArray::New();
  newPolys->AllocateEstimate(estimatedSize, 3);

  if (this->ComputeScalars)
  {
    newCellScalars = vtkFloatArray::New();
    newCellScalars->SetName("Scalars");
    newCellScalars->Allocate(estimatedSize, 3);
  }
  else
  {
    newCellScalars = nullptr;
  }

  if (this->ComputeAdjacentScalars)
  {
    newPointScalars = vtkFloatArray::New();
    newPointScalars->SetName("AdjacentScalars");
    newPointScalars->Allocate(estimatedSize, estimatedSize / 2);
  }
  else
  {
    newPointScalars = nullptr;
  }

  if (inScalars->GetNumberOfComponents() == 1)
  {
    vtkDiscreteMarchingCubesComputeGradientFunctor<1> functor;
    if (!vtkArrayDispatch::Dispatch::Execute(inScalars, functor, this, dims, this->Locator,
          newCellScalars, newPointScalars, newPolys, values, numContours))
    {
      functor(inScalars, this, dims, this->Locator, newCellScalars, newPointScalars, newPolys,
        values, numContours);
    }
  }

  else // multiple components - have to convert
  {
    vtkIdType dataSize = dims[0];
    dataSize *= dims[1]; // The "*=" ensures coercion to vtkIdType,
    dataSize *= dims[2]; // which might be wider than "int".
    vtkNew<vtkDoubleArray> image;
    image->SetNumberOfComponents(inScalars->GetNumberOfComponents());
    image->SetNumberOfTuples(image->GetNumberOfComponents() * dataSize);
    inScalars->GetTuples(0, dataSize, image);

    vtkDiscreteMarchingCubesComputeGradientFunctor<vtk::detail::DynamicTupleSize> functor;
    functor(image.Get(), this, dims, this->Locator, newCellScalars, newPointScalars, newPolys,
      values, numContours);
  }

  vtkDebugMacro(<< "Created: " << newPts->GetNumberOfPoints() << " points, "
                << newPolys->GetNumberOfCells() << " triangles");
  //
  // Update ourselves.  Because we don't know up front how many triangles
  // we've created, take care to reclaim memory.
  //
  output->SetPoints(newPts);
  newPts->Delete();

  output->SetPolys(newPolys);
  newPolys->Delete();

  if (newCellScalars)
  {
    output->GetCellData()->SetScalars(newCellScalars);
    newCellScalars->Delete();
  }
  if (newPointScalars)
  {
    int idx = output->GetPointData()->AddArray(newPointScalars);
    output->GetPointData()->SetActiveAttribute(idx, vtkDataSetAttributes::SCALARS);
    newPointScalars->Delete();
  }
  output->Squeeze();
  if (this->Locator)
  {
    this->Locator->Initialize(); // free storage
  }

  vtkImageTransform::TransformPointSet(input, output);

  return 1;
}
VTK_ABI_NAMESPACE_END