/*===================================================================
  Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv

  Authors belong to:
  - University of LYON              http://www.universite-lyon.fr/
  - Léon Bérard cancer center       http://www.centreleonberard.fr
  - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr

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

  It is distributed under dual licence

  - BSD        See included LICENSE.txt file
  - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
===========================================================================**/
#ifndef __clitkMultipleBSplineDeformableTransformInitializer_txx
#define __clitkMultipleBSplineDeformableTransformInitializer_txx
#include "clitkMultipleBSplineDeformableTransformInitializer.h"

namespace clitk
{

  template < class TTransform, class TImage >
  MultipleBSplineDeformableTransformInitializer<TTransform, TImage >
  ::MultipleBSplineDeformableTransformInitializer()
  {
    this->m_NumberOfControlPointsInsideTheImage.Fill( 5 );
    this->m_ControlPointSpacing.Fill(64.);
    this->m_ChosenSpacing.Fill(64.);
    m_NumberOfControlPointsIsGiven=false;
    m_ControlPointSpacingIsGiven=false;
    m_TransformRegionIsGiven=false;
    m_SamplingFactorIsGiven=false;
    m_SplineOrders.Fill(3);
    m_Parameters=NULL;
  }

  template < class TTransform, class TImage >
  void
  MultipleBSplineDeformableTransformInitializer<TTransform, TImage >
  ::InitializeTransform()
  {
    // Sanity check:
    // The image is required for the region and spacing
    if( ! this->m_Image )
      {
	itkExceptionMacro( "Reference Image has not been set" );
	return;
      }

    // A pointer to the transform is required
    if( ! this->m_Transform )
      {
	itkExceptionMacro( "Transform has not been set" );
	return;
      }

    // If the image come from a filter, then update that filter.
    if( this->m_Image->GetSource() )
      {
	this->m_Image->GetSource()->Update();
      }


    //--------------------------------------
    // Spacing & Size on image
    //--------------------------------------
    SpacingType fixedImageSpacing = m_Image->GetSpacing();
    SizeType fixedImageSize=m_Image->GetLargestPossibleRegion().GetSize();
    typename RegionType::SizeType   gridBorderSize;
    typename RegionType::SizeType   totalGridSize;

    // Only spacing given: adjust if necessary
    if (m_ControlPointSpacingIsGiven  && !m_NumberOfControlPointsIsGiven)
      {
	for(unsigned int r=0; r<InputDimension; r++)
	  {
	    // JV
	    m_ChosenSpacing[r]= m_ControlPointSpacing[r];
	    m_ControlPointSpacing[r]= ( itk::Math::Round<double>(m_ChosenSpacing[r]/fixedImageSpacing[r]) *fixedImageSpacing[r] ) ;
	    m_NumberOfControlPointsInsideTheImage[r] = ceil( (double)fixedImageSize[r]*fixedImageSpacing[r]/ m_ControlPointSpacing[r] );
	    if (  ( ceil( (double)fixedImageSize[r]*fixedImageSpacing[r]/ m_ControlPointSpacing[r] ) )
		 == ( (double)fixedImageSize[r]*fixedImageSpacing[r]/ m_ControlPointSpacing[r] ) )
	      {
		m_NumberOfControlPointsInsideTheImage[r] +=1;
		DD("Adding control point");
	      }
	  }
      }

    // Only number of CP given: adjust if necessary
    // JV -1 ?
    else if (!m_ControlPointSpacingIsGiven  && m_NumberOfControlPointsIsGiven)
      {
	for(unsigned int r=0; r<InputDimension; r++)
	  {
	    m_ChosenSpacing[r]=fixedImageSpacing[r]*( (double)(fixedImageSize[r])  /
						      (double)(m_NumberOfControlPointsInsideTheImage[r]) );
	    m_ControlPointSpacing[r]= fixedImageSpacing[r]* ceil( (double)(fixedImageSize[r] - 1)  /
								  (double)(m_NumberOfControlPointsInsideTheImage[r] - 1) );
	  }
      }

    // Both or none of Spacing and number of CP given: no adjustment adjust, just warnings
    else
      {
	for(unsigned int r=0; r<InputDimension; r++)
	  {
	    m_ChosenSpacing[r]= m_ControlPointSpacing[r];
	    if (m_NumberOfControlPointsInsideTheImage[r]*m_ControlPointSpacing[r]< fixedImageSize[r]*fixedImageSpacing[r])
	      {
		std::cout<<"WARNING: Specified control point region ("<<m_NumberOfControlPointsInsideTheImage[r]*m_ControlPointSpacing[r]
			 <<"mm) does not cover the fixed image region ("<< fixedImageSize[r]*fixedImageSpacing[r]
			 <<"mm) for dimension "<<r<<"!" <<std::endl
			 <<"Specify only --spacing or --control for automatic adjustment..."<<std::endl;
	      }
	    if (  fmod(m_ControlPointSpacing[r], fixedImageSpacing[r]) )
	      {
		std::cout<<"WARNING: Specified control point spacing for dimension "<<r
			 <<" does not allow exact representation of BLUT FFD!"<<std::endl
			 <<"Spacing ratio is non-integer: "<<m_ControlPointSpacing[r]/ fixedImageSpacing[r]<<std::endl
			 <<"Specify only --spacing or --control for automatic adjustment..."<<std::endl;
	      }
	  }
      }

    if (m_Verbose) std::cout<<"The chosen control point spacing "<<m_ChosenSpacing<<"..."<<std::endl;
    if (m_Verbose) std::cout<<"The control points spacing was adapted to "<<m_ControlPointSpacing<<"..."<<std::endl;
    if (m_Verbose) std::cout<<"The number of (internal) control points is "<<m_NumberOfControlPointsInsideTheImage<<"..."<<std::endl;


    //--------------------------------------
    // Border size
    //--------------------------------------
    for(unsigned int r=0; r<InputDimension; r++) gridBorderSize[r]=m_SplineOrders[r]; // Border for spline order = 3 ( 1 lower, 2 upper )
    totalGridSize = m_NumberOfControlPointsInsideTheImage + gridBorderSize;
    if (m_Verbose) std::cout<<"The total control point grid size was set to "<<totalGridSize<<"..."<<std::endl;
    RegionType gridRegion;
    gridRegion.SetSize(totalGridSize);


    //--------------------------------------
    // Origin
    //--------------------------------------
    typedef typename TransformType::OriginType OriginType;
    OriginType fixedImageOrigin, gridOrigin;

    // In case of non-zero index
    m_Image->TransformIndexToPhysicalPoint(m_Image->GetLargestPossibleRegion().GetIndex(),fixedImageOrigin);
    typename ImageType::DirectionType gridDirection = m_Image->GetDirection();
    SizeType orderShift;

    // Shift depends on order
    for(unsigned int r=0; r<InputDimension; r++)
      {
	orderShift[r] = m_SplineOrders[r] / 2;
	gridOrigin[r] = fixedImageOrigin[r]- (double) orderShift[r]* m_ControlPointSpacing[r];
      }
    if (m_Verbose) std::cout<<"The control point grid origin was set to "<<gridOrigin<<"..."<<std::endl;


    //--------------------------------------
    // LUT sampling factors
    //--------------------------------------
    for (unsigned int i=0; i< InputDimension; i++)
      {
	if (m_Verbose) std::cout<<"For dimension "<<i<<", the ideal sampling factor (if integer) is a multitude of "
				<< m_ControlPointSpacing[i]/ fixedImageSpacing[i]<<"..."<<std::endl;
	if ( !m_SamplingFactorIsGiven) m_SamplingFactors[i]= (int) ( m_ControlPointSpacing[i]/ fixedImageSpacing[i]);
	if (m_Verbose) std::cout<<"Setting sampling factor "<<i<<" to "<<m_SamplingFactors[i]<<"..."<<std::endl;
      }


    //--------------------------------------
    // Set
    //--------------------------------------
    this->m_Transform->SetSplineOrders(m_SplineOrders);
    this->m_Transform->SetGridRegion( gridRegion );
    this->m_Transform->SetGridOrigin( gridOrigin );
    this->m_Transform->SetGridSpacing( m_ControlPointSpacing );
    this->m_Transform->SetGridDirection( gridDirection );
    this->m_Transform->SetLUTSamplingFactors(m_SamplingFactors);

  }

  template < class TTransform, class TImage >
  void
  MultipleBSplineDeformableTransformInitializer<TTransform, TImage >
  ::SetInitialParameters( const std::string& s, ParametersType& params)
  {
    //---------------------------------------
    // Read Initial parameters
    //---------------------------------------
    typedef itk::ImageFileReader<CoefficientImageType> CoefficientReaderType;
    typename CoefficientReaderType::Pointer coeffReader = CoefficientReaderType::New();
    std::vector<typename CoefficientImageType::Pointer> coefficientImages;
    unsigned nLabels = m_Transform->GetnLabels();
    coefficientImages.resize(nLabels);

    int dotpos = s.length() - 1;
    while (dotpos >= 0 && s[dotpos] != '.')
      dotpos--;

    for (unsigned i = 0; i < nLabels; ++i)
    {
      std::ostringstream osfname;
      osfname << s.substr(0, dotpos) << '_' << i << s.substr(dotpos);
      coeffReader->SetFileName(osfname.str());
      if (m_Verbose)
        std::cout << "Reading initial coefficients from file: " << osfname.str() << "..." << std::endl;
      coeffReader->Update();
      coefficientImages[i] = coeffReader->GetOutput();
    }
    this->SetInitialParameters(coefficientImages, params);
  }

  template < class TTransform, class TImage >
  void
  MultipleBSplineDeformableTransformInitializer<TTransform, TImage >
  ::SetInitialParameters(std::vector<typename CoefficientImageType::Pointer> coefficientImages, ParametersType& params)
  {
    // Keep a reference
    m_Parameters=&params;

    // Resample
    typedef ResampleBSplineDeformableTransformImageFilter<CoefficientImageType,CoefficientImageType> ResamplerType;
    typename ResamplerType::Pointer resampler;
    unsigned int i = 0;
    for (unsigned l = 0; l < coefficientImages.size(); ++l)
    {
      typename ResamplerType::Pointer resampler = ResamplerType::New();
      resampler->SetInput(coefficientImages[l]);
      resampler->SetOutputParametersFromImage(m_Transform->GetCoefficientImages()[l]);
      if (m_Verbose)
        std::cout << "Resampling initial coefficients..." << std::endl;
      resampler->Update();

      // Copy parameters into the existing array, I know its crappy
      typedef itk::ImageRegionConstIterator<CoefficientImageType> Iterator;
      Iterator it(resampler->GetOutput(), resampler->GetOutput()->GetLargestPossibleRegion());
      it.GoToBegin();
      while (!it.IsAtEnd())
      {
        for (unsigned int j = 0; j < InputDimension; ++j)
          params[i + j] = it.Get()[j];
        ++it;
        i += InputDimension;
      }
    }

    // JV pass the reference to the array !!
    this->m_Transform->SetParameters(params);
  }



}  // namespace itk

#endif