#ifndef __itkMultiResolutionLCCDeformableRegistration_txx
#define __itkMultiResolutionLCCDeformableRegistration_txx
#include "itkMultiResolutionLCCDeformableRegistration.h"

#include "itkRecursiveGaussianImageFilter.h"
#include "itkRecursiveMultiResolutionPyramidImageFilter.h"
#include "itkImageRegionIterator.h"
#include "vnl/vnl_math.h"
#include "itkMultiplyByConstantImageFilter.h"
#include "itkNearestNeighborInterpolateImageFunction.h"

namespace itk {

// Default constructor
template <class TFixedImage, class TMovingImage, class TField, class TRealType>
MultiResolutionLCCDeformableRegistration<TFixedImage,TMovingImage,TField,TRealType>
::MultiResolutionLCCDeformableRegistration()
{
#if (ITK_VERSION_MAJOR < 4)
    this->SetNumberOfRequiredInputs(2);
#else
  //HACK:  This really should define the names of the required inputs.
  this->SetNumberOfIndexedInputs(2);
  // Primary input is optional in this filter
  this->RemoveRequiredInputName( "Primary" );
#endif

    m_RegistrationFilter = DefaultRegistrationType::New();

    ///\todo Choose the right type of pyramid
#if ( ITK_VERSION_MAJOR == 3 && ITK_VERSION_MINOR == 12 && ITK_VERSION_PATCH == 0 )
    // Work-around for http://public.kitware.com/Bug/view.php?id=503
    itkWarningMacro("This version of ITK has a bug in MultiResolutionPyramidImageFilter - using RecursiveMultiResolutionPyramidImageFilter instead");
    typedef RecursiveMultiResolutionPyramidImageFilter<FixedImageType, FloatImageType >  ActualFixedImagePyramidType;
    typedef RecursiveMultiResolutionPyramidImageFilter<MovingImageType, FloatImageType > ActualMovingImagePyramidType;
#else
    typedef MultiResolutionPyramidImageFilter<FixedImageType, FloatImageType >  ActualFixedImagePyramidType;
    typedef MultiResolutionPyramidImageFilter<MovingImageType, FloatImageType > ActualMovingImagePyramidType;
#endif


    m_RegularizationType = 1;

    m_HarmonicWeight   = 1e-4;
    m_BendingWeight    = 1e-5;

    // Initialize standard deviations for both Gaussian regularizations
    for( unsigned int i=0; i<ImageDimension; i++ )
    {
        this->m_StandardDeviations[i]            = 1.0;
        this->m_UpdateFieldStandardDeviations[i] = 1.0;
        this->m_SimilarityCriteriaStandardDeviations[i]=10.0;

    }
    this->m_StandardDeviationWorldUnit = true;

    m_SmoothVelocityField  = true;
    m_SmoothUpdateField    = false;

    m_MovingImagePyramid   = ActualMovingImagePyramidType::New();
    m_FixedImagePyramid    = ActualFixedImagePyramidType::New();
    m_FieldExpander        = FieldExpanderType::New();
    m_ImageExpander        = ImageExpanderType::New();
    m_InitialVelocityField = NULL;

    m_NumberOfLevels       = 3;
    m_FixedImagePyramid->SetNumberOfLevels( m_NumberOfLevels );
    m_MovingImagePyramid->SetNumberOfLevels( m_NumberOfLevels );

    m_NumberOfIterations.resize( m_NumberOfLevels );
    for( unsigned int ilevel = 0; ilevel < m_NumberOfLevels; ilevel++ )
        m_NumberOfIterations[ilevel] = 10;

    m_CurrentLevel         = 0;
    m_StopRegistrationFlag = false;
    m_Exponentiator        = FieldExponentiatorType::New();
}


// Set the moving image image.
template <class TFixedImage, class TMovingImage, class TField, class TRealType>
void
MultiResolutionLCCDeformableRegistration<TFixedImage,TMovingImage,TField,TRealType>
::SetMovingImage(
  const MovingImageType * ptr )
{
#if (ITK_VERSION_MAJOR < 4)
  this->ProcessObject::SetNthInput( 2, const_cast< MovingImageType * >( ptr ) );
#else
  this->ProcessObject::SetNthInput( 2, const_cast<MovingImageType *>( ptr ) );
#endif
}


// Get the moving image image.
template <class TFixedImage, class TMovingImage, class TField, class TRealType>
const typename MultiResolutionLCCDeformableRegistration<TFixedImage,TMovingImage,TField,TRealType>
::MovingImageType *
MultiResolutionLCCDeformableRegistration<TFixedImage,TMovingImage,TField,TRealType>
::GetMovingImage(void) const
{
  return dynamic_cast< const MovingImageType * >
#if (ITK_VERSION_MAJOR < 4)
  ( this->ProcessObject::GetInput( 2 ) );
#else
  ( this->ProcessObject::GetInput( 2 ) );
#endif
}

// Set the mask image.
template <class TFixedImage, class TMovingImage, class TField, class TRealType>
void
MultiResolutionLCCDeformableRegistration<TFixedImage,TMovingImage,TField,TRealType>
::UseMask(bool flag)
{
  m_UseMask=flag;
}

// Set the mask image.
template <class TFixedImage, class TMovingImage, class TField, class TRealType>
void
MultiResolutionLCCDeformableRegistration<TFixedImage,TMovingImage,TField,TRealType>
::SetMaskImage(
  const MovingImageType * ptr )
{
  m_MaskImage=MovingImageType::New();
  m_MaskImage=ptr;
}


// Get the moving image image.
template <class TFixedImage, class TMovingImage, class TField, class TRealType>
const typename MultiResolutionLCCDeformableRegistration<TFixedImage,TMovingImage,TField,TRealType>
::MovingImageType *
MultiResolutionLCCDeformableRegistration<TFixedImage,TMovingImage,TField,TRealType>
::GetMaskImage(void) const
{
  return m_MaskImage;
}


//Set Boundary Checking

template <class TFixedImage, class TMovingImage, class TField, class TRealType>
void
MultiResolutionLCCDeformableRegistration<TFixedImage,TMovingImage,TField,TRealType>
::SetBoundaryCheck(bool flag)
{
   this->m_BoundaryCheck=flag;

}


// Set the fixed image.
template <class TFixedImage, class TMovingImage, class TField, class TRealType>
void
MultiResolutionLCCDeformableRegistration<TFixedImage,TMovingImage,TField,TRealType>
::SetFixedImage(
  const FixedImageType * ptr )
{
#if (ITK_VERSION_MAJOR < 4)
  this->ProcessObject::SetNthInput( 1, const_cast< FixedImageType * >( ptr ) );
#else
  this->ProcessObject::SetNthInput( 1, const_cast<FixedImageType *>( ptr ) );
#endif
}


// Get the fixed image.
template <class TFixedImage, class TMovingImage, class TField, class TRealType>
const typename MultiResolutionLCCDeformableRegistration<TFixedImage,TMovingImage,TField,TRealType>
::FixedImageType *
MultiResolutionLCCDeformableRegistration<TFixedImage,TMovingImage,TField,TRealType>
::GetFixedImage(void) const
{
  return dynamic_cast< const FixedImageType * >
#if (ITK_VERSION_MAJOR < 4)
    ( this->ProcessObject::GetInput( 1 ) );
#else
    ( this->ProcessObject::GetInput( 1 ) );
#endif
}


template <class TFixedImage, class TMovingImage, class TField, class TRealType>
std::vector<SmartPointer<DataObject> >::size_type
MultiResolutionLCCDeformableRegistration<TFixedImage,TMovingImage,TField,TRealType>
::GetNumberOfValidRequiredInputs() const
{
  typename std::vector<SmartPointer<DataObject> >::size_type num = 0;

  if (this->GetFixedImage())
    {
    num++;
    }

  if (this->GetMovingImage())
    {
    num++;
    }
  
  return num;
}


// Set the number of multi-resolution levels
template <class TFixedImage, class TMovingImage, class TField, class TRealType>
void
MultiResolutionLCCDeformableRegistration<TFixedImage,TMovingImage,TField,TRealType>
::SetNumberOfLevels(
  unsigned int num )
{
  if( m_NumberOfLevels != num )
    {
    this->Modified();
    m_NumberOfLevels = num;
    m_NumberOfIterations.resize( m_NumberOfLevels );
    }

  if( m_MovingImagePyramid && m_MovingImagePyramid->GetNumberOfLevels() != num )
    {
    m_MovingImagePyramid->SetNumberOfLevels( m_NumberOfLevels );
    }
  if( m_FixedImagePyramid && m_FixedImagePyramid->GetNumberOfLevels() != num )
    {
    m_FixedImagePyramid->SetNumberOfLevels( m_NumberOfLevels );
    }  
    
}


template <class TFixedImage, class TMovingImage, class TField, class TRealType>
void
MultiResolutionLCCDeformableRegistration<TFixedImage,TMovingImage,TField,TRealType>::SetRegularizationType(unsigned int regularizationType)
{
    this->m_RegularizationType=regularizationType;
}

template <class TFixedImage, class TMovingImage, class TField, class TRealType>
unsigned int
MultiResolutionLCCDeformableRegistration<TFixedImage,TMovingImage,TField,TRealType>::GetRegularizationType(void)
{
    return(this->m_RegularizationType);
}


template <class TFixedImage, class TMovingImage, class TField, class TRealType>
void
MultiResolutionLCCDeformableRegistration<TFixedImage,TMovingImage,TField,TRealType>::SetHarmonicWeight(float harmonicWeight)
{
    this->m_HarmonicWeight=harmonicWeight;
}

template <class TFixedImage, class TMovingImage, class TField, class TRealType>
float
MultiResolutionLCCDeformableRegistration<TFixedImage,TMovingImage,TField,TRealType>::GetHarmonicWeight(void)
{
    return(this->m_HarmonicWeight);
}


template <class TFixedImage, class TMovingImage, class TField, class TRealType>
void
MultiResolutionLCCDeformableRegistration<TFixedImage,TMovingImage,TField,TRealType>::SetBendingWeight(float bendingWeight)
{
    this->m_BendingWeight=bendingWeight;
}

template <class TFixedImage, class TMovingImage, class TField, class TRealType>
float
MultiResolutionLCCDeformableRegistration<TFixedImage,TMovingImage,TField,TRealType>::GetBendingWeight(void)
{
    return(this->m_BendingWeight);
}



template <class TFixedImage, class TMovingImage, class TField, class TRealType>
void
MultiResolutionLCCDeformableRegistration<TFixedImage,TMovingImage,TField,TRealType>
::SetStandardDeviationsWorldUnit( double value )
{
    double tmp[ImageDimension];
    for( int i=0; i<ImageDimension; i++ )
        tmp[i] = value;
    SetStandardDeviationsWorldUnit( tmp );
}



template <class TFixedImage, class TMovingImage, class TField, class TRealType>
void
MultiResolutionLCCDeformableRegistration<TFixedImage,TMovingImage,TField,TRealType>
::SetStandardDeviationsWorldUnit( double value[] )
{
    // Check if the array has changed
    unsigned int i;
    for( i=0; i<ImageDimension; i++ )
    {
        if( value[i] != m_StandardDeviations[i] )
            break;
    }

    // Eventually set the values
    if( i<ImageDimension )
    {
        for( i=0; i<ImageDimension; i++ )
            m_StandardDeviations[i] = value[i];
        this->Modified();
    }

    // Eventually update the unit
    if ( !this->m_StandardDeviationWorldUnit )
    {
        this->m_StandardDeviationWorldUnit = true;
        this->Modified();
    }
}



template <class TFixedImage, class TMovingImage, class TField, class TRealType>
void
MultiResolutionLCCDeformableRegistration<TFixedImage,TMovingImage,TField,TRealType>
::SetStandardDeviationsVoxelUnit( double value )
{
    double tmp[ImageDimension];
    for( int i=0; i<ImageDimension; i++ )
        tmp[i] = value;
    SetStandardDeviationsVoxelUnit( tmp );
}



template <class TFixedImage, class TMovingImage, class TField, class TRealType>
void
MultiResolutionLCCDeformableRegistration<TFixedImage,TMovingImage,TField,TRealType>
::SetStandardDeviationsVoxelUnit( double value[] )
{
    // Check if the array has changed
    unsigned int i;
    for( i=0; i<ImageDimension; i++ )
    {
        if( value[i] != m_StandardDeviations[i] )
            break;
    }

    // Eventually set the values
    if( i<ImageDimension )
    {
        for( i=0; i<ImageDimension; i++ )
            m_StandardDeviations[i] = value[i];
        this->Modified();
    }

    // Eventually update the unit
    if ( this->m_StandardDeviationWorldUnit )
    {
        this->m_StandardDeviationWorldUnit = false;
        this->Modified();
    }
}



template <class TFixedImage, class TMovingImage, class TField, class TRealType>
void
MultiResolutionLCCDeformableRegistration<TFixedImage,TMovingImage,TField,TRealType>
::SetUpdateFieldStandardDeviationsWorldUnit( double value )
{
    double tmp[ImageDimension];
    for( int i=0; i<ImageDimension; i++ )
        tmp[i] = value;
    SetUpdateFieldStandardDeviationsWorldUnit( tmp );
}



template <class TFixedImage, class TMovingImage, class TField, class TRealType>
void
MultiResolutionLCCDeformableRegistration<TFixedImage,TMovingImage,TField,TRealType>
::SetUpdateFieldStandardDeviationsWorldUnit( double value[] )
{
    // Check if the array has changed
    unsigned int i;
    for( i=0; i<ImageDimension; i++ )
    {
        if( value[i] != m_UpdateFieldStandardDeviations[i] )
            break;
    }

    // Eventually set the values
    if( i<ImageDimension )
    {
        for( i=0; i<ImageDimension; i++ )
            m_UpdateFieldStandardDeviations[i] = value[i];
        this->Modified();
    }

    // Eventually update the unit
    if ( !this->m_StandardDeviationWorldUnit )
    {
        this->m_StandardDeviationWorldUnit = true;
        this->Modified();
    }
}



template <class TFixedImage, class TMovingImage, class TField, class TRealType>
void
MultiResolutionLCCDeformableRegistration<TFixedImage,TMovingImage,TField,TRealType>
::SetUpdateFieldStandardDeviationsVoxelUnit( double value )
{
    double tmp[ImageDimension];
    for( int i=0; i<ImageDimension; i++ )
        tmp[i] = value;
    SetUpdateFieldStandardDeviationsVoxelUnit( tmp );
}



template <class TFixedImage, class TMovingImage, class TField, class TRealType>
void
MultiResolutionLCCDeformableRegistration<TFixedImage,TMovingImage,TField,TRealType>
::SetUpdateFieldStandardDeviationsVoxelUnit( double value[] )
{
    // Check if the array has changed
    unsigned int i;
    for( i=0; i<ImageDimension; i++ )
    {
        if( value[i] != m_UpdateFieldStandardDeviations[i] )
            break;
    }

    // Eventually set the values
    if( i<ImageDimension )
    {
        for( i=0; i<ImageDimension; i++ )
            m_UpdateFieldStandardDeviations[i] = value[i];
        this->Modified();
    }

    // Eventually update the unit
    if ( this->m_StandardDeviationWorldUnit )
    {
        this->m_StandardDeviationWorldUnit = false;
        this->Modified();
    }
}


template <class TFixedImage, class TMovingImage, class TField, class TRealType>
void
MultiResolutionLCCDeformableRegistration<TFixedImage,TMovingImage,TField,TRealType>
::SetSimilarityCriteriaStandardDeviationsWorldUnit( double value )
{
    double tmp[ImageDimension];
    for( int i=0; i<ImageDimension; i++ )
        tmp[i] = value;
    SetSimilarityCriteriaStandardDeviationsWorldUnit( tmp );
}

template <class TFixedImage, class TMovingImage, class TField, class TRealType>
void
MultiResolutionLCCDeformableRegistration<TFixedImage,TMovingImage,TField,TRealType>
::SetSimilarityCriteriaStandardDeviationsWorldUnit( double value[] )
{
    // Check if the array has changed
    unsigned int i;
    for( i=0; i<ImageDimension; i++ )
    {
        if( value[i] != m_SimilarityCriteriaStandardDeviations[i] )
            break;
    }

    // Eventually set the values
    if( i<ImageDimension )
    {
        for( i=0; i<ImageDimension; i++ )
            m_SimilarityCriteriaStandardDeviations[i] = value[i];
        this->Modified();
    }

    // Eventually update the unit
    if ( !this->m_StandardDeviationWorldUnit )
    {
        this->m_StandardDeviationWorldUnit = true;
        this->Modified();
    }
}



template <class TFixedImage, class TMovingImage, class TField, class TRealType>
void
MultiResolutionLCCDeformableRegistration<TFixedImage,TMovingImage,TField,TRealType>
::SetSimilarityCriteriaStandardDeviationsVoxelUnit( double value )
{
    double tmp[ImageDimension];
    for( int i=0; i<ImageDimension; i++ )
        tmp[i] = value;
    SetSimilarityCriteriaStandardDeviationsVoxelUnit( tmp );
}

template <class TFixedImage, class TMovingImage, class TField, class TRealType>
void
MultiResolutionLCCDeformableRegistration<TFixedImage,TMovingImage,TField,TRealType>
::SetSimilarityCriteriaStandardDeviationsVoxelUnit( double value[] )
{
    // Check if the array has changed
    unsigned int i;
    for( i=0; i<ImageDimension; i++ )
    {
        if( value[i] != m_SimilarityCriteriaStandardDeviations[i] )
            break;
    }

    // Eventually set the values
    if( i<ImageDimension )
    {
        for( i=0; i<ImageDimension; i++ )
            m_SimilarityCriteriaStandardDeviations[i] = value[i];
        this->Modified();
    }

    // Eventually update the unit
    if ( this->m_StandardDeviationWorldUnit )
    {
        this->m_StandardDeviationWorldUnit = false;
        this->Modified();
    }
}

template <class TFixedImage, class TMovingImage, class TField, class TRealType>
void
MultiResolutionLCCDeformableRegistration<TFixedImage,TMovingImage,TField,TRealType>
::SetSigmaI( double value )
{
  this->m_SigmaI=value;   
}




// Standard PrintSelf method.
template <class TFixedImage, class TMovingImage, class TField, class TRealType>
void
MultiResolutionLCCDeformableRegistration<TFixedImage,TMovingImage,TField,TRealType>
::PrintSelf(std::ostream& os, Indent indent) const
{
  Superclass::PrintSelf(os, indent);
  os << indent << "NumberOfLevels: " << m_NumberOfLevels << std::endl;
  os << indent << "CurrentLevel: " << m_CurrentLevel << std::endl;

  os << indent << "NumberOfIterations: [";
  unsigned int ilevel;
  for( ilevel = 0; ilevel < m_NumberOfLevels - 1; ilevel++ )
    {
    os << m_NumberOfIterations[ilevel] << ", ";
    }
  os << m_NumberOfIterations[ilevel] << "]" << std::endl;
  
  os << indent << "RegistrationFilter: ";
  os << m_RegistrationFilter.GetPointer() << std::endl;
  os << indent << "MovingImagePyramid: ";
  os << m_MovingImagePyramid.GetPointer() << std::endl;
  os << indent << "FixedImagePyramid: ";
  os << m_FixedImagePyramid.GetPointer() << std::endl;

  os << indent << "FieldExpander: ";
  os << m_FieldExpander.GetPointer() << std::endl;

  os << indent << "StopRegistrationFlag: ";
  os << m_StopRegistrationFlag << std::endl;
  
  os << indent << "Exponentiator: ";
  os << m_Exponentiator << std::endl;

}

/*
 * Perform a the deformable registration using a multiresolution scheme
 * using an internal mini-pipeline
 *
 *  ref_pyramid ->  registrator  ->  field_expander --|| tempField
 * test_pyramid ->           |                              |
 *                           |                              |
 *                           --------------------------------
 *
 * A tempField image is used to break the cycle between the
 * registrator and field_expander.*/
template <class TFixedImage, class TMovingImage, class TField, class TRealType>
void
MultiResolutionLCCDeformableRegistration<TFixedImage,TMovingImage,TField,TRealType>
::GenerateData()
{
    // Check for NULL images and pointers
    MovingImageConstPointer movingImage = this->GetMovingImage();
    FixedImageConstPointer  fixedImage = this->GetFixedImage();

    if( !movingImage || !fixedImage )
    {
        itkExceptionMacro( << "Fixed and/or moving image not set" );
    }

    if( !m_MovingImagePyramid || !m_FixedImagePyramid )
    {
        itkExceptionMacro( << "Fixed and/or moving pyramid not set" );
    }

    if( !m_RegistrationFilter )
    {
        itkExceptionMacro( << "Registration filter not set" );
    }
#if (ITK_VERSION_MAJOR < 4)
    if( this->m_InitialVelocityField && this->GetInput(0) )
#else
  if( this->m_InitialVelocityField && this->GetInput(VELOCITYFIELD_IMAGE_CODE) )
#endif
    {
        itkExceptionMacro( << "Only one initial velocity can be given. "
                           << "SetInitialVelocityField should not be used in "
                           << "cunjunction with SetArbitraryInitialVelocityField "
                           << "or SetInput.");
    }

    // Create the image pyramids.
    m_MovingImagePyramid->SetInput( movingImage );
    m_MovingImagePyramid->UpdateLargestPossibleRegion();

    m_FixedImagePyramid->SetInput( fixedImage );
    m_FixedImagePyramid->UpdateLargestPossibleRegion();

    // Initializations
    m_CurrentLevel = 0;
    m_StopRegistrationFlag = false;

    unsigned int movingLevel = std::min( (int) m_CurrentLevel,
                                             (int) m_MovingImagePyramid->GetNumberOfLevels() );

    unsigned int fixedLevel = std::min( (int) m_CurrentLevel,
                                            (int) m_FixedImagePyramid->GetNumberOfLevels() );

    VelocityFieldPointer tempField = nullptr;

    VelocityFieldPointer inputPtr =
#if (ITK_VERSION_MAJOR < 4)
            const_cast< VelocityFieldType * >( this->GetInput(0) );
#else
    const_cast<VelocityFieldType *>( this->GetInput(VELOCITYFIELD_IMAGE_CODE) );
#endif

    if ( this->m_InitialVelocityField )
    {
        tempField = this->m_InitialVelocityField;
    }
    else if( inputPtr )
    {
        // Arbitrary initial velocity field is set.
        // smooth it and resample

        // First smooth it
        tempField = inputPtr;

        typedef RecursiveGaussianImageFilter< VelocityFieldType, VelocityFieldType> GaussianFilterType;
        typename GaussianFilterType::Pointer smoother = GaussianFilterType::New();

        for (unsigned int dim=0; dim<VelocityFieldType::ImageDimension; ++dim)
        {
            // sigma accounts for the subsampling of the pyramid
            double sigma = 0.5 * static_cast<float>( m_FixedImagePyramid->GetSchedule()[fixedLevel][dim] );

            // but also for a possible discrepancy in the spacing
            sigma *= fixedImage->GetSpacing()[dim] / inputPtr->GetSpacing()[dim];

            smoother->SetInput( tempField );
            smoother->SetSigma( sigma );
            smoother->SetDirection( dim );

            smoother->Update();

            tempField = smoother->GetOutput();
            tempField->DisconnectPipeline();
        }


        // Now resample
        m_FieldExpander->SetInput( tempField );
	//std::cout<<"Fixed Level "<<fixedLevel<<std::endl;
        typename FloatImageType::Pointer fi = m_FixedImagePyramid->GetOutput( fixedLevel );
        m_FieldExpander->SetSize(             fi->GetLargestPossibleRegion().GetSize() );
        m_FieldExpander->SetOutputStartIndex( fi->GetLargestPossibleRegion().GetIndex() );
        m_FieldExpander->SetOutputOrigin(     fi->GetOrigin() );
        m_FieldExpander->SetOutputSpacing(    fi->GetSpacing());
        m_FieldExpander->SetOutputDirection(  fi->GetDirection());
        m_FieldExpander->UpdateLargestPossibleRegion();

        tempField = m_FieldExpander->GetOutput();
        tempField->DisconnectPipeline();
    }

    if (m_UseMask)
      {
       m_ImageExpander->SetInput( m_MaskImage );
       m_ImageExpander->SetShrinkFactors(pow(2.0,(int)(m_MovingImagePyramid->GetNumberOfLevels()-movingLevel-1)));
       m_ImageExpander->Update();
       m_RegistrationFilter->UseMask(m_UseMask);
       m_RegistrationFilter->SetMaskImage(m_ImageExpander->GetOutput() );
      }
    else m_RegistrationFilter->UseMask(m_UseMask);

    bool lastShrinkFactorsAllOnes = false;

    m_RegistrationFilter->SetRegularizationType(this->m_RegularizationType);

    if (m_RegularizationType==0)
      {
        // Mode for smoothing or not the velocity field
        if ( this->m_SmoothVelocityField )
           m_RegistrationFilter->SmoothVelocityFieldOn();
        else
           m_RegistrationFilter->SmoothVelocityFieldOff();

        // Mode for smoothing or not the velocity field
        if ( this->m_SmoothUpdateField )
          m_RegistrationFilter->SmoothUpdateFieldOn();
        else
          m_RegistrationFilter->SmoothUpdateFieldOff();

        // Sets the filter standard deviation for global field regularization (voxel unit only)
        if ( !m_StandardDeviationWorldUnit )
           m_RegistrationFilter->SetStandardDeviationsVoxelUnit( m_StandardDeviations );

        // Sets the filter standard deviation for update field regularization (voxel unit only)
        if ( !m_StandardDeviationWorldUnit )
           m_RegistrationFilter->SetUpdateFieldStandardDeviationsVoxelUnit( m_UpdateFieldStandardDeviations );
       }

    else if (m_RegularizationType==1)
       {
        // Set penalization for the harmonic energy term
        m_RegistrationFilter->SetHarmonicWeight(this->m_HarmonicWeight);

        // Set penalization for the bending energy term
        m_RegistrationFilter->SetBendingWeight(this->m_BendingWeight);
       }

        // Sets the filter standard deviation for the similarity criterion (voxel unit only)
    m_RegistrationFilter->SetSimilarityCriteriaStandardDeviationsVoxelUnit( m_SimilarityCriteriaStandardDeviations );


    m_RegistrationFilter->SetSigmaI( this->m_SigmaI );

    m_RegistrationFilter->SetBoundaryCheck(this->m_BoundaryCheck);
    // Loop
    while ( !this->Halt() )
    {

        if( tempField.IsNull() )
        {
            m_RegistrationFilter->SetInitialVelocityField( NULL );
        }
        else
        {
            // Resample the field to be the same size as the fixed image at the current level
            m_FieldExpander->SetInput( tempField );

            typename FloatImageType::Pointer fi = m_FixedImagePyramid->GetOutput( fixedLevel );
            m_FieldExpander->SetSize(             fi->GetLargestPossibleRegion().GetSize() );
            m_FieldExpander->SetOutputStartIndex( fi->GetLargestPossibleRegion().GetIndex() );
            m_FieldExpander->SetOutputOrigin(     fi->GetOrigin() );
            m_FieldExpander->SetOutputSpacing(    fi->GetSpacing());
            m_FieldExpander->SetOutputDirection(  fi->GetDirection());

            m_FieldExpander->UpdateLargestPossibleRegion();
            tempField = m_FieldExpander->GetOutput();
            tempField->DisconnectPipeline();

	    
            m_RegistrationFilter->SetInitialVelocityField( tempField );

            // Resample the mask image to be the same size as the fixed image at the current level
            if (m_UseMask)

            {
              m_ImageExpander->SetInput( m_MaskImage );

    /*        m_ImageExpander->SetSize(             fi->GetLargestPossibleRegion().GetSize() );
              m_ImageExpander->SetOutputStartIndex( fi->GetLargestPossibleRegion().GetIndex() );
              m_ImageExpander->SetOutputOrigin(     fi->GetOrigin() );
              m_ImageExpander->SetOutputSpacing(    fi->GetSpacing());
              m_ImageExpander->SetOutputDirection(  fi->GetDirection());
*/
              m_ImageExpander->SetShrinkFactors(pow(2.0,(int)(m_MovingImagePyramid->GetNumberOfLevels()-movingLevel-1)));
              m_ImageExpander->UpdateLargestPossibleRegion();
           
              m_RegistrationFilter->UseMask( m_UseMask );
              m_RegistrationFilter->SetMaskImage( m_ImageExpander->GetOutput() );
            }
            else
              m_RegistrationFilter->UseMask( m_UseMask );

        }


        // Setup registration filter and pyramids
        m_RegistrationFilter->SetMovingImage(        m_MovingImagePyramid->GetOutput(movingLevel) );
        m_RegistrationFilter->SetFixedImage(         m_FixedImagePyramid->GetOutput(fixedLevel) );
        m_RegistrationFilter->SetNumberOfIterations( m_NumberOfIterations[m_CurrentLevel] );


        // Resolution
        double resolution = 1;
        for (int i=0; i<ImageDimension; i++)
        {
            double s1  = fixedImage->GetLargestPossibleRegion().GetSize()[i];
            double s2  = m_FixedImagePyramid->GetOutput(movingLevel)->GetLargestPossibleRegion().GetSize()[i];
            resolution = std::max( resolution, s1/s2 );
        }


        // Sets the filter standard deviation for global field regularization (world unit only)
        if ( m_StandardDeviationWorldUnit )
        {
            double sd[ImageDimension];
            for (int i=0; i<ImageDimension; i++)
              sd[i] = m_StandardDeviations[i];

            m_RegistrationFilter->SetStandardDeviationsWorldUnit( sd );
        }


        // Sets the filter standard deviation for update field regularization (world unit only)
        if ( m_StandardDeviationWorldUnit )
        {
            double sd[ImageDimension];
            for (int i=0; i<ImageDimension; i++)
              sd[i] = m_UpdateFieldStandardDeviations[i];

            m_RegistrationFilter->SetUpdateFieldStandardDeviationsWorldUnit( sd );
        }

     // Sets the filter standard deviation for the similarity criterion (world unit only)
        if ( m_StandardDeviationWorldUnit )
        {
            double sd[ImageDimension];
            for (int i=0; i<ImageDimension; i++)
                sd[i] = m_SimilarityCriteriaStandardDeviations[i] * resolution;

            m_RegistrationFilter->SetSimilarityCriteriaStandardDeviationsWorldUnit( sd );
        }

        // cache shrink factors for computing the next expand factors.
        lastShrinkFactorsAllOnes = true;
        for( unsigned int idim = 0; idim < ImageDimension; idim++ )
        {
            if ( m_FixedImagePyramid->GetSchedule()[fixedLevel][idim] > 1 )
            {
                lastShrinkFactorsAllOnes = false;
                break;
            }
        }

        // compute new velocity field
        m_RegistrationFilter->UpdateLargestPossibleRegion();
       // std::cout<<"Smoothing velocity!"<<std::endl;
       // m_RegistrationFilter->SmoothVelocityField();

        //Workaround for keep multithreading while regularing in Fourier domain
        tempField = m_RegistrationFilter->GetOutput();

        tempField->DisconnectPipeline();

	
	
        // Increment level counter.
        m_CurrentLevel++;
        movingLevel = std::min( (int) m_CurrentLevel,
                                    (int) m_MovingImagePyramid->GetNumberOfLevels() );
        fixedLevel = std::min( (int) m_CurrentLevel,
                                   (int) m_FixedImagePyramid->GetNumberOfLevels() );

        // Invoke an iteration event.
        this->InvokeEvent( IterationEvent() );

        // We can release data from pyramid which are no longer required.
        if ( movingLevel > 0 )
        {
            m_MovingImagePyramid->GetOutput( movingLevel - 1 )->ReleaseData();
        }
        if( fixedLevel > 0 )
        {
            m_FixedImagePyramid->GetOutput( fixedLevel - 1 )->ReleaseData();
        }

    } // while not Halt()

    if( !lastShrinkFactorsAllOnes )
    {
        // Some of the last shrink factors are not one
        // graft the output of the expander filter to
        // to output of this filter

        // resample the field to the same size as the fixed image
        m_FieldExpander->SetInput( tempField );
        m_FieldExpander->SetSize(
                fixedImage->GetLargestPossibleRegion().GetSize() );
        m_FieldExpander->SetOutputStartIndex(
                fixedImage->GetLargestPossibleRegion().GetIndex() );
        m_FieldExpander->SetOutputOrigin( fixedImage->GetOrigin() );
        m_FieldExpander->SetOutputSpacing( fixedImage->GetSpacing());
        m_FieldExpander->SetOutputDirection( fixedImage->GetDirection());

        m_FieldExpander->UpdateLargestPossibleRegion();
////MINE
	typedef typename itk::MultiplyByConstantImageFilter<VelocityFieldType,float, VelocityFieldType> VelocityMultiplierType;
	    typename VelocityMultiplierType::Pointer VelocityMultiplier=VelocityMultiplierType::New();
	    
	    VelocityMultiplier->SetInput(m_FieldExpander->GetOutput() );
	    VelocityMultiplier->SetConstant(2);
	    VelocityMultiplier->UpdateLargestPossibleRegion();	


        this->GraftOutput( VelocityMultiplier->GetOutput() );
    }
    else
    {
        // all the last shrink factors are all ones
        // graft the output of registration filter to
        // to output of this filter

	typedef typename itk::MultiplyByConstantImageFilter<VelocityFieldType,float, VelocityFieldType> VelocityMultiplierType;
	    typename VelocityMultiplierType::Pointer VelocityMultiplier=VelocityMultiplierType::New();
	    
	    VelocityMultiplier->SetInput(tempField);
	    VelocityMultiplier->SetConstant(2);
	    VelocityMultiplier->UpdateLargestPossibleRegion();	

            tempField = VelocityMultiplier->GetOutput();
            tempField->DisconnectPipeline();

        this->GraftOutput( tempField );
    }

    // Release memory
    m_FieldExpander->SetInput( NULL );
    m_FieldExpander->GetOutput()->ReleaseData();
    m_RegistrationFilter->SetInput( NULL );
    m_RegistrationFilter->GetOutput()->ReleaseData();

}


template <class TFixedImage, class TMovingImage, class TField, class TRealType>
void
MultiResolutionLCCDeformableRegistration<TFixedImage,TMovingImage,TField,TRealType>
::StopRegistration()
{
  m_RegistrationFilter->StopRegistration();
  m_StopRegistrationFlag = true;
}

template <class TFixedImage, class TMovingImage, class TField, class TRealType>
bool
MultiResolutionLCCDeformableRegistration<TFixedImage,TMovingImage,TField,TRealType>
::Halt()
{
  // Halt the registration after the user-specified number of levels
  if (m_NumberOfLevels != 0)
    {
    this->UpdateProgress( static_cast<float>( m_CurrentLevel ) /
                          static_cast<float>( m_NumberOfLevels ) );
    }

  if ( m_CurrentLevel >= m_NumberOfLevels )
    {
    return true;
    }
  if ( m_StopRegistrationFlag )
    {
    return true;
    }
  else
    { 
    return false; 
    }

}


template <class TFixedImage, class TMovingImage, class TField, class TRealType>
void
MultiResolutionLCCDeformableRegistration<TFixedImage,TMovingImage,TField,TRealType>
::GenerateOutputInformation()
{

  typename DataObject::Pointer output;

#if (ITK_VERSION_MAJOR < 4)
  if( this->GetInput(0) )
#else
  if( this->GetInput(0) )
#endif
    {
    // Initial velocity field is set.
    // Copy information from initial field.
    this->Superclass::GenerateOutputInformation();

    }
  else if( this->GetFixedImage() )
    {
    // Initial deforamtion field is not set. 
    // Copy information from the fixed image.
    for (unsigned int idx = 0; idx < 
           this->GetNumberOfOutputs(); ++idx )
      {
      output = this->GetOutput(idx);
      if (output)
        {
        output->CopyInformation(this->GetFixedImage());
        }  
      }

    }

}


template <class TFixedImage, class TMovingImage, class TField, class TRealType>
void
MultiResolutionLCCDeformableRegistration<TFixedImage,TMovingImage,TField,TRealType>
::GenerateInputRequestedRegion()
{

  // call the superclass's implementation
  Superclass::GenerateInputRequestedRegion();

  // request the largest possible region for the moving image
  MovingImagePointer movingPtr = 
    const_cast< MovingImageType * >( this->GetMovingImage() );
  if( movingPtr )
    {
    movingPtr->SetRequestedRegionToLargestPossibleRegion();
    }
  
  // just propagate up the output requested region for
  // the fixed image and initial velocity field.
  VelocityFieldPointer inputPtr = 
#if (ITK_VERSION_MAJOR < 4)
    const_cast< VelocityFieldType * >( this->GetInput(0) );
#else
    const_cast<VelocityFieldType *>( this->GetInput(VELOCITYFIELD_IMAGE_CODE) );
#endif
  VelocityFieldPointer outputPtr = this->GetOutput();
  FixedImagePointer fixedPtr = 
    const_cast< FixedImageType *>( this->GetFixedImage() );

  if( inputPtr )
    {
    inputPtr->SetRequestedRegion( outputPtr->GetRequestedRegion() );
    }

  if( fixedPtr )
    {
    fixedPtr->SetRequestedRegion( outputPtr->GetRequestedRegion() );
    }

}


template <class TFixedImage, class TMovingImage, class TField, class TRealType>
void
MultiResolutionLCCDeformableRegistration<TFixedImage,TMovingImage,TField,TRealType>
::EnlargeOutputRequestedRegion(
  DataObject * ptr )
{
  // call the superclass's implementation
  Superclass::EnlargeOutputRequestedRegion( ptr );

  // set the output requested region to largest possible.
  VelocityFieldType * outputPtr;
  outputPtr = dynamic_cast<VelocityFieldType*>( ptr );

  if( outputPtr )
    {
    outputPtr->SetRequestedRegionToLargestPossibleRegion();
    }

}


template <class TFixedImage, class TMovingImage, class TField, class TRealType>
typename MultiResolutionLCCDeformableRegistration<TFixedImage,TMovingImage,TField,TRealType>
::DeformationFieldPointer
MultiResolutionLCCDeformableRegistration<TFixedImage,TMovingImage,TField,TRealType>
::GetDeformationField()
{
  //std::cout<<"MultiResolutionLCCDeformableRegistration::GetDeformationField"<<std::endl;
  m_Exponentiator->SetInput( this->GetVelocityField() );
  m_Exponentiator->ComputeInverseOff();
  m_Exponentiator->Update();
  DeformationFieldPointer field = m_Exponentiator->GetOutput();
  field->DisconnectPipeline();
  return field;
}


template <class TFixedImage, class TMovingImage, class TField, class TRealType>
typename MultiResolutionLCCDeformableRegistration<TFixedImage,TMovingImage,TField,TRealType>
::DeformationFieldPointer
MultiResolutionLCCDeformableRegistration<TFixedImage,TMovingImage,TField,TRealType>
::GetInverseDeformationField()
{
  //std::cout<<"MultiResolutionLCCDeformableRegistration::GetInverseDeformationField"<<std::endl;
  m_Exponentiator->SetInput( this->GetVelocityField() );
  m_Exponentiator->ComputeInverseOn();
  m_Exponentiator->Update();
  DeformationFieldPointer field = m_Exponentiator->GetOutput();
  field->DisconnectPipeline();
  // Reset compute inverse back to off to avoid some broder effects
  m_Exponentiator->ComputeInverseOff();
  return field;
}


} // end namespace itk

#endif