itkOptImageToImageMetric.txx

Go to the documentation of this file.

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

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkOptImageToImageMetric.txx,v $
  Language:  C++
  Date:      $Date: 2010-03-27 20:54:43 $
  Version:   $Revision: 1.40 $

  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 __itkOptImageToImageMetric_txx
#define __itkOptImageToImageMetric_txx


#include "itkOptImageToImageMetric.h"
#include "itkImageRandomConstIteratorWithIndex.h"
#include "itkMersenneTwisterRandomVariateGenerator.h"

namespace itk
{

template <class TFixedImage, class TMovingImage> 
ImageToImageMetric<TFixedImage,TMovingImage>
::ImageToImageMetric()
{
  m_NumberOfFixedImageSamples = 50000;
  m_UseAllPixels = false;
  m_UseSequentialSampling = false;
  m_UseFixedImageIndexes = false;
  m_UseFixedImageSamplesIntensityThreshold = false;
  m_FixedImageSamplesIntensityThreshold = 0;
  m_ReseedIterator = false;
  m_RandomSeed = -1;

  m_TransformIsBSpline    = false;
  m_NumBSplineWeights     = 0;
  m_BSplineTransform      = NULL;

  m_Threader = MultiThreaderType::New();
  m_ThreaderParameter.metric = this;
  m_ThreaderNumberOfMovingImageSamples = NULL;
  m_WithinThreadPreProcess = false;
  m_WithinThreadPostProcess = false;

  m_FixedImage    = 0; // has to be provided by the user.
  m_FixedImageMask = 0;

  m_MovingImage   = 0; // has to be provided by the user.
  m_MovingImageMask = 0;
  m_NumberOfPixelsCounted = 0;

  m_Transform         = NULL; // has to be provided by the user.
  m_ThreaderTransform = NULL; // constructed at initialization.

  m_Interpolator  = 0; // has to be provided by the user.

  m_GradientImage = 0; // will receive the output of the filter;
  m_ComputeGradient = true; // metric computes gradient by default
  m_GradientImage = NULL; // computed at initialization

  m_InterpolatorIsBSpline = false;
  m_BSplineInterpolator = NULL;
  m_DerivativeCalculator = NULL;

  m_NumberOfThreads = m_Threader->GetNumberOfThreads();

  this->m_ThreaderBSplineTransformWeights = NULL;
  this->m_ThreaderBSplineTransformIndices = NULL;

  this->m_UseCachingOfBSplineWeights = true;

  /* if 100% backward compatible, we should include this...but...
  typename BSplineTransformType::Pointer transformer =
           BSplineTransformType::New();
  this->SetTransform (transformer);

  typename BSplineInterpolatorType::Pointer interpolator =
           BSplineInterpolatorType::New();
  this->SetInterpolator (interpolator);
  */
}

template <class TFixedImage, class TMovingImage> 
ImageToImageMetric<TFixedImage,TMovingImage>
::~ImageToImageMetric()
{

  if(m_ThreaderNumberOfMovingImageSamples != NULL)
    {
    delete [] m_ThreaderNumberOfMovingImageSamples;
    }
  m_ThreaderNumberOfMovingImageSamples = NULL;

  if(m_ThreaderTransform != NULL)
    {
    delete [] m_ThreaderTransform;
    }
  m_ThreaderTransform = NULL;

  if( this->m_ThreaderBSplineTransformWeights != NULL )
    {
    delete [] this->m_ThreaderBSplineTransformWeights;
    }
  this->m_ThreaderBSplineTransformWeights = NULL;
  
  if( this->m_ThreaderBSplineTransformIndices != NULL )
    {
    delete [] this->m_ThreaderBSplineTransformIndices;
    }
  this->m_ThreaderBSplineTransformIndices = NULL;
}

template <class TFixedImage, class TMovingImage> 
void
ImageToImageMetric<TFixedImage,TMovingImage>
::SetNumberOfThreads( unsigned int numberOfThreads )
{
  m_Threader->SetNumberOfThreads( numberOfThreads);
  m_NumberOfThreads = m_Threader->GetNumberOfThreads();
}

template <class TFixedImage, class TMovingImage> 
void
ImageToImageMetric<TFixedImage,TMovingImage>
::SetTransformParameters( const ParametersType & parameters ) const
{
  if( !m_Transform )
    {
    itkExceptionMacro(<<"Transform has not been assigned");
    }
  m_Transform->SetParameters( parameters );

  m_Parameters = parameters;
}

template <class TFixedImage, class TMovingImage> 
void
ImageToImageMetric<TFixedImage,TMovingImage>
::SetNumberOfFixedImageSamples( unsigned long numSamples )
{
  if( numSamples != m_NumberOfFixedImageSamples )
    {
    m_NumberOfFixedImageSamples = numSamples;
    if( m_NumberOfFixedImageSamples != this->m_FixedImageRegion.GetNumberOfPixels() )
      {
      this->SetUseAllPixels( false );
      }
    this->Modified();
    }
}

template <class TFixedImage, class TMovingImage> 
void
ImageToImageMetric<TFixedImage,TMovingImage>
::SetFixedImageIndexes( const FixedImageIndexContainer & indexes )
{
  this->SetUseFixedImageIndexes( true );
  m_NumberOfFixedImageSamples = indexes.size();
  m_FixedImageIndexes.resize( m_NumberOfFixedImageSamples );
  for(unsigned int i=0; i<m_NumberOfFixedImageSamples; i++)
    {
    m_FixedImageIndexes[i] = indexes[i];
    }
}

template <class TFixedImage, class TMovingImage> 
void
ImageToImageMetric<TFixedImage,TMovingImage>
::SetUseFixedImageIndexes( bool useIndexes )
{
  if( useIndexes != m_UseFixedImageIndexes )
    {
    m_UseFixedImageIndexes = useIndexes;
    if( m_UseFixedImageIndexes )
      {
      this->SetUseAllPixels( false );
      }
    else
      {
      this->Modified();
      }
    }
}

template <class TFixedImage, class TMovingImage> 
void
ImageToImageMetric<TFixedImage,TMovingImage>
::SetFixedImageSamplesIntensityThreshold( const FixedImagePixelType & thresh )
{
  if( thresh != m_FixedImageSamplesIntensityThreshold )
    {
    m_FixedImageSamplesIntensityThreshold = thresh;
    this->SetUseFixedImageSamplesIntensityThreshold( true );
    this->Modified();
    }
}

template <class TFixedImage, class TMovingImage> 
void
ImageToImageMetric<TFixedImage,TMovingImage>
::SetUseFixedImageSamplesIntensityThreshold( bool useThresh )
{
  if( useThresh != m_UseFixedImageSamplesIntensityThreshold )
    {
    m_UseFixedImageSamplesIntensityThreshold = useThresh;
    if( m_UseFixedImageSamplesIntensityThreshold )
      {
      this->SetUseAllPixels( false );
      }
    else
      {
      this->Modified();
      }
    }
}

template <class TFixedImage, class TMovingImage> 
void
ImageToImageMetric<TFixedImage,TMovingImage>
::SetFixedImageRegion( const FixedImageRegionType reg )
{
  if( reg != m_FixedImageRegion )
    {
    m_FixedImageRegion = reg;
    if( this->GetUseAllPixels() )
      {
      this->SetNumberOfFixedImageSamples( this->m_FixedImageRegion.GetNumberOfPixels() );
      }
    }
}

template <class TFixedImage, class TMovingImage> 
void
ImageToImageMetric<TFixedImage,TMovingImage>
::SetUseAllPixels( bool useAllPixels )
{
  if( useAllPixels != m_UseAllPixels )
    {
    m_UseAllPixels = useAllPixels;
    if( m_UseAllPixels )
      {
      this->SetUseFixedImageSamplesIntensityThreshold( false );
      this->SetNumberOfFixedImageSamples( this->m_FixedImageRegion.GetNumberOfPixels() );
      this->SetUseSequentialSampling( true );
      }
    else
      {
      this->SetUseSequentialSampling( false );
      this->Modified();
      }
    }
}


template <class TFixedImage, class TMovingImage> 
void
ImageToImageMetric<TFixedImage,TMovingImage>
::SetUseSequentialSampling( bool useSequential )
{
  if( useSequential != m_UseSequentialSampling )
    {
    m_UseSequentialSampling = useSequential;
    if( !m_UseSequentialSampling )
      {
      this->SetUseAllPixels( false );
      }
    else
      {
      this->Modified();
      }
    }
}

template <class TFixedImage, class TMovingImage> 
void
ImageToImageMetric<TFixedImage,TMovingImage>
::Initialize(void) throw ( ExceptionObject )
{

  if( !m_Transform )
    {
    itkExceptionMacro(<<"Transform is not present");
    }
  m_NumberOfParameters = m_Transform->GetNumberOfParameters();

  if( !m_Interpolator )
    {
    itkExceptionMacro(<<"Interpolator is not present");
    }

  if( !m_MovingImage )
    {
    itkExceptionMacro(<<"MovingImage is not present");
    }

  if( !m_FixedImage )
    {
    itkExceptionMacro(<<"FixedImage is not present");
    }

  if( m_FixedImageRegion.GetNumberOfPixels() == 0 )
    {
    itkExceptionMacro(<<"FixedImageRegion is empty");
    }

  // If the image is provided by a source, update the source.
  if( m_MovingImage->GetSource() )
    {
    m_MovingImage->GetSource()->Update();
    }

  // If the image is provided by a source, update the source.
  if( m_FixedImage->GetSource() )
    {
    m_FixedImage->GetSource()->Update();
    }

  // Make sure the FixedImageRegion is within the FixedImage buffered region
  if ( !m_FixedImageRegion.Crop( m_FixedImage->GetBufferedRegion() ) )
    {
    itkExceptionMacro(
      <<"FixedImageRegion does not overlap the fixed image buffered region" );
    }

  m_Interpolator->SetInputImage( m_MovingImage );
 
  if ( m_ComputeGradient )
    {
    ComputeGradient();
    }

  // If there are any observers on the metric, call them to give the
  // user code a chance to set parameters on the metric
  this->InvokeEvent( InitializeEvent() );

}


template <class TFixedImage, class TMovingImage> 
void
ImageToImageMetric<TFixedImage,TMovingImage>
::MultiThreadingInitialize(void) throw ( ExceptionObject )
{

  m_Threader->SetNumberOfThreads( m_NumberOfThreads );

  if(m_ThreaderNumberOfMovingImageSamples != NULL)
    {
    delete [] m_ThreaderNumberOfMovingImageSamples;
    }
  m_ThreaderNumberOfMovingImageSamples = new unsigned int[m_NumberOfThreads-1];

  // Allocate the array of transform clones to be used in every thread
  if(m_ThreaderTransform != NULL)
    {
    delete [] m_ThreaderTransform;
    }
  m_ThreaderTransform = new TransformPointer[m_NumberOfThreads-1];
  for( unsigned int ithread=0; ithread < m_NumberOfThreads-1; ++ithread)
    {
    // Create a copy of the main transform to be used in this thread.
    LightObject::Pointer anotherTransform = this->m_Transform->CreateAnother();
    // This static_cast should always work since the pointer was created by
    // CreateAnother() called from the transform itself.
    TransformType * transformCopy = static_cast< TransformType * >( anotherTransform.GetPointer() );
    transformCopy->SetFixedParameters( this->m_Transform->GetFixedParameters() );
    transformCopy->SetParameters( this->m_Transform->GetParameters() );
    this->m_ThreaderTransform[ithread] = transformCopy;
    }

  m_FixedImageSamples.resize( m_NumberOfFixedImageSamples );
  if( m_UseSequentialSampling )
    {
    // 
    // Take all the pixels within the fixed image region)
    // to create the sample points list.
    // 
    SampleFullFixedImageRegion( m_FixedImageSamples );
    }
  else
    {
    if( m_UseFixedImageIndexes )
      {
      //
      //  Use the list of indexes passed to the SetFixedImageIndexes
      //  member function .
      //
      SampleFixedImageIndexes( m_FixedImageSamples );
      }
    else
      {
      // 
      // Uniformly sample the fixed image (within the fixed image region)
      // to create the sample points list.
      // 
      SampleFixedImageRegion( m_FixedImageSamples );
      }
    }

  //  
  //  Check if the interpolator is of type BSplineInterpolateImageFunction.
  //  If so, we can make use of its EvaluateDerivatives method.
  //  Otherwise, we instantiate an external central difference
  //  derivative calculator.
  //  
  m_InterpolatorIsBSpline = true;

  BSplineInterpolatorType * testPtr = dynamic_cast<BSplineInterpolatorType *>(
    this->m_Interpolator.GetPointer() );
  if ( !testPtr )
    {
    m_InterpolatorIsBSpline = false;

    m_DerivativeCalculator = DerivativeFunctionType::New();

#ifdef ITK_USE_ORIENTED_IMAGE_DIRECTION
    m_DerivativeCalculator->UseImageDirectionOn();
#endif

    m_DerivativeCalculator->SetInputImage( this->m_MovingImage );

    m_BSplineInterpolator = NULL;
    itkDebugMacro( "Interpolator is not BSpline" );
    } 
  else
    {
    m_BSplineInterpolator = testPtr;
    m_BSplineInterpolator->SetNumberOfThreads( m_NumberOfThreads );

#ifdef ITK_USE_ORIENTED_IMAGE_DIRECTION
    m_BSplineInterpolator->UseImageDirectionOn();
#endif

    m_DerivativeCalculator = NULL;
    itkDebugMacro( "Interpolator is BSpline" );
    }

  //  
  //  Check if the transform is of type BSplineDeformableTransform.
  //  
  //  If so, several speed up features are implemented.
  //  [1] Precomputing the results of bulk transform for each sample point.
  //  [2] Precomputing the BSpline weights for each sample point,
  //      to be used later to directly compute the deformation vector
  //  [3] Precomputing the indices of the parameters within the 
  //      the support region of each sample point.
  //  
  m_TransformIsBSpline = true;

  BSplineTransformType * testPtr2 = dynamic_cast<BSplineTransformType *>(
    this->m_Transform.GetPointer() );
  if( !testPtr2 )
    {
    m_TransformIsBSpline = false;
    m_BSplineTransform = NULL;
    itkDebugMacro( "Transform is not BSplineDeformable" );
    }
  else
    {
    m_BSplineTransform = testPtr2;
    m_NumBSplineWeights = m_BSplineTransform->GetNumberOfWeights();
    itkDebugMacro( "Transform is BSplineDeformable" );
    }

  if( this->m_TransformIsBSpline )
    {
    // First, deallocate memory that may have been used from previous run of the Metric
    this->m_BSplineTransformWeightsArray.SetSize( 1, 1 );
    this->m_BSplineTransformIndicesArray.SetSize( 1, 1 );
    this->m_BSplinePreTransformPointsArray.resize( 1 );
    this->m_WithinBSplineSupportRegionArray.resize( 1 );
    this->m_BSplineTransformWeights.SetSize( 1 );
    this->m_BSplineTransformIndices.SetSize( 1 );

    if( this->m_ThreaderBSplineTransformWeights != NULL )
      {
      delete [] this->m_ThreaderBSplineTransformWeights;
      }
    this->m_ThreaderBSplineTransformWeights = NULL;
    
    if( this->m_ThreaderBSplineTransformIndices != NULL )
      {
      delete [] this->m_ThreaderBSplineTransformIndices;
      }
    this->m_ThreaderBSplineTransformIndices = NULL;

    if( this->m_UseCachingOfBSplineWeights )
      {
      m_BSplineTransformWeightsArray.SetSize( 
        m_NumberOfFixedImageSamples, m_NumBSplineWeights );
      m_BSplineTransformIndicesArray.SetSize( 
        m_NumberOfFixedImageSamples, m_NumBSplineWeights );
      m_BSplinePreTransformPointsArray.resize( m_NumberOfFixedImageSamples );
      m_WithinBSplineSupportRegionArray.resize( m_NumberOfFixedImageSamples );

      this->PreComputeTransformValues();
      }
    else
      {
      this->m_BSplineTransformWeights.SetSize( this->m_NumBSplineWeights );
      this->m_BSplineTransformIndices.SetSize( this->m_NumBSplineWeights );

      this->m_ThreaderBSplineTransformWeights = new BSplineTransformWeightsType[m_NumberOfThreads-1];
      this->m_ThreaderBSplineTransformIndices = new BSplineTransformIndexArrayType[m_NumberOfThreads-1];

      for( unsigned int ithread=0; ithread < m_NumberOfThreads-1; ++ithread)
        {
        this->m_ThreaderBSplineTransformWeights[ithread].SetSize( this->m_NumBSplineWeights );
        this->m_ThreaderBSplineTransformIndices[ithread].SetSize( this->m_NumBSplineWeights );
        }
      }

    for ( unsigned int j = 0; j < FixedImageDimension; j++ )
      {
      this->m_BSplineParametersOffset[j] = j * this->m_BSplineTransform->GetNumberOfParametersPerDimension();
      }
    }

}


template < class TFixedImage, class TMovingImage >
void
ImageToImageMetric<TFixedImage,TMovingImage>
::SampleFixedImageIndexes( FixedImageSampleContainer & samples ) const
{
  typename FixedImageSampleContainer::iterator iter;

  unsigned long len = m_FixedImageIndexes.size();
  if( len != m_NumberOfFixedImageSamples 
       || samples.size() != m_NumberOfFixedImageSamples )
    {
    throw ExceptionObject(__FILE__, __LINE__, 
       "Index list size does not match desired number of samples" );
    }

  iter=samples.begin();
  for(unsigned long i=0; i<len; i++)
    {
    // Get sampled index
    FixedImageIndexType index = m_FixedImageIndexes[i];
    // Translate index to point
    m_FixedImage->TransformIndexToPhysicalPoint( index, (*iter).point );

    // Get sampled fixed image value
    (*iter).value = m_FixedImage->GetPixel( index );
    (*iter).valueIndex = 0;

    ++iter;
    }
}

template < class TFixedImage, class TMovingImage >
void
ImageToImageMetric<TFixedImage,TMovingImage>
::SampleFixedImageRegion( FixedImageSampleContainer & samples ) const
{
  if( samples.size() != m_NumberOfFixedImageSamples )
    {
    throw ExceptionObject(__FILE__, __LINE__, 
       "Sample size does not match desired number of samples" );
    }

  // Set up a random interator within the user specified fixed image region.
  typedef ImageRandomConstIteratorWithIndex<FixedImageType> RandomIterator;
  RandomIterator randIter( m_FixedImage, GetFixedImageRegion() );

  typename FixedImageSampleContainer::iterator iter;
  typename FixedImageSampleContainer::const_iterator end=samples.end();

  if( m_FixedImageMask.IsNotNull()
      || m_UseFixedImageSamplesIntensityThreshold )
    {
    InputPointType inputPoint;

    iter=samples.begin();
    unsigned long int samplesFound = 0;
    randIter.SetNumberOfSamples( m_NumberOfFixedImageSamples * 1000 );
    randIter.GoToBegin();
    while( iter != end )
      {
      if( randIter.IsAtEnd() )
        {
        // Must be a small mask since after many random samples we don't
        // have enough to fill the desired number.   So, we will replicate
        // the samples we've found so far to fill-in the desired number
        // of samples
        unsigned long int count = 0;
        while( iter != end )
          {
          (*iter).point = samples[count].point;
          (*iter).value = samples[count].value;
          (*iter).valueIndex = 0;
          ++count;
          if(count >= samplesFound)
            {
            count = 0;
            }
          ++iter;
          }
        break;
        }
      
      // Get sampled index
      FixedImageIndexType index = randIter.GetIndex();
      // Check if the Index is inside the mask, translate index to point
      m_FixedImage->TransformIndexToPhysicalPoint( index, inputPoint );

      if( m_FixedImageMask.IsNotNull() )
        {
        double val;
        if( m_FixedImageMask->ValueAt( inputPoint, val ) )
          {
          if( val == 0 )
            {
            ++randIter; // jump to another random position
            continue;
            }
          }
        else
          {
          ++randIter; // jump to another random position
          continue;
          }
        }

      if( m_UseFixedImageSamplesIntensityThreshold &&
          randIter.Get() < m_FixedImageSamplesIntensityThreshold )
        {
        ++randIter;
        continue;
        }

      // Translate index to point
      (*iter).point = inputPoint;
      // Get sampled fixed image value
      (*iter).value = randIter.Get();
      (*iter).valueIndex = 0;

      ++samplesFound;
      ++randIter;
      ++iter;
      }
    }
  else
    {
    randIter.SetNumberOfSamples( m_NumberOfFixedImageSamples );
    randIter.GoToBegin();
    for( iter=samples.begin(); iter != end; ++iter )
      {
      // Get sampled index
      FixedImageIndexType index = randIter.GetIndex();
      // Translate index to point
      m_FixedImage->TransformIndexToPhysicalPoint( index,
                                                (*iter).point );
      // Get sampled fixed image value
      (*iter).value = randIter.Get();
      (*iter).valueIndex = 0;

      // Jump to random position
      ++randIter;
      }
    }
}

template < class TFixedImage, class TMovingImage >
void
ImageToImageMetric<TFixedImage,TMovingImage>
::SampleFullFixedImageRegion( FixedImageSampleContainer& samples ) const
{

  if( samples.size() != m_NumberOfFixedImageSamples )
    {
    throw ExceptionObject(__FILE__, __LINE__, 
       "Sample size does not match desired number of samples" );
    }

  // Set up a region interator within the user specified fixed image region.
  typedef ImageRegionConstIteratorWithIndex<FixedImageType> RegionIterator;
  RegionIterator regionIter( m_FixedImage, GetFixedImageRegion() );

  regionIter.GoToBegin();

  typename FixedImageSampleContainer::iterator iter;
  typename FixedImageSampleContainer::const_iterator end=samples.end();

  if( m_FixedImageMask.IsNotNull()
      || m_UseFixedImageSamplesIntensityThreshold )
    {
    InputPointType inputPoint;

    // repeat until we get enough samples to fill the array
    iter=samples.begin();
    while( iter != end )
      {
      // Get sampled index
      FixedImageIndexType index = regionIter.GetIndex();
      // Check if the Index is inside the mask, translate index to point
      m_FixedImage->TransformIndexToPhysicalPoint( index, inputPoint );
 
      if( m_FixedImageMask.IsNotNull() )
        {
        // If not inside the mask, ignore the point
        if( !m_FixedImageMask->IsInside( inputPoint ) )
          {
          ++regionIter; // jump to next pixel
          if( regionIter.IsAtEnd() )
            {
            regionIter.GoToBegin();
            }
          continue;
          }
        }

      if( m_UseFixedImageSamplesIntensityThreshold &&
          regionIter.Get() < m_FixedImageSamplesIntensityThreshold )
        {
        ++regionIter; // jump to next pixel
        if( regionIter.IsAtEnd() )
          {
          regionIter.GoToBegin();
          }
        continue;
        }
 
      // Translate index to point
      (*iter).point = inputPoint;
      // Get sampled fixed image value
      (*iter).value = regionIter.Get();
      (*iter).valueIndex = 0;
 
      ++regionIter;
      if( regionIter.IsAtEnd() )
        {
        regionIter.GoToBegin();
        }
      ++iter;
      }
    }
  else // not restricting sample throwing to a mask
    {
    for( iter=samples.begin(); iter != end; ++iter )
      {
      // Get sampled index
      FixedImageIndexType index = regionIter.GetIndex();

      // Translate index to point
      m_FixedImage->TransformIndexToPhysicalPoint( index,
                                               (*iter).point );
      // Get sampled fixed image value
      (*iter).value = regionIter.Get();
      (*iter).valueIndex = 0;

      ++regionIter;
      if( regionIter.IsAtEnd() )
        {
        regionIter.GoToBegin();
        }
      }
    }
}

template <class TFixedImage, class TMovingImage> 
void
ImageToImageMetric<TFixedImage,TMovingImage>
::ComputeGradient() 
{
  GradientImageFilterPointer gradientFilter = GradientImageFilterType::New();

  gradientFilter->SetInput( m_MovingImage );

  const typename MovingImageType::SpacingType & spacing = m_MovingImage
                                                          ->GetSpacing();
  double maximumSpacing=0.0;
  for(unsigned int i=0; i<MovingImageDimension; i++)
    {
    if( spacing[i] > maximumSpacing )
      {
      maximumSpacing = spacing[i];
      }
    }
  gradientFilter->SetSigma( maximumSpacing );
  gradientFilter->SetNormalizeAcrossScale( true );
  gradientFilter->SetNumberOfThreads( m_NumberOfThreads );
  
#ifdef ITK_USE_ORIENTED_IMAGE_DIRECTION
  gradientFilter->SetUseImageDirection( true );
#endif
  
  gradientFilter->Update();
  
  m_GradientImage = gradientFilter->GetOutput();
}

// Method to reinitialize the seed of the random number generator
template < class TFixedImage, class TMovingImage  > void
ImageToImageMetric<TFixedImage,TMovingImage>
::ReinitializeSeed()
{
  Statistics::MersenneTwisterRandomVariateGenerator::GetInstance()->SetSeed();
}

// Method to reinitialize the seed of the random number generator
template < class TFixedImage, class TMovingImage  > void
ImageToImageMetric<TFixedImage,TMovingImage>
::ReinitializeSeed(int seed)
{
  Statistics::MersenneTwisterRandomVariateGenerator::GetInstance()->SetSeed(
                                                                         seed);
}


template < class TFixedImage, class TMovingImage >
void
ImageToImageMetric<TFixedImage,TMovingImage>
::PreComputeTransformValues()
{
  // Note: This code is specific to the b-spline deformable transform.

  // Unfortunately, the BSplineDeformableTransform stores a 
  // pointer to parameters passed to SetParameters(). Since
  // we're creating a dummy set of parameters below on the
  // stack, this can cause a crash if the transform's 
  // parameters are not later reset with a more properly 
  // scoped set of parameters. In addition, we're overwriting
  // any previously set parameters. In order to be kinder,
  // we'll save a pointer to the current set of parameters 
  // and restore them after we're done.

  // Note the address operator.
  // const TransformParametersType* previousParameters = & m_Transform->GetParameters();

  // Create all zero dummy transform parameters
  ParametersType dummyParameters( m_NumberOfParameters );
  dummyParameters.Fill( 0.0 );
  m_Transform->SetParameters( dummyParameters );

  // Cycle through each sampled fixed image point
  BSplineTransformWeightsType weights( m_NumBSplineWeights );
  BSplineTransformIndexArrayType indices( m_NumBSplineWeights );
  bool valid;
  MovingImagePointType mappedPoint;

  // Declare iterators for iteration over the sample container
  typename FixedImageSampleContainer::const_iterator fiter;
  typename FixedImageSampleContainer::const_iterator fend = 
                                                     m_FixedImageSamples.end();
  unsigned long counter = 0;

  for( fiter = m_FixedImageSamples.begin(); fiter != fend; ++fiter, counter++ )
    {
    m_BSplineTransform->TransformPoint( m_FixedImageSamples[counter].point,
                                        mappedPoint, weights, indices, valid );

    for( unsigned long k = 0; k < m_NumBSplineWeights; k++ )
      {
      m_BSplineTransformWeightsArray[counter][k] = weights[k];
      m_BSplineTransformIndicesArray[counter][k] = indices[k];
      }

    m_BSplinePreTransformPointsArray[counter]      = mappedPoint;
    m_WithinBSplineSupportRegionArray[counter]     = valid;
    }

  // Restore the previous parameters.
  // m_Transform->SetParameters( *previousParameters );
}


template < class TFixedImage, class TMovingImage >
void
ImageToImageMetric<TFixedImage,TMovingImage>
::TransformPoint( unsigned int sampleNumber, 
                  MovingImagePointType & mappedPoint,
                  bool & sampleOk,
                  double & movingImageValue,
                  unsigned int threadID ) const
{
  sampleOk = true;
  TransformType * transform;
  
  if( threadID > 0 )
    {
    transform = this->m_ThreaderTransform[threadID-1];
    }
  else
    {
    transform = this->m_Transform;
    }

  if ( !m_TransformIsBSpline )
    {
    // Use generic transform to compute mapped position
    mappedPoint = transform->TransformPoint( m_FixedImageSamples[sampleNumber].point );
    sampleOk = true;
    }
  else
    {
    if( this->m_UseCachingOfBSplineWeights )
      {
      sampleOk = m_WithinBSplineSupportRegionArray[sampleNumber];

      if(sampleOk)
        {
        // If the transform is BSplineDeformable, we can use the precomputed
        // weights and indices to obtained the mapped position
        const WeightsValueType * weights = 
                                     m_BSplineTransformWeightsArray[sampleNumber];
        const IndexValueType   * indices = 
                                     m_BSplineTransformIndicesArray[sampleNumber];

        for( unsigned int j = 0; j < FixedImageDimension; j++ )
          {
          mappedPoint[j] = m_BSplinePreTransformPointsArray[sampleNumber][j];
          }
    
        for ( unsigned int k = 0; k < m_NumBSplineWeights; k++ )
          {
          for ( unsigned int j = 0; j < FixedImageDimension; j++ )
            {
            mappedPoint[j] += weights[k] * m_Parameters[ indices[k] 
                                                 + m_BSplineParametersOffset[j] ];
            }
          }
        }
      }
    else
      {
      BSplineTransformWeightsType    * weightsHelper;
      BSplineTransformIndexArrayType * indicesHelper;

      if( threadID > 0 )
        {
        weightsHelper = &(this->m_ThreaderBSplineTransformWeights[threadID-1]);
        indicesHelper = &(this->m_ThreaderBSplineTransformIndices[threadID-1]);
        }
      else
        {
        weightsHelper = &(this->m_BSplineTransformWeights);
        indicesHelper = &(this->m_BSplineTransformIndices);
        }

      // If not caching values, we invoke the Transform to recompute the
      // mapping of the point.
      this->m_BSplineTransform->TransformPoint( 
        this->m_FixedImageSamples[sampleNumber].point,
        mappedPoint, *weightsHelper, *indicesHelper, sampleOk);
      }
    }
  
  if(sampleOk)
    {
    // If user provided a mask over the Moving image
    if ( m_MovingImageMask )
      {
      // Check if mapped point is within the support region of the moving image
      // mask
      sampleOk = sampleOk && m_MovingImageMask->IsInside( mappedPoint );
      }

  
    if( m_InterpolatorIsBSpline )
      {
      // Check if mapped point inside image buffer
      sampleOk = sampleOk && m_BSplineInterpolator->IsInsideBuffer( mappedPoint );
      if( sampleOk )
        {
        movingImageValue = m_BSplineInterpolator->Evaluate( mappedPoint, threadID );
        }
      }
    else
      {
      // Check if mapped point inside image buffer
      sampleOk = sampleOk && m_Interpolator->IsInsideBuffer( mappedPoint );
      if( sampleOk )
        {
        movingImageValue = m_Interpolator->Evaluate( mappedPoint );
        }
      }
    }
}


template < class TFixedImage, class TMovingImage >
void
ImageToImageMetric<TFixedImage,TMovingImage>
::TransformPointWithDerivatives( unsigned int sampleNumber, 
                  MovingImagePointType& mappedPoint,
                  bool& sampleOk,
                  double& movingImageValue,
                  ImageDerivativesType & movingImageGradient,
                  unsigned int threadID ) const
{
  TransformType * transform;
  sampleOk = true;
  
  if( threadID > 0 )
    {
    transform = this->m_ThreaderTransform[threadID-1];
    }
  else
    {
    transform = this->m_Transform;
    }

  if ( !m_TransformIsBSpline )
    {
    // Use generic transform to compute mapped position
    mappedPoint = transform->TransformPoint( m_FixedImageSamples[sampleNumber].point );
    sampleOk = true;
    }
  else
    {
    if( this->m_UseCachingOfBSplineWeights )
      {
      sampleOk = m_WithinBSplineSupportRegionArray[sampleNumber];

      if(sampleOk)
        {
        // If the transform is BSplineDeformable, we can use the precomputed
        // weights and indices to obtained the mapped position
        const WeightsValueType * weights = 
                                     m_BSplineTransformWeightsArray[sampleNumber];
        const IndexValueType   * indices = 
                                     m_BSplineTransformIndicesArray[sampleNumber];

        for( unsigned int j = 0; j < FixedImageDimension; j++ )
          {
          mappedPoint[j] = m_BSplinePreTransformPointsArray[sampleNumber][j];
          }
    
        for ( unsigned int k = 0; k < m_NumBSplineWeights; k++ )
          {
          for ( unsigned int j = 0; j < FixedImageDimension; j++ )
            {
            mappedPoint[j] += weights[k] * m_Parameters[ indices[k] 
                                                 + m_BSplineParametersOffset[j] ];
            }
          }
        }
      }
    else
      {
      BSplineTransformWeightsType    * weightsHelper;
      BSplineTransformIndexArrayType * indicesHelper;

      if( threadID > 0 )
        {
        weightsHelper = &(this->m_ThreaderBSplineTransformWeights[threadID-1]);
        indicesHelper = &(this->m_ThreaderBSplineTransformIndices[threadID-1]);
        }
      else
        {
        weightsHelper = &(this->m_BSplineTransformWeights);
        indicesHelper = &(this->m_BSplineTransformIndices);
        }

      // If not caching values, we invoke the Transform to recompute the
      // mapping of the point.
      this->m_BSplineTransform->TransformPoint( 
        this->m_FixedImageSamples[sampleNumber].point,
        mappedPoint, *weightsHelper, *indicesHelper, sampleOk);
      }
    }
  
  if(sampleOk)
    {
    // If user provided a mask over the Moving image
    if ( m_MovingImageMask )
      {
      // Check if mapped point is within the support region of the moving image
      // mask
      sampleOk = sampleOk && m_MovingImageMask->IsInside( mappedPoint );
      }
 
  
    if( m_InterpolatorIsBSpline )
      {
      // Check if mapped point inside image buffer
      sampleOk = sampleOk && m_BSplineInterpolator->IsInsideBuffer( mappedPoint );
      if( sampleOk )
        {
        this->m_BSplineInterpolator->EvaluateValueAndDerivative(mappedPoint,
                                                          movingImageValue,
                                                          movingImageGradient,
                                                          threadID);
        }
      }
    else
      {
      // Check if mapped point inside image buffer
      sampleOk = sampleOk && m_Interpolator->IsInsideBuffer( mappedPoint );
      if( sampleOk )
        {
        this->ComputeImageDerivatives( mappedPoint, movingImageGradient, threadID );
        movingImageValue = this->m_Interpolator->Evaluate( mappedPoint );
        }
      }
    }
}

template < class TFixedImage, class TMovingImage >
void
ImageToImageMetric<TFixedImage,TMovingImage>
::ComputeImageDerivatives( const MovingImagePointType & mappedPoint, 
                           ImageDerivativesType & gradient,
                           unsigned int threadID) const
{
  
  if( m_InterpolatorIsBSpline )
    {
    // Computed moving image gradient using derivative BSpline kernel.
    gradient = m_BSplineInterpolator->EvaluateDerivative( mappedPoint,
                                                          threadID );
    }
  else
    {
    if ( m_ComputeGradient )
      {
      ContinuousIndex<double, MovingImageDimension> tempIndex;
      m_MovingImage->TransformPhysicalPointToContinuousIndex( mappedPoint,
                                                              tempIndex );
      MovingImageIndexType mappedIndex;
      mappedIndex.CopyWithRound( tempIndex );
      gradient = m_GradientImage->GetPixel( mappedIndex );
      }
    else
      {
      // if not using the gradient image
      gradient = m_DerivativeCalculator->Evaluate( mappedPoint );
      }
    }

}

template < class TFixedImage, class TMovingImage  >
void
ImageToImageMetric<TFixedImage,TMovingImage>
::GetValueMultiThreadedPreProcessInitiate( void ) const
{
  this->SynchronizeTransforms();

  m_Threader->SetSingleMethod(GetValueMultiThreadedPreProcess,
                              (void *)(&m_ThreaderParameter));
  m_Threader->SingleMethodExecute();
}

template < class TFixedImage, class TMovingImage  >
void
ImageToImageMetric<TFixedImage,TMovingImage>
::GetValueMultiThreadedInitiate( void ) const
{

  this->SynchronizeTransforms();

  m_Threader->SetSingleMethod(GetValueMultiThreaded,
                              const_cast<void *>(static_cast<const void *>(&m_ThreaderParameter)));
  m_Threader->SingleMethodExecute();

  for( unsigned int threadID = 0; threadID<m_NumberOfThreads-1; threadID++ )
    {
    this->m_NumberOfPixelsCounted += m_ThreaderNumberOfMovingImageSamples[threadID];
    }
}

template < class TFixedImage, class TMovingImage  >
void
ImageToImageMetric<TFixedImage,TMovingImage>
::GetValueMultiThreadedPostProcessInitiate( void ) const
{
  m_Threader->SetSingleMethod(GetValueMultiThreadedPostProcess,
                              const_cast<void *>(static_cast<const void *>(&m_ThreaderParameter)));
  m_Threader->SingleMethodExecute();
}

template < class TFixedImage, class TMovingImage  >
ITK_THREAD_RETURN_TYPE
ImageToImageMetric<TFixedImage,TMovingImage>
::GetValueMultiThreadedPreProcess( void * arg ) 
{
  int threadID;
  MultiThreaderParameterType * mtParam;

  threadID = ((MultiThreaderType::ThreadInfoStruct *)(arg))->ThreadID;

  mtParam = (MultiThreaderParameterType *)
            (((MultiThreaderType::ThreadInfoStruct *)(arg))->UserData);

  mtParam->metric->GetValueThreadPreProcess(threadID, false);

  return ITK_THREAD_RETURN_VALUE;
}


template < class TFixedImage, class TMovingImage  >
ITK_THREAD_RETURN_TYPE
ImageToImageMetric<TFixedImage,TMovingImage>
::GetValueMultiThreaded( void * arg ) 
{
  int threadID;
  MultiThreaderParameterType * mtParam;

  threadID = ((MultiThreaderType::ThreadInfoStruct *)(arg))->ThreadID;

  mtParam = (MultiThreaderParameterType *)
            (((MultiThreaderType::ThreadInfoStruct *)(arg))->UserData);

  mtParam->metric->GetValueThread(threadID);

  return ITK_THREAD_RETURN_VALUE;
}
template < class TFixedImage, class TMovingImage  >
ITK_THREAD_RETURN_TYPE
ImageToImageMetric<TFixedImage,TMovingImage>
::GetValueMultiThreadedPostProcess( void * arg ) 
{
  int threadID;
  MultiThreaderParameterType * mtParam;

  threadID = ((MultiThreaderType::ThreadInfoStruct *)(arg))->ThreadID;

  mtParam = (MultiThreaderParameterType *)
            (((MultiThreaderType::ThreadInfoStruct *)(arg))->UserData);

  mtParam->metric->GetValueThreadPostProcess(threadID, false);

  return ITK_THREAD_RETURN_VALUE;
}


template < class TFixedImage, class TMovingImage  >
void
ImageToImageMetric<TFixedImage,TMovingImage>
::GetValueThread( unsigned int threadID ) const
{
  // Figure out how many samples to process
  int chunkSize = m_NumberOfFixedImageSamples / m_NumberOfThreads;

  // Skip to this thread's samples to process
  unsigned int fixedImageSample = threadID * chunkSize;

  if(threadID == m_NumberOfThreads - 1)
    {
    chunkSize = m_NumberOfFixedImageSamples 
                              - ((m_NumberOfThreads-1) 
                                 * chunkSize);
    }

  int numSamples = 0;

  if(m_WithinThreadPreProcess)
    {
    this->GetValueThreadPreProcess(threadID, true);
    }

  // Process the samples
  MovingImagePointType mappedPoint;
  bool sampleOk;
  double movingImageValue;
  for( int count=0; count < chunkSize; ++count, ++fixedImageSample )
    {
    // Get moving image value
    this->TransformPoint( fixedImageSample, mappedPoint, sampleOk, movingImageValue,
                          threadID );

    if( sampleOk )
      {
      // CALL USER FUNCTION
      if(GetValueThreadProcessSample(threadID, fixedImageSample,
                                     mappedPoint, movingImageValue))
        {
        ++numSamples;
        }
      }
    }

  if(threadID > 0)
    {
    m_ThreaderNumberOfMovingImageSamples[threadID-1] = numSamples;
    }
  else
    {
    m_NumberOfPixelsCounted = numSamples;
    }

  if(m_WithinThreadPostProcess)
    {
    this->GetValueThreadPostProcess(threadID, true);
    }
}

template < class TFixedImage, class TMovingImage  >
void
ImageToImageMetric<TFixedImage,TMovingImage>
::GetValueAndDerivativeMultiThreadedPreProcessInitiate( void ) const
{
  this->SynchronizeTransforms();

  m_Threader->SetSingleMethod(GetValueAndDerivativeMultiThreadedPreProcess,
                              (void *)(&m_ThreaderParameter));
  m_Threader->SingleMethodExecute();
}

template < class TFixedImage, class TMovingImage  >
void
ImageToImageMetric<TFixedImage,TMovingImage>
::GetValueAndDerivativeMultiThreadedInitiate( void ) const
{
  this->SynchronizeTransforms();

  m_Threader->SetSingleMethod(GetValueAndDerivativeMultiThreaded,
                              const_cast<void *>(static_cast<const void *>(&m_ThreaderParameter)));
  m_Threader->SingleMethodExecute();

  for( unsigned int threadID = 0; threadID<m_NumberOfThreads-1; threadID++ )
    {
    this->m_NumberOfPixelsCounted += m_ThreaderNumberOfMovingImageSamples[threadID];
    }
}

template < class TFixedImage, class TMovingImage  >
void
ImageToImageMetric<TFixedImage,TMovingImage>
::GetValueAndDerivativeMultiThreadedPostProcessInitiate( void ) const
{
  m_Threader->SetSingleMethod(GetValueAndDerivativeMultiThreadedPostProcess,
                              (void *)(&m_ThreaderParameter));
  m_Threader->SingleMethodExecute();
}

template < class TFixedImage, class TMovingImage  >
ITK_THREAD_RETURN_TYPE
ImageToImageMetric<TFixedImage,TMovingImage>
::GetValueAndDerivativeMultiThreadedPreProcess( void * arg ) 
{
  int threadID;
  MultiThreaderParameterType * mtParam;

  threadID = ((MultiThreaderType::ThreadInfoStruct *)(arg))->ThreadID;

  mtParam = (MultiThreaderParameterType *)
            (((MultiThreaderType::ThreadInfoStruct *)(arg))->UserData);

  mtParam->metric->GetValueAndDerivativeThreadPreProcess(threadID, false);

  return ITK_THREAD_RETURN_VALUE;
}

template < class TFixedImage, class TMovingImage  >
ITK_THREAD_RETURN_TYPE
ImageToImageMetric<TFixedImage,TMovingImage>
::GetValueAndDerivativeMultiThreaded( void * arg ) 
{
  int threadID;
  MultiThreaderParameterType * mtParam;

  threadID = ((MultiThreaderType::ThreadInfoStruct *)(arg))->ThreadID;

  mtParam = (MultiThreaderParameterType *)
            (((MultiThreaderType::ThreadInfoStruct *)(arg))->UserData);

  mtParam->metric->GetValueAndDerivativeThread(threadID);

  return ITK_THREAD_RETURN_VALUE;
}

template < class TFixedImage, class TMovingImage  >
ITK_THREAD_RETURN_TYPE
ImageToImageMetric<TFixedImage,TMovingImage>
::GetValueAndDerivativeMultiThreadedPostProcess( void * arg ) 
{
  int threadID;
  MultiThreaderParameterType * mtParam;

  threadID = ((MultiThreaderType::ThreadInfoStruct *)(arg))->ThreadID;

  mtParam = (MultiThreaderParameterType *)
            (((MultiThreaderType::ThreadInfoStruct *)(arg))->UserData);

  mtParam->metric->GetValueAndDerivativeThreadPostProcess(threadID, false);

  return ITK_THREAD_RETURN_VALUE;
}

template < class TFixedImage, class TMovingImage  >
void
ImageToImageMetric<TFixedImage,TMovingImage>
::GetValueAndDerivativeThread( unsigned int threadID ) const
{
  // Figure out how many samples to process
  int chunkSize = m_NumberOfFixedImageSamples / m_NumberOfThreads;

  // Skip to this thread's samples to process
  unsigned int fixedImageSample = threadID * chunkSize;

  if(threadID == m_NumberOfThreads - 1)
    {
    chunkSize = m_NumberOfFixedImageSamples 
                              - ((m_NumberOfThreads-1) 
                                 * chunkSize);
    }

  int numSamples = 0;

  if(m_WithinThreadPreProcess)
    {
    this->GetValueAndDerivativeThreadPreProcess(threadID, true);
    }

  // Process the samples
  MovingImagePointType mappedPoint;
  bool sampleOk;
  double movingImageValue;
  ImageDerivativesType movingImageGradientValue;
  for( int count=0; count < chunkSize; ++count, ++fixedImageSample )
    {
    // Get moving image value
    TransformPointWithDerivatives( fixedImageSample, mappedPoint, sampleOk,
                                   movingImageValue, movingImageGradientValue,
                                   threadID );

    if( sampleOk )
      {
      // CALL USER FUNCTION
      if( this->GetValueAndDerivativeThreadProcessSample( threadID,
                                                          fixedImageSample,
                                                          mappedPoint,
                                                          movingImageValue,
                                                          movingImageGradientValue ))
        {
        ++numSamples;
        }
      }
    }

  if(threadID > 0)
    {
    m_ThreaderNumberOfMovingImageSamples[threadID-1] = numSamples;
    }
  else
    {
    m_NumberOfPixelsCounted = numSamples;
    }

  if(m_WithinThreadPostProcess)
    {
    this->GetValueAndDerivativeThreadPostProcess(threadID, true);
    }
}


template <class TFixedImage, class TMovingImage> 
void
ImageToImageMetric<TFixedImage,TMovingImage>
::PrintSelf(std::ostream& os, Indent indent) const
{
  Superclass::PrintSelf( os, indent );

  os << indent << "NumberOfFixedImageSamples: ";
  os << m_NumberOfFixedImageSamples << std::endl;

  os << indent << "FixedImageSamplesIntensityThreshold: "
     << static_cast<typename NumericTraits<FixedImagePixelType>::PrintType>(m_FixedImageSamplesIntensityThreshold)
     << std::endl;

  os << indent << "UseFixedImageSamplesIntensityThreshold: ";
  os << m_UseFixedImageSamplesIntensityThreshold << std::endl;

  if( m_UseFixedImageIndexes )
    {
    os << indent << "Use Fixed Image Indexes: True" << std::endl;
    os << indent << "Number of Fixed Image Indexes = "
                 << m_FixedImageIndexes.size() << std::endl;
    }
  else
    {
    os << indent << "Use Fixed Image Indexes: False" << std::endl;
    }

  if( m_UseSequentialSampling )
    {
    os << indent << "Use Sequential Sampling: True" << std::endl;
    }
  else
    {
    os << indent << "Use Sequential Sampling: False" << std::endl;
    }

  os << indent << "UseAllPixels: ";
  os << m_UseAllPixels << std::endl;

  os << indent << "Threader: " << m_Threader << std::endl;
  os << indent << "Number of Threads: " << m_NumberOfThreads << std::endl;
  os << indent << "ThreaderParameter: " << std::endl;
  os << indent << "ThreaderNumberOfMovingImageSamples: " << std::endl;
  if( m_ThreaderNumberOfMovingImageSamples )
    {
    for(unsigned int i=0; i<m_NumberOfThreads-1; i++)
      {
      os << "  Thread[" << i << "]= " << (unsigned int)m_ThreaderNumberOfMovingImageSamples[i] << std::endl;
      }
    }

  os << indent << "ComputeGradient: "
     << static_cast<typename NumericTraits<bool>::PrintType>(m_ComputeGradient)
     << std::endl;
  os << indent << "Moving Image: " << m_MovingImage.GetPointer()  << std::endl;
  os << indent << "Fixed  Image: " << m_FixedImage.GetPointer()   << std::endl;
  os << indent << "Gradient Image: " << m_GradientImage.GetPointer() 
     << std::endl;
  os << indent << "Transform:    " << m_Transform.GetPointer()    << std::endl;
  os << indent << "Interpolator: " << m_Interpolator.GetPointer() << std::endl;
  os << indent << "FixedImageRegion: " << m_FixedImageRegion << std::endl;
  os << indent << "Moving Image Mask: " << m_MovingImageMask.GetPointer() 
     << std::endl;
  os << indent << "Fixed Image Mask: " << m_FixedImageMask.GetPointer() 
     << std::endl;
  os << indent << "Number of Moving Image Samples: " << m_NumberOfPixelsCounted 
     << std::endl;

  os << indent << "UseCachingOfBSplineWeights: ";
  os << this->m_UseCachingOfBSplineWeights << std::endl;
}

template <class TFixedImage, class TMovingImage> 
void
ImageToImageMetric<TFixedImage,TMovingImage>
::SynchronizeTransforms() const
{
  for( unsigned int threadID = 0; threadID<m_NumberOfThreads-1; threadID++ )
    {
    this->m_ThreaderTransform[threadID]->SetFixedParameters( this->m_Transform->GetFixedParameters() );
    this->m_ThreaderTransform[threadID]->SetParameters( this->m_Transform->GetParameters() );
    }
}

} // end namespace itk

#endif