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

  Program:   Visualization Toolkit
  Module:    vtkBezierContourLineInterpolator.cxx

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

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

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

#include "vtkContourRepresentation.h"
#include "vtkIntArray.h"
#include "vtkMath.h"
#include "vtkObjectFactory.h"
#include "vtkRenderer.h"

vtkStandardNewMacro(vtkBezierContourLineInterpolator);

//------------------------------------------------------------------------------
vtkBezierContourLineInterpolator::vtkBezierContourLineInterpolator()
{
  this->MaximumCurveError = 0.005;
  this->MaximumCurveLineSegments = 100;
}

//------------------------------------------------------------------------------
vtkBezierContourLineInterpolator::~vtkBezierContourLineInterpolator() = default;

//------------------------------------------------------------------------------
int vtkBezierContourLineInterpolator::InterpolateLine(
  vtkRenderer* vtkNotUsed(ren), vtkContourRepresentation* rep, int idx1, int idx2)
{
  int maxRecursion = 0;
  int tmp = 3;

  while (2 * tmp < this->MaximumCurveLineSegments)
  {
    tmp *= 2;
    maxRecursion++;
  }

  if (maxRecursion == 0)
  {
    return 1;
  }

  // There are four control points with 3 components each, plus one
  // value for the recursion depth of this point
  double* controlPointsStack = new double[(3 * 4 + 1) * (maxRecursion + 1)];
  int stackCount = 0;

  double slope1[3];
  double slope2[3];

  rep->GetNthNodeSlope(idx1, slope1);
  rep->GetNthNodeSlope(idx2, slope2);

  controlPointsStack[0] = 0;
  double* p1 = controlPointsStack + 1;
  double* p2 = controlPointsStack + 4;
  double* p3 = controlPointsStack + 7;
  double* p4 = controlPointsStack + 10;

  rep->GetNthNodeWorldPosition(idx1, p1);
  rep->GetNthNodeWorldPosition(idx2, p4);

  double distance = sqrt(vtkMath::Distance2BetweenPoints(p1, p4));

  p2[0] = p1[0] + .333 * distance * slope1[0];
  p2[1] = p1[1] + .333 * distance * slope1[1];
  p2[2] = p1[2] + .333 * distance * slope1[2];

  p3[0] = p4[0] - .333 * distance * slope2[0];
  p3[1] = p4[1] - .333 * distance * slope2[1];
  p3[2] = p4[2] - .333 * distance * slope2[2];

  stackCount++;

  while (stackCount)
  {
    // process last point on stack
    int recursionLevel = static_cast<int>(controlPointsStack[13 * (stackCount - 1)]);

    p1 = controlPointsStack + 13 * (stackCount - 1) + 1;
    p2 = controlPointsStack + 13 * (stackCount - 1) + 4;
    p3 = controlPointsStack + 13 * (stackCount - 1) + 7;
    p4 = controlPointsStack + 13 * (stackCount - 1) + 10;

    double totalDist = 0;
    totalDist += sqrt(vtkMath::Distance2BetweenPoints(p1, p2));
    totalDist += sqrt(vtkMath::Distance2BetweenPoints(p2, p3));
    totalDist += sqrt(vtkMath::Distance2BetweenPoints(p3, p4));

    distance = sqrt(vtkMath::Distance2BetweenPoints(p1, p4));

    if (recursionLevel >= maxRecursion || distance == 0 ||
      (totalDist - distance) / distance < this->MaximumCurveError)
    {
      rep->AddIntermediatePointWorldPosition(idx1, p2);
      rep->AddIntermediatePointWorldPosition(idx1, p3);

      if (stackCount > 1)
      {
        rep->AddIntermediatePointWorldPosition(idx1, p4);
      }
      stackCount--;
    }
    else
    {
      double p12[3], p23[3], p34[3], p123[3], p234[3], p1234[3];

      this->ComputeMidpoint(p1, p2, p12);
      this->ComputeMidpoint(p2, p3, p23);
      this->ComputeMidpoint(p3, p4, p34);
      this->ComputeMidpoint(p12, p23, p123);
      this->ComputeMidpoint(p23, p34, p234);
      this->ComputeMidpoint(p123, p234, p1234);

      // add these two points to the stack
      controlPointsStack[13 * (stackCount - 1)] = recursionLevel + 1;
      controlPointsStack[13 * (stackCount)] = recursionLevel + 1;

      double* newp1 = controlPointsStack + 13 * (stackCount) + 1;
      double* newp2 = controlPointsStack + 13 * (stackCount) + 4;
      double* newp3 = controlPointsStack + 13 * (stackCount) + 7;
      double* newp4 = controlPointsStack + 13 * (stackCount) + 10;

      newp1[0] = p1[0];
      newp1[1] = p1[1];
      newp1[2] = p1[2];

      newp2[0] = p12[0];
      newp2[1] = p12[1];
      newp2[2] = p12[2];

      newp3[0] = p123[0];
      newp3[1] = p123[1];
      newp3[2] = p123[2];

      newp4[0] = p1234[0];
      newp4[1] = p1234[1];
      newp4[2] = p1234[2];

      p1[0] = p1234[0];
      p1[1] = p1234[1];
      p1[2] = p1234[2];

      p2[0] = p234[0];
      p2[1] = p234[1];
      p2[2] = p234[2];

      p3[0] = p34[0];
      p3[1] = p34[1];
      p3[2] = p34[2];

      stackCount++;
    }
  }

  delete[] controlPointsStack;

  return 1;
}

//------------------------------------------------------------------------------
void vtkBezierContourLineInterpolator::GetSpan(
  int nodeIndex, vtkIntArray* nodeIndices, vtkContourRepresentation* rep)
{
  int start = nodeIndex - 2;
  int end = nodeIndex - 1;
  int index[2];

  // Clear the array
  nodeIndices->Reset();
  nodeIndices->Squeeze();
  nodeIndices->SetNumberOfComponents(2);

  for (int i = 0; i < 4; i++)
  {
    index[0] = start++;
    index[1] = end++;

    if (rep->GetClosedLoop())
    {
      if (index[0] < 0)
      {
        index[0] += rep->GetNumberOfNodes();
      }
      if (index[1] < 0)
      {
        index[1] += rep->GetNumberOfNodes();
      }
      if (index[0] >= rep->GetNumberOfNodes())
      {
        index[0] -= rep->GetNumberOfNodes();
      }
      if (index[1] >= rep->GetNumberOfNodes())
      {
        index[1] -= rep->GetNumberOfNodes();
      }
    }

    if (index[0] >= 0 && index[0] < rep->GetNumberOfNodes() && index[1] >= 0 &&
      index[1] < rep->GetNumberOfNodes())
    {
      nodeIndices->InsertNextTypedTuple(index);
    }
  }
}

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

  os << indent << "Maximum Curve Error: " << this->MaximumCurveError << "\n";
  os << indent << "Maximum Curve Line Segments: " << this->MaximumCurveLineSegments << "\n";
}