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

  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 "vtkMath.h"
#include "vtkObjectFactory.h"
#include "vtkRenderer.h"
#include "vtkIntArray.h"

vtkStandardNewMacro(vtkBezierContourLineInterpolator);

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

//----------------------------------------------------------------------
vtkBezierContourLineInterpolator::~vtkBezierContourLineInterpolator()
{
}

//----------------------------------------------------------------------
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->InsertNextTupleValue( 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";
}