//#define TESTBSPLINEFLOOR

#include "itkBSplineDeformableTransform.h"
#include "itkNumericTraits.h"
#include <iostream>

int main( int argc, char ** argv )
{
  const unsigned int Dimension = 3;
  typedef double CoordRepType;

  typedef itk::BSplineDeformableTransform<CoordRepType, Dimension, 3>
    TransformType;
  typedef TransformType::ParametersType ParametersType;
  typedef TransformType::SizeType SizeType;
  typedef TransformType::RegionType RegionType;
  typedef TransformType::OriginType OriginType;
  typedef TransformType::SpacingType SpacingType;
  typedef TransformType::InputPointType PointType;
  typedef PointType::VectorType VectorType;

  TransformType::Pointer transform = TransformType::New();
  SizeType size;
  OriginType origin;
  SpacingType spacing;
  RegionType region;
  ParametersType parameters;
  PointType inpoint;
  PointType outpoint;
  VectorType vec;

  for (unsigned int i =0; i < Dimension ; ++i)
  {
    size[i]=100;
    origin[i]=0.0;
    spacing[i]=1.0;
  }
  region.SetSize( size);
  transform->SetGridRegion( region );
  transform->SetGridOrigin( origin );
  transform->SetGridSpacing( spacing );
  unsigned int N = transform->GetNumberOfParameters();
  parameters.SetSize(N);
  parameters.Fill(0.0);
  for ( unsigned int p = 0; p < N; ++p )
  {
    parameters[p]=static_cast<double>(p);
  }
  
  transform->SetParameters( parameters );
    
  // the following explains the problem
  int testint = ::BSplineFloor(97.999999);
  std::cerr << "::BSplineFloor(97.999999) = " << testint << std::endl;
    
  // consequently, this loop crashes at i=1
  inpoint[0]=70.0;
  inpoint[1]=50.0;
  inpoint[2]=100.0;
  for ( unsigned int i = 0; i <= 100; i++)
  {
    inpoint[2] = 97.900000 + 0.099999*static_cast<double>(i);
    std::cerr << "i =" << i << " inpoint " << inpoint[2] << " -> ";     
    outpoint = transform->TransformPoint(inpoint);
    vec=outpoint-inpoint;
    std::cerr << outpoint[2] << "\t" << "vec=" << vec << std::endl;
  }
    
  return 0;

} // end main