text: /*=========================================================================
   *
   *  Copyright Insight Software Consortium
   *
   *  Licensed under the Apache License, Version 2.0 (the "License");
   *  you may not use this file except in compliance with the License.
   *  You may obtain a copy of the License at
   *
   *         http://www.apache.org/licenses/LICENSE-2.0.txt
   *
   *  Unless required by applicable law or agreed to in writing, software
   *  distributed under the License is distributed on an "AS IS" BASIS,
   *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
   *  See the License for the specific language governing permissions and
   *  limitations under the License.
   *
   *=========================================================================*/
  #ifndef itkImageToImageMetricv4_hxx
  #define itkImageToImageMetricv4_hxx
  
  #include "itkImageToImageMetricv4.h"
  #include "itkPixelTraits.h"
  #include "itkDisplacementFieldTransform.h"
  #include "itkCompositeTransform.h"
  #include "itkLinearInterpolateImageFunction.h"
  #include "itkIdentityTransform.h"
  
  namespace itk
  {
  
  template<typename TFixedImage,typename TMovingImage,typename TVirtualImage, typename TInternalComputationValueType, typename TMetricTraits>
  ImageToImageMetricv4<TFixedImage, TMovingImage, TVirtualImage, TInternalComputationValueType, TMetricTraits>
  ::ImageToImageMetricv4()
  {
    /* Interpolators. Default to linear. */
    typedef LinearInterpolateImageFunction< FixedImageType,
                                            CoordinateRepresentationType >
                                                    FixedLinearInterpolatorType;
    typedef LinearInterpolateImageFunction< MovingImageType,
                                            CoordinateRepresentationType >
                                                    MovingLinearInterpolatorType;
    this->m_FixedInterpolator  = FixedLinearInterpolatorType::New();
    this->m_MovingInterpolator = MovingLinearInterpolatorType::New();
  
    /* Setup default gradient filter. It gets initialized with default
     * parameters during Initialize. */
    this->m_DefaultFixedImageGradientFilter  = DefaultFixedImageGradientFilter::New();
    this->m_DefaultMovingImageGradientFilter = DefaultMovingImageGradientFilter::New();
    this->m_FixedImageGradientFilter         = this->m_DefaultFixedImageGradientFilter;
    this->m_MovingImageGradientFilter        = this->m_DefaultMovingImageGradientFilter;
  
    /* Interpolators for image gradient filters */
    this->m_FixedImageGradientInterpolator  = FixedImageGradientInterpolatorType::New();
    this->m_MovingImageGradientInterpolator = MovingImageGradientInterpolatorType::New();
  
    /* Setup default gradient image function */
    this->m_DefaultFixedImageGradientCalculator = DefaultFixedImageGradientCalculator::New();
    this->m_DefaultFixedImageGradientCalculator->UseImageDirectionOn();
    this->m_FixedImageGradientCalculator = this->m_DefaultFixedImageGradientCalculator;
  
    this->m_DefaultMovingImageGradientCalculator = DefaultMovingImageGradientCalculator::New();
    this->m_DefaultMovingImageGradientCalculator->UseImageDirectionOn();
    this->m_MovingImageGradientCalculator = this->m_DefaultMovingImageGradientCalculator;
  
    /* Setup default options assuming dense-sampling */
    this->m_UseFixedImageGradientFilter  = true;
    this->m_UseMovingImageGradientFilter = true;
    this->m_UseFixedSampledPointSet      = false;
  
    this->m_FloatingPointCorrectionResolution = 1e6;
    this->m_UseFloatingPointCorrection = false;
  
    this->m_HaveMadeGetValueWarning = false;
    this->m_NumberOfSkippedFixedSampledPoints = 0;
  
    this->m_Value = NumericTraits<MeasureType>::max();
    this->m_DerivativeResult = ITK_NULLPTR;
    this->m_ComputeDerivative = false;
  }
  
  template<typename TFixedImage,typename TMovingImage,typename TVirtualImage, typename TInternalComputationValueType, typename TMetricTraits>
  ImageToImageMetricv4<TFixedImage, TMovingImage, TVirtualImage, TInternalComputationValueType, TMetricTraits>
  ::~ImageToImageMetricv4()
  {
  }
  
  template<typename TFixedImage,typename TMovingImage,typename TVirtualImage, typename TInternalComputationValueType, typename TMetricTraits>
  void
  ImageToImageMetricv4<TFixedImage, TMovingImage, TVirtualImage, TInternalComputationValueType, TMetricTraits>
  ::Initialize()
  {
    itkDebugMacro("Initialize entered");
  
    /* Verify things are connected */
    if ( this->m_FixedImage.IsNull() )
      {
      itkExceptionMacro(<< "FixedImage is not present");
      }
    if ( this->m_MovingImage.IsNull() )
      {
      itkExceptionMacro(<< "MovingImage is not present");
      }
    if ( this->m_FixedTransform.IsNull() )
      {
      itkExceptionMacro(<< "FixedTransform is not present");
      }
    if ( this->m_MovingTransform.IsNull() )
      {
      itkExceptionMacro(<< "MovingTransform is not present");
      }
  
    // If the image is provided by a source, update the source.
    if ( this->m_MovingImage->GetSource() )
      {
      this->m_MovingImage->GetSource()->Update();
      }
  
    // If the image is provided by a source, update the source.
    if ( this->m_FixedImage->GetSource() )
      {
      this->m_FixedImage->GetSource()->Update();
      }
  
    /* If a virtual image has not been set or created,
     * create one from fixed image settings */
    if( ! this->m_UserHasSetVirtualDomain )
      {
      /* Instantiate a virtual image, but do not call Allocate to allocate
       * the data, to save memory. We don't need data. We'll simply be iterating
       * over the image to get indecies and transform to points.
       * Note that it will be safer to have a dedicated VirtualImage class
       * that prevents accidental access of data. */
      /* Just copy information from fixed image */
      VirtualImagePointer image = VirtualImageType::New();
      image->CopyInformation( this->m_FixedImage );
      /* CopyInformation does not copy buffered region */
      image->SetBufferedRegion( this->m_FixedImage->GetBufferedRegion() );
      image->SetRequestedRegion( this->m_FixedImage->GetRequestedRegion() );
      this->SetVirtualDomainFromImage( image );
      }
  
    /*
     * Superclass Initialize.
     * Requires the above actions to already have been taken.
     */
    Superclass::Initialize();
  
    /* Map the fixed samples into the virtual domain and store in
     * a searpate point set. */
    if( this->m_UseFixedSampledPointSet )
      {
      this->MapFixedSampledPointSetToVirtual();
      }
  
    /* Inititialize interpolators. */
    itkDebugMacro("Initialize Interpolators");
    this->m_FixedInterpolator->SetInputImage( this->m_FixedImage );
    this->m_MovingInterpolator->SetInputImage( this->m_MovingImage );
  
    /* Setup for image gradient calculations. */
    if( ! this->m_UseFixedImageGradientFilter )
      {
      itkDebugMacro("Initialize FixedImageGradientCalculator");
      this->m_FixedImageGradientImage = ITK_NULLPTR;
      this->m_FixedImageGradientCalculator->SetInputImage(this->m_FixedImage);
      }
    if( ! this->m_UseMovingImageGradientFilter )
      {
      itkDebugMacro("Initialize MovingImageGradientCalculator");
      this->m_MovingImageGradientImage = ITK_NULLPTR;
      this->m_MovingImageGradientCalculator->SetInputImage(this->m_MovingImage);
      }
  
    /* Initialize default gradient image filters. */
    itkDebugMacro("InitializeDefaultFixedImageGradientFilter");
    this->InitializeDefaultFixedImageGradientFilter();
    itkDebugMacro("InitializeDefaultMovingImageGradientFilter");
    this->InitializeDefaultMovingImageGradientFilter();
  
    /* If user set to use a pre-calculated fixed gradient image,
     * and the metric is set to use fixed image gradients,
     * then we need to calculate the gradient image.
     * We only need to compute once. */
    if ( this->GetGradientSourceIncludesFixed() && this->m_UseFixedImageGradientFilter )
      {
      itkDebugMacro("Initialize: ComputeFixedImageGradientFilterImage");
      this->ComputeFixedImageGradientFilterImage();
      }
  
    /* Compute gradient image for moving image. */
    if( this->GetGradientSourceIncludesMoving() && this->m_UseMovingImageGradientFilter )
      {
      itkDebugMacro("Initialize: ComputeMovingImageGradientFilterImage");
      this->ComputeMovingImageGradientFilterImage();
      }
  }
  
  template<typename TFixedImage,typename TMovingImage,typename TVirtualImage, typename TInternalComputationValueType, typename TMetricTraits>
  typename ImageToImageMetricv4<TFixedImage, TMovingImage, TVirtualImage, TInternalComputationValueType, TMetricTraits>::MeasureType
  ImageToImageMetricv4<TFixedImage, TMovingImage, TVirtualImage, TInternalComputationValueType, TMetricTraits>
  ::GetValue() const
  {
    this->m_ComputeDerivative = false;
  
    DerivativeType derivative;
    this->m_DerivativeResult = &derivative;
    this->InitializeForIteration();
  
    // Do the threaded processing using the appropriate
    // GetValueAndDerivativeThreader. Results get written to
    // member vars.
    this->GetValueAndDerivativeExecute();
  
    return this->m_Value;
  }
  
  template<typename TFixedImage,typename TMovingImage,typename TVirtualImage, typename TInternalComputationValueType, typename TMetricTraits>
  void
  ImageToImageMetricv4<TFixedImage, TMovingImage, TVirtualImage, TInternalComputationValueType, TMetricTraits>
  ::GetDerivative( DerivativeType & derivative ) const
  {
    MeasureType value;
    this->GetValueAndDerivative( value, derivative );
  }
  
  template<typename TFixedImage,typename TMovingImage,typename TVirtualImage, typename TInternalComputationValueType, typename TMetricTraits>
  void
  ImageToImageMetricv4<TFixedImage, TMovingImage, TVirtualImage, TInternalComputationValueType, TMetricTraits>
  ::GetValueAndDerivative( MeasureType & value, DerivativeType & derivative ) const
  {
    this->m_ComputeDerivative = true;
  
    this->m_DerivativeResult = &derivative;
    this->InitializeForIteration();
  
    // Do the threaded processing using the appropriate
    // GetValueAndDerivativeThreader. Results get written to
    // member vars.
    this->GetValueAndDerivativeExecute();
  
    value = this->m_Value;
  }
  
  template<typename TFixedImage,typename TMovingImage,typename TVirtualImage, typename TInternalComputationValueType, typename TMetricTraits>
  void
  ImageToImageMetricv4<TFixedImage, TMovingImage, TVirtualImage, TInternalComputationValueType, TMetricTraits>
  ::GetValueAndDerivativeExecute() const
  {
    if( this->m_UseFixedSampledPointSet ) // sparse sampling
      {
      SizeValueType numberOfPoints = this->GetNumberOfDomainPoints();
      if( numberOfPoints < 1 )
        {
        itkExceptionMacro("VirtualSampledPointSet must have 1 or more points.");
        }
      typename ImageToImageMetricv4GetValueAndDerivativeThreader< ThreadedIndexedContainerPartitioner, Self >::DomainType range;
      range[0] = 0;
      range[1] = numberOfPoints - 1;
      this->m_SparseGetValueAndDerivativeThreader->Execute( const_cast< Self* >(this), range );
      }
    else // dense sampling
      {
      this->m_DenseGetValueAndDerivativeThreader->Execute( const_cast< Self* >(this), this->GetVirtualRegion() );
      }
  }
  
  template<typename TFixedImage,typename TMovingImage,typename TVirtualImage, typename TInternalComputationValueType, typename TMetricTraits>
  void
  ImageToImageMetricv4<TFixedImage, TMovingImage, TVirtualImage, TInternalComputationValueType, TMetricTraits>
  ::InitializeForIteration() const
  {
    if( this->m_ComputeDerivative )
      {
      /* This size always comes from the active transform */
      const NumberOfParametersType globalDerivativeSize = this->GetNumberOfParameters();
      if( this->m_DerivativeResult->GetSize() != globalDerivativeSize )
        {
        this->m_DerivativeResult->SetSize( globalDerivativeSize );
        }
      /* Clear derivative final result. */
      this->m_DerivativeResult->Fill( NumericTraits< DerivativeValueType >::ZeroValue() );
      }
  }
  
  template<typename TFixedImage,typename TMovingImage,typename TVirtualImage, typename TInternalComputationValueType, typename TMetricTraits>
  bool
  ImageToImageMetricv4<TFixedImage, TMovingImage, TVirtualImage, TInternalComputationValueType, TMetricTraits>
  ::TransformAndEvaluateFixedPoint(
                           const VirtualPointType & virtualPoint,
                           FixedImagePointType & mappedFixedPoint,
                           FixedImagePixelType & mappedFixedPixelValue ) const
  {
    bool pointIsValid = true;
    mappedFixedPixelValue = NumericTraits<FixedImagePixelType>::ZeroValue();
  
    // map the point into fixed space
    this->LocalTransformPoint(virtualPoint,mappedFixedPoint);
  
    // check against the mask if one is assigned
    if ( this->m_FixedImageMask )
      {
      // Check if mapped point is within the support region of the fixed image
      // mask
      pointIsValid = this->m_FixedImageMask->IsInside( mappedFixedPoint );
      if( ! pointIsValid )
        {
        return pointIsValid;
        }
      }
  
    // Check if mapped point is inside image buffer
    pointIsValid = this->m_FixedInterpolator->IsInsideBuffer(mappedFixedPoint);
  
    // Evaluate
    if( pointIsValid )
      {
      mappedFixedPixelValue = this->m_FixedInterpolator->Evaluate(mappedFixedPoint);
      }
  
    return pointIsValid;
  }
  
  template<typename TFixedImage,typename TMovingImage,typename TVirtualImage, typename TInternalComputationValueType, typename TMetricTraits>
  bool
  ImageToImageMetricv4<TFixedImage, TMovingImage, TVirtualImage, TInternalComputationValueType, TMetricTraits>
  ::TransformAndEvaluateMovingPoint(
                           const VirtualPointType & virtualPoint,
                           MovingImagePointType & mappedMovingPoint,
                           MovingImagePixelType & mappedMovingPixelValue ) const
  {
    bool pointIsValid = true;
    mappedMovingPixelValue = NumericTraits<MovingImagePixelType>::ZeroValue();
  
    // map the point into moving space
  
    // Before transforming points, we should convert their types from the ImagePointType (aka Point<double, dim>)
    // to TransformPointType (aka Point<ScalarType, dim>).
    typename MovingTransformType::OutputPointType localVirtualPoint;
    typename MovingTransformType::OutputPointType localMappedMovingPoint;
  
    localVirtualPoint.CastFrom(virtualPoint);
    localMappedMovingPoint.CastFrom(mappedMovingPoint);
  
    localMappedMovingPoint = this->m_MovingTransform->TransformPoint( localVirtualPoint );
    mappedMovingPoint.CastFrom(localMappedMovingPoint);
  
    // check against the mask if one is assigned
    if ( this->m_MovingImageMask )
      {
      // Check if mapped point is within the support region of the fixed image
      // mask
      pointIsValid = this->m_MovingImageMask->IsInside( mappedMovingPoint );
      if( ! pointIsValid )
        {
        return pointIsValid;
        }
      }
  
    // Check if mapped point is inside image buffer
    pointIsValid = this->m_MovingInterpolator->IsInsideBuffer(mappedMovingPoint);
  
    // Evaluate
    if( pointIsValid )
      {
      mappedMovingPixelValue = this->m_MovingInterpolator->Evaluate( mappedMovingPoint );
      }
  
    return pointIsValid;
  }
  
  template<typename TFixedImage,typename TMovingImage,typename TVirtualImage, typename TInternalComputationValueType, typename TMetricTraits>
  void
  ImageToImageMetricv4<TFixedImage, TMovingImage, TVirtualImage, TInternalComputationValueType, TMetricTraits>
  ::ComputeFixedImageGradientAtPoint( const FixedImagePointType & mappedPoint, FixedImageGradientType & gradient ) const
  {
    if ( this->m_UseFixedImageGradientFilter )
      {
      if( ! this->GetGradientSourceIncludesFixed() )
        {
        itkExceptionMacro("Attempted to retrieve fixed image gradient from gradient image filter, "
                          "but GradientSource does not include 'fixed', and thus the gradient image has not been calculated.");
        }
      gradient = m_FixedImageGradientInterpolator->Evaluate( mappedPoint );
      }
    else
      {
      // if not using the gradient image
      gradient = this->m_FixedImageGradientCalculator->Evaluate( mappedPoint );
      }
  }
  
  template<typename TFixedImage,typename TMovingImage,typename TVirtualImage, typename TInternalComputationValueType, typename TMetricTraits>
  void
  ImageToImageMetricv4<TFixedImage, TMovingImage, TVirtualImage, TInternalComputationValueType, TMetricTraits>
  ::ComputeMovingImageGradientAtPoint( const MovingImagePointType & mappedPoint, MovingImageGradientType & gradient ) const
  {
    if ( this->m_UseMovingImageGradientFilter )
      {
      if( ! this->GetGradientSourceIncludesMoving() )
        {
        itkExceptionMacro("Attempted to retrieve moving image gradient from gradient image filter, "
                          "but GradientSource does not include 'moving', and thus the gradient image has not been calculated.");
        }
      gradient = m_MovingImageGradientInterpolator->Evaluate( mappedPoint );
      }
    else
      {
      // if not using the gradient image
      gradient = this->m_MovingImageGradientCalculator->Evaluate(mappedPoint);
      }
  }
  
  template<typename TFixedImage,typename TMovingImage,typename TVirtualImage, typename TInternalComputationValueType, typename TMetricTraits>
  void
  ImageToImageMetricv4<TFixedImage, TMovingImage, TVirtualImage, TInternalComputationValueType, TMetricTraits>
  ::ComputeFixedImageGradientFilterImage()
  {
    this->m_FixedImageGradientFilter->SetInput( this->m_FixedImage );
    this->m_FixedImageGradientFilter->Update();
    this->m_FixedImageGradientImage = this->m_FixedImageGradientFilter->GetOutput();
    this->m_FixedImageGradientInterpolator->SetInputImage( this->m_FixedImageGradientImage );
  }
  
  template<typename TFixedImage,typename TMovingImage,typename TVirtualImage, typename TInternalComputationValueType, typename TMetricTraits>
  void
  ImageToImageMetricv4<TFixedImage, TMovingImage, TVirtualImage, TInternalComputationValueType, TMetricTraits>
  ::ComputeMovingImageGradientFilterImage() const
  {
    this->m_MovingImageGradientFilter->SetInput( this->m_MovingImage );
    this->m_MovingImageGradientFilter->Update();
    this->m_MovingImageGradientImage = this->m_MovingImageGradientFilter->GetOutput();
    this->m_MovingImageGradientInterpolator->SetInputImage( this->m_MovingImageGradientImage );
  }
  
  template<typename TFixedImage,typename TMovingImage,typename TVirtualImage, typename TInternalComputationValueType, typename TMetricTraits>
  void
  ImageToImageMetricv4<TFixedImage, TMovingImage, TVirtualImage, TInternalComputationValueType, TMetricTraits>
  ::InitializeDefaultFixedImageGradientFilter()
  {
    const typename FixedImageType::SpacingType & spacing = this->m_FixedImage->GetSpacing();
    double maximumSpacing = 0.0;
    for ( ImageDimensionType i = 0; i < FixedImageDimension; i++ )
      {
      if ( spacing[i] > maximumSpacing )
        {
        maximumSpacing = spacing[i];
        }
      }
    this->m_DefaultFixedImageGradientFilter->SetSigma( maximumSpacing );
    this->m_DefaultFixedImageGradientFilter->SetNormalizeAcrossScale( true );
    this->m_DefaultFixedImageGradientFilter->SetNumberOfThreads( this->GetMaximumNumberOfThreads() );
    this->m_DefaultFixedImageGradientFilter->SetUseImageDirection( true );
  }
  
  template<typename TFixedImage,typename TMovingImage,typename TVirtualImage, typename TInternalComputationValueType, typename TMetricTraits>
  void
  ImageToImageMetricv4<TFixedImage, TMovingImage, TVirtualImage, TInternalComputationValueType, TMetricTraits>
  ::InitializeDefaultMovingImageGradientFilter()
  {
    const typename MovingImageType::SpacingType & spacing = this->m_MovingImage->GetSpacing();
    double maximumSpacing = 0.0;
    for ( ImageDimensionType i = 0; i < MovingImageDimension; i++ )
      {
      if ( spacing[i] > maximumSpacing )
        {
        maximumSpacing = spacing[i];
        }
      }
    this->m_DefaultMovingImageGradientFilter->SetSigma(maximumSpacing);
    this->m_DefaultMovingImageGradientFilter->SetNormalizeAcrossScale(true);
    this->m_DefaultMovingImageGradientFilter->SetNumberOfThreads(this->GetMaximumNumberOfThreads());
    this->m_DefaultMovingImageGradientFilter->SetUseImageDirection(true);
  }
  
  template<typename TFixedImage,typename TMovingImage,typename TVirtualImage, typename TInternalComputationValueType, typename TMetricTraits>
  void
  ImageToImageMetricv4<TFixedImage, TMovingImage, TVirtualImage, TInternalComputationValueType, TMetricTraits>
  ::SetMaximumNumberOfThreads( const ThreadIdType number )
  {
    if( number != this->m_SparseGetValueAndDerivativeThreader->GetMaximumNumberOfThreads() )
      {
      this->m_SparseGetValueAndDerivativeThreader->SetMaximumNumberOfThreads( number );
      this->Modified();
      }
    if( number != this->m_DenseGetValueAndDerivativeThreader->GetMaximumNumberOfThreads() )
      {
      this->m_DenseGetValueAndDerivativeThreader->SetMaximumNumberOfThreads( number );
      this->Modified();
      }
  }
  
  template<typename TFixedImage,typename TMovingImage,typename TVirtualImage, typename TInternalComputationValueType, typename TMetricTraits>
  ThreadIdType
  ImageToImageMetricv4<TFixedImage, TMovingImage, TVirtualImage, TInternalComputationValueType, TMetricTraits>
  ::GetMaximumNumberOfThreads() const
  {
    if( this->m_UseFixedSampledPointSet )
      {
      return this->m_SparseGetValueAndDerivativeThreader->GetMaximumNumberOfThreads();
      }
    return  this->m_DenseGetValueAndDerivativeThreader->GetMaximumNumberOfThreads();
  }
  
  template<typename TFixedImage,typename TMovingImage,typename TVirtualImage, typename TInternalComputationValueType, typename TMetricTraits>
  ThreadIdType
  ImageToImageMetricv4<TFixedImage, TMovingImage, TVirtualImage, TInternalComputationValueType, TMetricTraits>
  ::GetNumberOfThreadsUsed() const
  {
    if( this->m_UseFixedSampledPointSet )
      {
      return this->m_SparseGetValueAndDerivativeThreader->GetNumberOfThreadsUsed();
      }
    else
      {
      return this->m_DenseGetValueAndDerivativeThreader->GetNumberOfThreadsUsed();
      }
  }
  
  template<typename TFixedImage,typename TMovingImage,typename TVirtualImage, typename TInternalComputationValueType, typename TMetricTraits>
  void
  ImageToImageMetricv4<TFixedImage, TMovingImage, TVirtualImage, TInternalComputationValueType, TMetricTraits>
  ::MapFixedSampledPointSetToVirtual()
  {
    this->m_VirtualSampledPointSet = VirtualPointSetType::New();
    this->m_VirtualSampledPointSet->Initialize();
  
    typedef typename FixedSampledPointSetType::PointsContainer PointsContainer;
    typename PointsContainer::ConstPointer
      points = this->m_FixedSampledPointSet->GetPoints();
    if( points.IsNull() )
      {
      itkExceptionMacro("Fixed Sample point set is empty.");
      }
    typename PointsContainer::ConstIterator fixedIt = points->Begin();
  
    typename FixedTransformType::InverseTransformBasePointer
      inverseTransform = this->m_FixedTransform->GetInverseTransform();
    if( inverseTransform.IsNull() )
      {
      itkExceptionMacro("Unable to get inverse transform for mapping sampled "
                        " point set.");
      }
  
    this->m_NumberOfSkippedFixedSampledPoints = 0;
    SizeValueType virtualIndex = 0;
    while( fixedIt != points->End() )
      {
      typename FixedSampledPointSetType::PointType point = inverseTransform->TransformPoint( fixedIt.Value() );
      typename VirtualImageType::IndexType tempIndex;
      /* Verify that the point is valid. We may be working with a resized virtual domain,
       * and a fixed sampled point list that was created before the resizing. */
      if( this->TransformPhysicalPointToVirtualIndex( point, tempIndex ) )
        {
        this->m_VirtualSampledPointSet->SetPoint( virtualIndex, point );
        virtualIndex++;
        }
      else
        {
        this->m_NumberOfSkippedFixedSampledPoints++;
        }
      ++fixedIt;
      }
    if( this->m_VirtualSampledPointSet->GetNumberOfPoints() == 0 )
      {
      itkExceptionMacro("The virtual sampled point set has zero points because "
                        "no fixed sampled points were within the virtual "
                        "domain after mapping. There are no points to evaulate.");
      }
  }
  
  template<typename TFixedImage,typename TMovingImage,typename TVirtualImage, typename TInternalComputationValueType, typename TMetricTraits>
  SizeValueType
  ImageToImageMetricv4<TFixedImage, TMovingImage, TVirtualImage, TInternalComputationValueType, TMetricTraits>
  ::GetNumberOfDomainPoints() const
  {
    if( this->m_UseFixedSampledPointSet )
      {
      //The virtual sampled point set holds the actual points
      // over which we're evaluating over.
      return this->m_VirtualSampledPointSet->GetNumberOfPoints();
      }
    else
      {
      typename VirtualImageType::RegionType region = this->GetVirtualRegion();
      return region.GetNumberOfPixels();
      }
  }
  
  template<typename TFixedImage,typename TMovingImage,typename TVirtualImage, typename TInternalComputationValueType, typename TMetricTraits>
  void
  ImageToImageMetricv4<TFixedImage, TMovingImage, TVirtualImage, TInternalComputationValueType, TMetricTraits>
  ::PrintSelf(std::ostream& os, Indent indent) const
  {
    Superclass::PrintSelf(os, indent);
  
    os << indent << "ImageToImageMetricv4: " << std::endl
       << indent << "GetUseFixedImageGradientFilter: " << this->GetUseFixedImageGradientFilter() << std::endl
       << indent << "GetUseMovingImageGradientFilter: " << this->GetUseMovingImageGradientFilter() << std::endl
       << indent << "UseFloatingPointCorrection: " << this->GetUseFloatingPointCorrection() << std::endl
       << indent << "FloatingPointCorrectionResolution: " << this->GetFloatingPointCorrectionResolution() << std::endl;
  
    itkPrintSelfObjectMacro( FixedImage );
    itkPrintSelfObjectMacro( MovingImage );
    itkPrintSelfObjectMacro( FixedTransform );
    itkPrintSelfObjectMacro( MovingTransform );
    itkPrintSelfObjectMacro( FixedImageMask );
    itkPrintSelfObjectMacro( MovingImageMask );
  
  }
  
  }//namespace itk
  
  #endif
  