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

#include "vtkCellArray.h"
#include "vtkCellArrayIterator.h"
#include "vtkCellData.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkMergePoints.h"
#include "vtkPointData.h"
#include "vtkPolyData.h"
#include "vtkSmartPointer.h"
#include "vtkTriangle.h"

VTK_ABI_NAMESPACE_BEGIN
vtkStandardNewMacro(vtkAdaptiveSubdivisionFilter);
vtkCxxSetObjectMacro(vtkAdaptiveSubdivisionFilter, Locator, vtkIncrementalPointLocator);

//------------------------------------------------------------------------------
// Construct object
vtkAdaptiveSubdivisionFilter::vtkAdaptiveSubdivisionFilter()
{
  this->MaximumEdgeLength = 1.0;
  this->MaximumTriangleArea = 1.0;
  this->MaximumNumberOfTriangles = VTK_ID_MAX;
  this->MaximumNumberOfPasses = VTK_ID_MAX;
  this->Locator = nullptr;
  this->OutputPointsPrecision = DEFAULT_PRECISION;
}

//------------------------------------------------------------------------------
// Construct object with number of subdivisions set to 1.
vtkAdaptiveSubdivisionFilter::~vtkAdaptiveSubdivisionFilter()
{
  this->SetLocator(nullptr);
}

//------------------------------------------------------------------------------
void vtkAdaptiveSubdivisionFilter::CreateDefaultLocator()
{
  if (this->Locator == nullptr)
  {
    this->Locator = vtkMergePoints::New();
    this->Locator->Register(this);
    this->Locator->Delete();
  }
}

//------------------------------------------------------------------------------
// Overload standard modified time function.
vtkMTimeType vtkAdaptiveSubdivisionFilter::GetMTime()
{
  vtkMTimeType mTime = this->Superclass::GetMTime();
  vtkMTimeType time;

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

  return mTime;
}

// Helper functions and structures
namespace
{

// There are eight possible subdivision cases (each of the three edges may
// or may not be subdivided). Case 0 just outputs the original triangle;
// the other cases output between 2 and four triangles. Note that when
// three triangles are generated, then the diagonal of the quadrilateral
// produced can go one of two ways. The tetCases is set up so that the two
// triangles forming the quad are the last two triangles and can be
// adjusted as necessary.
int CASE_MASK[3] = { 1, 2, 4 };
vtkIdType tessCases[16][13] = {
  { 1, 0, 1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, // case 0
  { 2, 0, 3, 2, 3, 1, 2, 0, 0, 0, 0, 0, 0 }, // case 1
  { 2, 0, 1, 4, 4, 2, 0, 0, 0, 0, 0, 0, 0 }, // case 2
  { 3, 3, 1, 4, 3, 4, 2, 2, 0, 3, 0, 0, 0 }, // case 3
  { 2, 0, 1, 5, 5, 1, 2, 0, 0, 0, 0, 0, 0 }, // case 4
  { 3, 0, 3, 5, 5, 3, 1, 1, 2, 5, 0, 0, 0 }, // case 5
  { 3, 5, 4, 2, 0, 1, 4, 4, 5, 0, 0, 0, 0 }, // case 6
  { 4, 0, 3, 5, 3, 1, 4, 5, 3, 4, 5, 4, 2 }, // case 7
  { 1, 0, 1, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, // case 0a
  { 2, 0, 3, 2, 3, 1, 2, 0, 0, 0, 0, 0, 0 }, // case 1a
  { 2, 0, 1, 4, 4, 2, 0, 0, 0, 0, 0, 0, 0 }, // case 2a
  { 3, 3, 1, 4, 0, 3, 4, 4, 2, 0, 0, 0, 0 }, // case 3a
  { 2, 0, 1, 5, 5, 1, 2, 0, 0, 0, 0, 0, 0 }, // case 4a
  { 3, 0, 3, 5, 3, 1, 2, 2, 5, 3, 0, 0, 0 }, // case 5a
  { 3, 4, 2, 5, 5, 0, 1, 1, 4, 5, 0, 0, 0 }, // case 6a
  { 4, 0, 3, 5, 3, 1, 4, 5, 3, 4, 5, 4, 2 }, // case 7a
};

// This method assumes that the diagonal of the quadrilateral formed by
// triangles 2 & 3 may be "swapped" to produce a better triangulation. It
// assumes a lot about the ordering of the connectivity array (subTess).
vtkIdType* SelectTessellation(unsigned char subCase, vtkIdType* ptIds, vtkPoints* newPts)
{
  // If no choice in triangulation return the table entry
  if (tessCases[subCase][0] != 3)
  {
    return tessCases[subCase];
  }

  // Else select best triangulation based on diagonal length
  vtkIdType* subTess = tessCases[subCase];
  double x0[3], x1[3], x2[3], x3[3];
  newPts->GetPoint(ptIds[subTess[4]], x0);
  newPts->GetPoint(ptIds[subTess[6]], x1);
  newPts->GetPoint(ptIds[subTess[5]], x2);
  newPts->GetPoint(ptIds[subTess[8]], x3);

  if (vtkMath::Distance2BetweenPoints(x0, x1) <= vtkMath::Distance2BetweenPoints(x2, x3))
  {
    return tessCases[subCase];
  }
  else
  {
    return tessCases[subCase + 8]; // alternate triangulation
  }
}

} // anonymous namespace

//------------------------------------------------------------------------------
// This uses a very simple, serial implementation that makes repeated passes
// over the triangles using a swap buffer approach.
int vtkAdaptiveSubdivisionFilter::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 and check its validity
  vtkPolyData* input = vtkPolyData::SafeDownCast(inInfo->Get(vtkDataObject::DATA_OBJECT()));
  vtkPolyData* output = vtkPolyData::SafeDownCast(outInfo->Get(vtkDataObject::DATA_OBJECT()));

  vtkIdType numPts = input->GetNumberOfPoints();
  vtkCellArray* inTris = input->GetPolys();
  vtkIdType numTris = inTris->GetNumberOfCells();
  if (numPts < 1 || numTris < 1)
  {
    vtkDebugMacro(<< "No data to subdivide!");
    return 1;
  }
  vtkPointData* inPointData = input->GetPointData();
  vtkCellData* inCellData = input->GetCellData();

  if (inTris->IsHomogeneous() != 3)
  {
    vtkDebugMacro(<< "Filter operates only on triangles!");
    return 1;
  }

  // Need a locator
  if (!this->Locator)
  {
    this->CreateDefaultLocator();
  }

  // The first thing is to take the existing points and push them into the
  // incremental point locator. We know that we are going to use the original
  // points. Note that points are only created and are not swapped as each
  // pass is invoked.
  vtkPoints* inPts = input->GetPoints();
  vtkPoints* newPts = vtkPoints::New();
  vtkPointData *swapPointData, *newPointData = vtkPointData::New();
  newPointData->CopyAllocate(inPointData);

  // set precision for the points in the output
  if (this->OutputPointsPrecision == vtkAlgorithm::DEFAULT_PRECISION)
  {
    newPts->SetDataType(inPts->GetDataType());
  }
  else if (this->OutputPointsPrecision == vtkAlgorithm::SINGLE_PRECISION)
  {
    newPts->SetDataType(VTK_FLOAT);
  }
  else if (this->OutputPointsPrecision == vtkAlgorithm::DOUBLE_PRECISION)
  {
    newPts->SetDataType(VTK_DOUBLE);
  }
  this->Locator->InitPointInsertion(newPts, input->GetBounds(), input->GetNumberOfPoints());
  // Load in the already existing points. Also load in the point data
  // associated with the existing points.
  for (vtkIdType ptId = 0; ptId < numPts; ++ptId)
  {
    this->Locator->InsertNextPoint(inPts->GetPoint(ptId));
    newPointData->CopyData(inPointData, ptId, ptId);
  }

  // This is a multipass algorithm. From a list of triangles, check each
  // against the edge length and area criteria. If necessary, break the
  // triangle (using a case table) into smaller triangles by inserting one or
  // more points on edges (the edge is broken at its midpoint). The new
  // triangles are placed into a new list which serves as the starting point
  // for the next pass. An important note: triangles are split independently
  // without neighbor "links" (i.e.,cell links) and new points are merged
  // into the locator. Since the algorithm treats edges on triangles in an
  // identical way, the end result is that triangle neighbors remain
  // compatible (due to coincident point merging).
  auto cellIter = vtk::TakeSmartPointer(inTris->NewIterator());
  vtkCellArray *swapTris, *newTris = vtkCellArray::New();
  newTris->AllocateEstimate(2 * numTris, 3);
  vtkCellData *swapCellData, *newCellData = vtkCellData::New();
  newCellData->CopyAllocate(inCellData);

  int i;
  double area, eLengths[3];
  double maxLen2 = this->MaximumEdgeLength * this->MaximumEdgeLength;
  double maxArea = this->MaximumTriangleArea;
  double x[6][3]; // three vertices plus potential mid-edge points
  const vtkIdType* tri;
  vtkIdType triId;
  vtkIdType newId;
  vtkIdType passNum;
  vtkIdType totalTriangles = 0;
  bool changesMade;
  bool abort = false;

  for (passNum = 0, changesMade = true; passNum < this->MaximumNumberOfPasses &&
       totalTriangles < this->MaximumNumberOfTriangles && changesMade && !abort;
       ++passNum)
  {
    changesMade = false;
    for (cellIter->GoToFirstCell(); !cellIter->IsDoneWithTraversal(); cellIter->GoToNextCell())
    {
      if (this->CheckAbort())
      {
        abort = true;
        break;
      }
      triId = cellIter->GetCurrentCellId();
      {
        vtkIdType unused;
        cellIter->GetCurrentCell(unused, tri);
      }

      newPts->GetPoint(tri[0], x[0]);
      newPts->GetPoint(tri[1], x[1]);
      newPts->GetPoint(tri[2], x[2]);
      eLengths[0] = vtkMath::Distance2BetweenPoints(x[0], x[1]);
      eLengths[1] = vtkMath::Distance2BetweenPoints(x[1], x[2]);
      eLengths[2] = vtkMath::Distance2BetweenPoints(x[2], x[0]);
      area = vtkTriangle::TriangleArea(x[0], x[1], x[2]);

      // Various subdivision cases are possible
      unsigned char subCase = 0;
      if (area > maxArea)
      {
        subCase = 7;
      }
      else
      {
        for (i = 0; i < 3; ++i)
        {
          if (eLengths[i] > maxLen2)
          {
            subCase |= CASE_MASK[i];
          }
        }
      } // determine edges to divide

      // If not just outputting original triangle then changes are made
      if (subCase > 0)
      {
        changesMade = true;
      }

      // Now create new points and triangles dividing edges as appropriate.
      double xNew[3];
      vtkIdType ptIds[6];
      ptIds[0] = tri[0];
      ptIds[1] = tri[1];
      ptIds[2] = tri[2];
      for (i = 0; i < 3; ++i)
      {
        if (subCase & CASE_MASK[i]) // ith edge needs subdivision
        {
          xNew[0] = 0.5 * (x[i][0] + x[(i + 1) % 3][0]);
          xNew[1] = 0.5 * (x[i][1] + x[(i + 1) % 3][1]);
          xNew[2] = 0.5 * (x[i][2] + x[(i + 1) % 3][2]);
          if ((ptIds[3 + i] = this->Locator->IsInsertedPoint(xNew)) < 0)
          {
            ptIds[3 + i] = this->Locator->InsertNextPoint(xNew);
            newPointData->InterpolateEdge(inPointData, ptIds[3 + i], tri[i], tri[(i + 1) % 3], 0.5);
          }
        }
      }

      // The tessellation may vary based on geometric concerns (selecting best
      // diagonal during triangulation of quadrilateral)
      vtkIdType newTIds[3], *subTess;
      subTess = SelectTessellation(subCase, ptIds, newPts);
      vtkIdType numTessTris = *subTess++;

      for (i = 0; i < numTessTris; ++i, subTess += 3)
      {
        newTIds[0] = ptIds[subTess[0]];
        newTIds[1] = ptIds[subTess[1]];
        newTIds[2] = ptIds[subTess[2]];
        newId = newTris->InsertNextCell(3, newTIds);
        newCellData->CopyData(inCellData, triId, newId);
        if (++totalTriangles >= this->MaximumNumberOfTriangles)
        {
          break;
        }
      }
    } // for all triangles in this pass

    // Prepare for the next pass, which means swapping input and output.
    // Remember that the initial pass uses the filter input; subsequent passes
    // cannot modify the input to a new cell array must be created to support
    // the swapping.
    if (passNum == 0)
    {
      inTris = vtkCellArray::New();
      inCellData = vtkCellData::New();
      inCellData->CopyAllocate(newCellData);

      inPointData = vtkPointData::New();
      inPointData->CopyAllocate(newPointData);
    }

    // Prepare for new triangles
    swapTris = newTris;
    newTris = inTris;
    inTris = swapTris;
    cellIter = vtk::TakeSmartPointer(inTris->NewIterator());

    numTris = inTris->GetNumberOfCells();
    newTris->Reset();
    newTris->AllocateEstimate(2 * numTris, 3);

    // Prepare for new cell data
    swapCellData = newCellData;
    newCellData = inCellData;
    inCellData = swapCellData;

    // Prepare for new point data
    numPts = newPts->GetNumberOfPoints();
    swapPointData = newPointData;
    newPointData = inPointData;
    inPointData = swapPointData;
    for (vtkIdType ptId = 0; ptId < numPts; ++ptId)
    {
      newPointData->CopyData(inPointData, ptId, ptId);
    }
  } // for another pass

  // Configure output and clean up
  output->SetPoints(newPts);
  newPts->Delete();
  output->GetPointData()->ShallowCopy(inPointData);
  newPointData->Delete();

  output->SetPolys(inTris);
  inTris->Delete();
  newTris->Delete();
  output->GetCellData()->ShallowCopy(inCellData);
  inCellData->Delete();
  inPointData->Delete();
  newCellData->Delete();

  return 1;
}

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

  os << indent << "Maximum Edge Length: " << this->MaximumEdgeLength << "\n";
  os << indent << "Maximum Triangle Area: " << this->MaximumTriangleArea << "\n";
  os << indent << "Maximum Number Of Triangles: " << this->MaximumNumberOfTriangles << "\n";
  os << indent << "Maximum Number Of Passes: " << this->MaximumNumberOfPasses << "\n";

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

  os << indent << "Precision of the output points: " << this->OutputPointsPrecision << "\n";
}
VTK_ABI_NAMESPACE_END