// SOURCE: CATMULL-ROM CUBIC SPLINE
// AUTHOR: SEAN M. KESSLER (Adopted from Numerical Recipes in 'C')
// DATE: 12/28/1993
#include <analytic/catmull.hpp>

CatmullRom::CatmullRom(void)
{
}

CatmullRom::CatmullRom(const CatmullRom &someCatmullRom)
{ // no implementation
  *this=someCatmullRom;
}

inline
CatmullRom::~CatmullRom()
{
}

inline
CatmullRom &CatmullRom::operator=(const CatmullRom &/*someCatmullRom*/)
{ // no implementation
  return *this;
}

BOOL CatmullRom::performSpline(GlobalData<FloatPairs> &sourcePairs,GlobalData<FloatPairs> &destPairs)
{
  double a0,a1,a2,a3;
  double dx,dx1,dx2;
  double dy,dy1,dy2;
  double endPointOne,endPointTwo,resamplingPos;
  double xPoint;
  short clampOne,clampTwo;
  short direction;
  short destSize((short)destPairs.size());
  short sourceSize((short)sourcePairs.size());
  short inputIndex,index;

  if(sourceSize<2||destSize<2)return FALSE;
  if(sourcePairs[0].column()<sourcePairs[1].column())
  {
    if(destPairs[0].column()<sourcePairs[0].column()||
       destPairs[destSize-1].column()>sourcePairs[sourceSize-1].column())direction=0;
      else direction=1;
  }
  else
  {
    if(destPairs[0].column()>sourcePairs[0].column()||
       destPairs[destSize-1].column()<sourcePairs[sourceSize-1].column())direction=0;
      else direction=-1;
  }
  if(!direction)return FALSE;
  if(1==direction)endPointTwo=destPairs[0].column()-1;
  else endPointTwo=destPairs[0].column()+1;
  for(index=0;index<destSize;index++)
  {
    resamplingPos=destPairs[index].column();
    if((1==direction&&resamplingPos>endPointTwo)||(-1==direction&&resamplingPos<endPointTwo))
    {
      for(inputIndex=0;inputIndex<sourceSize&&resamplingPos>sourcePairs[inputIndex].column();inputIndex++);
      if(resamplingPos<sourcePairs[inputIndex].column())inputIndex--;
      if(inputIndex<0)inputIndex=0;
      endPointOne=sourcePairs[inputIndex].column();
      endPointTwo=sourcePairs[inputIndex+1].column();
      clampOne=FloatPairs::fpmax(inputIndex-1,0);
      clampTwo=FloatPairs::fpmin(inputIndex+2,sourceSize-1);
      dx=1.0/(endPointTwo-endPointOne);
      dx1=1.0/(endPointTwo-sourcePairs[clampOne].column());
      dx2=1.0/(sourcePairs[clampTwo].column()-endPointOne);
      dy=(sourcePairs[inputIndex+1].row()-sourcePairs[inputIndex].row())*dx;
      dy1=(sourcePairs[inputIndex+1].row()-sourcePairs[clampOne].row())*dx1;
      dy2=(sourcePairs[clampTwo].row()-sourcePairs[inputIndex].row())*dx2;
      a0=sourcePairs[inputIndex].row();
      a1=dy1;
      a2=dx*(3*dy-2*dy1-dy2);
      a3=dx*dx*(-2*dy+dy1+dy2);
    }
    xPoint=resamplingPos-endPointOne;
    destPairs[index].row(((a3*xPoint+a2)*xPoint+a1)*xPoint+a0);
  }
  return TRUE;
}