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

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkAnytimeImageToImageMetric.txx,v $
  Language:  C++
  Date:      $Date: 2007/09/28 21:38:43 $
  Version:   $Revision: 1.32 $

  Copyright (c) Insight Software Consortium. All rights reserved.
  See ITKCopyright.txt or http://www.itk.org/HTML/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 notices for more information.

=========================================================================*/
#ifndef _itkModImageToImageMetric_txx
#define _itkModImageToImageMetric_txx

#include "itkModImageToImageMetric.h"
#include "itkFastOrientedImage.h"


namespace itk
{

/**
 * Compute image derivatives using a central difference function
 * if we are not using a BSplineInterpolator, which includes
 * derivatives.
 */
//const FixedImageIndexType& fixedIndex,
//		const FixedImagePointType& fixedPoint,
#define FIXGRADIENT \
	if(fastOrientedImage) { \
		gradient=fastOrientedImage->TransformImageGradientToSpatialGradient(gradient); \
	}


//  virtual void ComputeImageDerivatives( const MovingImagePointType& mappedPoint,
//                                        ImageDerivativesType& gradient ) const;
//		::ComputeImageDerivatives2(   typename FixedImageSpatialSampleContainer::const_iterator fiter,
//		const MovingImagePointType& mappedPoint, 
//		ParametersType& parametergradient, const MovingImageContinuousIndexType& tempIndex) const
	template < class TFixedImage, class TMovingImage >
		void
		ModImageToImageMetric<TFixedImage,TMovingImage>
		::ComputeImageDerivatives( const MovingImagePointType& mappedPoint,
                                        GradientPixelType& gradient ) const	{
	    //const TransformJacobianType & jacobian; 
		//GradientPixelType gradient;
	    const unsigned int ParametersDimension=this->m_Transform->GetNumberOfParameters();	
		const unsigned int ImageDimension = FixedImageType::ImageDimension;
		const FastOrientedImage<FixedImagePixelType,FixedImageType::ImageDimension > * fastOrientedImage=
			dynamic_cast<const FastOrientedImage<FixedImagePixelType,FixedImageType::ImageDimension > *>(this->m_FixedImage.GetPointer());
	
		switch (m_GradientMethod) {
		case GradientImage: 
			{
			// Get the gradient by NearestNeighboorInterpolation: 
			// which is equivalent to round up the point components.
      typedef typename Superclass::OutputPointType OutputPointType;
      typedef typename OutputPointType::CoordRepType CoordRepType;
      typedef ContinuousIndex<CoordRepType,MovingImageType::ImageDimension>
        MovingImageContinuousIndexType;

      MovingImageContinuousIndexType tempIndex;
      this->m_MovingImage->TransformPhysicalPointToContinuousIndex( mappedPoint, tempIndex );


			typename MovingImageType::IndexType mappedIndex; 
			for( unsigned int j = 0; j < MovingImageType::ImageDimension; j++ )
			{
				mappedIndex[j] = static_cast<long>( vnl_math_rnd( tempIndex[j] ) );
			}
			gradient = this->GetGradientImage()->GetPixel( mappedIndex );
			FIXGRADIENT;
			}
			break;

		case BSpline:
			{
			//  if( m_InterpolatorIsBSpline )
			//  {
			// Computed moving image gradient using derivative BSpline kernel.
			gradient = m_BSplineInterpolator->EvaluateDerivative( mappedPoint );
			FIXGRADIENT;
			}
		break;
		case Calculator:
			{   
			// For all generic interpolator use central differencing.
			gradient = m_DerivativeCalculator->Evaluate( mappedPoint );
			FIXGRADIENT;
			}
			break;
		default:
			itkExceptionMacro( "Unknown method");

		}

	}



/**
 * Initialize
 */
template <class TFixedImage, class TMovingImage> 
void
ModImageToImageMetric<TFixedImage,TMovingImage>
::Initialize(void) throw ( ExceptionObject )
{
	if(m_GradientMethod==GradientImage) {
		this->m_ComputeGradient=true;
	} else {
		this->m_ComputeGradient=false;
	}
	if(m_GradientMethod==BSpline) {
		m_InterpolatorIsBSpline = true;

		BSplineInterpolatorType * testPtr = dynamic_cast<BSplineInterpolatorType *>(
			this->m_Interpolator.GetPointer() );
		if ( !testPtr )
		{
			m_InterpolatorIsBSpline = false;
			m_GradientMethod=Calculator;
			std::cerr<<"Interpolator not BSpline, using calculator for derivative"<<std::endl;
			itkWarningMacro("Interpolator not BSpline, using calculator for derivative");
			m_BSplineInterpolator = NULL;
			itkDebugMacro( "Interpolator is not BSpline" );
		} 
		else
		{
			m_BSplineInterpolator = testPtr;
			itkDebugMacro( "Interpolator is BSpline" );
		}
	}

	if(m_GradientMethod==Calculator) {
		m_DerivativeCalculator = DerivativeFunctionType::New();
		m_DerivativeCalculator->SetInputImage( this->m_MovingImage );
	}
	
	this->Superclass::Initialize();


}


  


} // end namespace itk

#endif