/*
* Copyright (c) ICG. All rights reserved.
*
* Institute for Computer Graphics and Vision
* Graz, University of Technology / Austria
*
*
* 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.
*
*
* Project     : projects
* Module      : MIPViewer - plugins
* Class       : $RCSfile: DownAndUpsampleUtilities.h,v $
* Language    : C++
* Description : 
*
* Author     : Martin Urschler
* EMail      : urschler@icg.tu-graz.ac.at
* Date       : $Date: 2007-07-04 14:10:29 $
* Version    : $Revision: 1.15 $
* Full Id    : $Id: DownAndUpsampleUtilities.h,v 1.15 2007-07-04 14:10:29 urschler Exp $
*
*/


#ifndef __DOWN_AND_UPSAMPLE_UTILITIES_H__
#define __DOWN_AND_UPSAMPLE_UTILITIES_H__


#include "itkImage.h"
#include "itkImageRegionIterator.h"
#include "itkImageLinearIteratorWithIndex.h"


namespace DownAndUpsampleUtilities
{

   enum _PyramidDownsampleMode {
      EVEN_TO_EVEN,
      EVEN_TO_ODD,
      ODD_TO_EVEN,
      ODD_TO_ODD
   };
  
   typedef enum _PyramidDownsampleMode PyramidDownsampleMode;
   
   typedef std::vector<PyramidDownsampleMode> PyramidDownsampleModeNTuple;
   
   typedef std::vector<PyramidDownsampleModeNTuple> 
      PyramidDownsampleModeSequence;
   
   

template< typename TInputImageType, typename TRealOutputImageType >
typename TRealOutputImageType::Pointer
lineBasedConvolutionWithShrinking(
   const typename TInputImageType::Pointer& inputImage,
   const unsigned int dim,
   const float *const operatorKernel,
   const unsigned int operatorSize,
   const unsigned int operatorRadius )
{  
   typename TInputImageType::SizeType originalSize = 
      inputImage->GetLargestPossibleRegion().GetSize();
      
//   typename TInputImageType::IndexType index;
//   for (unsigned int y=0; y<originalSize[1]; ++y)
//   {
//      for (unsigned int x=0; x<originalSize[0]; ++x)
//      {
//         index[0] = x;
//         index[1] = y;
//         std::cout << x << ";" << y << ":" << inputImage->GetPixel(index) << 
//            " ... ";
//      }
//      std::cout << std::endl;
//   }
      
   typename TRealOutputImageType::SizeType shrinkedSize = originalSize; 
   shrinkedSize[dim] = 1+(shrinkedSize[dim]-1)/2; 
   
   typename TRealOutputImageType::SpacingType shrinkedSpacing =
      inputImage->GetSpacing();
   shrinkedSpacing[dim] *= 2.0;
   
//   std::cout << "intermediate original size: " << originalSize << std::endl;
//   std::cout << "intermediate shrinked size: " << shrinkedSize << std::endl;
 
   // allocate the intermediate image
   typename TRealOutputImageType::RegionType shrinkedRegion =
      inputImage->GetLargestPossibleRegion();
   shrinkedRegion.SetSize( shrinkedSize );
   
   // create the downsampled image
   typename TRealOutputImageType::Pointer shrinkedImage = 
      TRealOutputImageType::New(); 
   shrinkedImage->SetRegions( shrinkedRegion );
   shrinkedImage->SetSpacing( shrinkedSpacing ); 
   shrinkedImage->SetOrigin( inputImage->GetOrigin() ); 
   shrinkedImage->Allocate();
   
   
   // create a scanline buffer for scanlines over current dimension
   
   const unsigned int sizeOfInputScanlineBuffer = 
      originalSize[dim] + 2 * operatorRadius;
   float *scanlineBuffer = new float[sizeOfInputScanlineBuffer];
   
   const unsigned int sizeOfResultBuffer = shrinkedSize[dim];
   float *resultBuffer = new float[sizeOfResultBuffer];

   itk::ImageLinearIteratorWithIndex<TInputImageType> inIter(
      inputImage, inputImage->GetLargestPossibleRegion() );
   inIter.SetDirection(dim);
   
   itk::ImageLinearIteratorWithIndex<TRealOutputImageType> outIter(
      shrinkedImage, shrinkedImage->GetLargestPossibleRegion() );
   outIter.SetDirection(dim);
   
   while (!inIter.IsAtEnd())
   {
      // 1. initialize the scanline buffer
      
      inIter.GoToBeginOfLine();
      
      // fill the left padding values
      unsigned int counter = 0;
      while (counter < operatorRadius)
      {
         scanlineBuffer[counter++] = inIter.Get();
      }
      
      // fetch the values from the iterator
      while (!inIter.IsAtEndOfLine())
      {
         scanlineBuffer[counter++] = inIter.Get();
         ++inIter;
      }
      
      // fill the right padding values
      const unsigned int lastIndex = 
         sizeOfInputScanlineBuffer-operatorRadius-1;
      while (counter < sizeOfInputScanlineBuffer)
      {
         scanlineBuffer[counter++] = scanlineBuffer[lastIndex];
      }
      
      // 2. initialize the result buffer
      
      std::fill( resultBuffer, &resultBuffer[sizeOfResultBuffer], 0.0f );
      
      // 3. now perform the convolution
      
      for (unsigned int outBufferPos=0; outBufferPos<sizeOfResultBuffer; 
           ++outBufferPos)
      {
         unsigned int inBufferStartPos = 2*outBufferPos;
//         std::cout << "inBufferPos start: " << inBufferStartPos << std::endl;
         
         for(unsigned int inBufferPos = 0; inBufferPos < operatorSize; 
             ++inBufferPos)
         {
            resultBuffer[outBufferPos] +=
               operatorKernel[inBufferPos%operatorSize] * 
               scanlineBuffer[inBufferStartPos+inBufferPos];
               
//            std::cout << "inBufferPos: " << inBufferStartPos+inBufferPos << 
//               " scanline: " << 
//               scanlineBuffer[inBufferStartPos+inBufferPos] << 
//               " outBufferPos " << outBufferPos 
//               << " intermed result: " << resultBuffer[outBufferPos] << 
//               std::endl; 
         }
         
//         std::cout << "result: " << resultBuffer[outBufferPos] << std::endl;
      }
      
      // 4. store the result in the result image
      
      outIter.GoToBeginOfLine();
      counter = 0;
      while (!outIter.IsAtEndOfLine())
      {
         outIter.Set( resultBuffer[counter++] );
         ++outIter;
      }
      
      inIter.NextLine();
      outIter.NextLine();
   }
   
   delete [] scanlineBuffer;
   delete [] resultBuffer;
   
//   for (unsigned int y=0; y<shrinkedSize[1]; ++y)
//   {
//      for (unsigned int x=0; x<shrinkedSize[0]; ++x)
//      {
//         index[0] = x;
//         index[1] = y;
//         std::cout << x << ";" << y << ":" << shrinkedImage->GetPixel(index) << 
//            " ... ";
//      }
//      std::cout << std::endl;
//   }
   
   return shrinkedImage;
}




// this method shrinks an image to half by first applying a smoothing with 
// a simple gauss approximation (binomial kernels) and sampling every second 
// pixel
// the odd 1D kernel is: 1/4 * (1 2 1)
// the even 1D kernel is: 1/8 * (1 3 3 1)
// for smoothing the borders are clamped

template< typename TImageType >
typename TImageType::Pointer shrinkImageToHalfByBinomialSmoothingImpl( 
   const typename TImageType::Pointer image,
   PyramidDownsampleModeNTuple& downsampleModeNTuple )
{
   typename TImageType::SizeType originalImageSize =
      image->GetLargestPossibleRegion().GetSize();

   std::cout << "originalImageSize: " << originalImageSize << std::endl;
   std::cout << "originalImageSpacing: " << image->GetSpacing() << std::endl;
   
   const unsigned int Dimensions = 
      itk::GetImageDimension<TImageType>::ImageDimension;
      
   // main loop
   
   typename TImageType::Pointer shrinkedImage = 0;
   
   for (unsigned int dim = 0; dim < Dimensions; ++dim)
   {
      // first we make a decision which of the four possible cases 
      // we have
      
      PyramidDownsampleMode downsampleMode;
      
      if (originalImageSize[dim] % 2 == 0)
      {
         if ( (originalImageSize[dim] / 2) % 2 == 0 )
            downsampleMode = EVEN_TO_EVEN;
         else
            downsampleMode = EVEN_TO_ODD;
      }
      else
      {
         if ( (( (originalImageSize[dim]-1) / 2) + 1) % 2 == 0 )
            downsampleMode = ODD_TO_EVEN;
         else
            downsampleMode = ODD_TO_ODD;
      }
      
      downsampleModeNTuple.push_back(downsampleMode);
      
      // now we perform the actual downsampling
      
      if (originalImageSize[dim] % 2 == 0) // the even case
      {
         const unsigned int operatorSize = 4;
         const float operatorKernel[operatorSize] =  
            { 0.125f, 0.375f, 0.375f, 0.125f };
         //const unsigned int operatorSize = 6;
         //const float operatorKernel[operatorSize] =  // kernel 1 5 10 10 5 1
         //   { 0.03125f, 0.15625f, 0.3125f, 0.3125f, 0.15625f, 0.03125f };
         
         const unsigned int operatorRadius = operatorSize/2 - 1;
         
         if (dim == 0)
         {
            shrinkedImage = lineBasedConvolutionWithShrinking<
               TImageType,TImageType>( 
                  image, dim, operatorKernel, operatorSize, 
                  operatorRadius );
         }
         else
         {
            shrinkedImage = lineBasedConvolutionWithShrinking<
               TImageType,TImageType>( 
                  shrinkedImage, dim, operatorKernel, operatorSize, 
                  operatorRadius );
         } 
      }
      else // the odd case
      {
         const unsigned int operatorSize = 3;
         const float operatorKernel[operatorSize] = 
            { 0.25f, 0.5f, 0.25f };
         //const unsigned int operatorSize = 5;
         //const float operatorKernel[operatorSize] =  // kernel 1 4 6 4 1
         //   { 0.0625f, 0.25f, 0.375f, 0.25f, 0.0625f };
         
         const unsigned int operatorRadius = operatorSize/2;
         
         if (dim == 0)
         {
            shrinkedImage = lineBasedConvolutionWithShrinking<
               TImageType,TImageType>( 
                  image, dim, operatorKernel, operatorSize, 
                  operatorRadius );
         }
         else
         {
            shrinkedImage = lineBasedConvolutionWithShrinking<
               TImageType,TImageType>( 
                  shrinkedImage, dim, operatorKernel, operatorSize, 
                  operatorRadius );
         } 
      }
   }

   std::cout << "shrinkedImageSize: " << 
      shrinkedImage->GetLargestPossibleRegion().GetSize() << std::endl;
   
   return shrinkedImage;
}

// FIXXME: the SameType considerations are still wrong!!!

template< typename TInputImageType, typename TOutputImageType >
typename TOutputImageType::Pointer shrinkImageToHalfByBinomialSmoothing( 
   const typename TInputImageType::Pointer image,
   PyramidDownsampleModeNTuple& downsampleModeNTuple )
{
   if (SameType<TInputImageType,TOutputImageType>::value == true)
   {
      return shrinkImageToHalfByBinomialSmoothingImpl<TOutputImageType>(
         image, downsampleModeNTuple);
   }
   else
   {
      typename TInputImageType::Pointer shrinkedImage = 
         shrinkImageToHalfByBinomialSmoothingImpl<TInputImageType>(
            image, downsampleModeNTuple);
   
      // create final result image
      ITK_IMAGE_ALLOCATE_MACRO_TN( TOutputImageType, result,
                                   TInputImageType, shrinkedImage );
   
      ITK_IMAGE_COPY_MACRO_TN( TInputImageType, shrinkedImage, 
                               TOutputImageType, result );

      return result;
   }
}

//template< typename TImageType >
//typename TImageType::Pointer shrinkImageToHalfByBinomialSmoothing( 
//   const typename TImageType::Pointer image,
//   PyramidDownsampleModeNTuple& downsampleModeNTuple )
//{
//   return shrinkImageToHalfByBinomialSmoothingImpl<TImageType>(
//      image, downsampleModeNTuple);
//}


template< typename TInputImageType, typename TRealOutputImageType >
typename TRealOutputImageType::Pointer
lineBasedInterpolationWithDoublingOddCase(
   const typename TInputImageType::Pointer& inputImage,
   const typename TRealOutputImageType::SizeType& doubledSize,
   const unsigned int dim )
{  
   typename TInputImageType::SizeType originalSize = 
      inputImage->GetLargestPossibleRegion().GetSize();
      
//   typename TInputImageType::IndexType index;
//   for (unsigned int y=0; y<originalSize[1]; ++y)
//   {
//      for (unsigned int x=0; x<originalSize[0]; ++x)
//      {
//         index[0] = x;
//         index[1] = y;
//         std::cout << x << ";" << y << ":" << inputImage->GetPixel(index) << 
//            " ... ";
//      }
//      std::cout << std::endl;
//   }
   
   typename TRealOutputImageType::SpacingType doubledSpacing =
      inputImage->GetSpacing();
   doubledSpacing[dim] *= 0.5;
   
   //std::cout << "intermediate original size: " << originalSize << std::endl;
   //std::cout << "intermediate doubled size: " << doubledSize << std::endl;
 
   // allocate the intermediate image
   
   typename TRealOutputImageType::RegionType doubledRegion =
      inputImage->GetLargestPossibleRegion();
   doubledRegion.SetSize( doubledSize );
   
   // create the upsampled image
   typename TRealOutputImageType::Pointer doubledImage = 
      TRealOutputImageType::New(); 
   doubledImage->SetRegions( doubledRegion );
   doubledImage->SetSpacing( doubledSpacing ); 
   doubledImage->SetOrigin( inputImage->GetOrigin() ); 
   doubledImage->Allocate();
   
   
   // create the result buffer that is used for the calculations
   
   const unsigned int sizeOfResultBuffer = doubledSize[dim];
   float *resultBuffer = new float[sizeOfResultBuffer];

   itk::ImageLinearIteratorWithIndex<TInputImageType> inIter(
      inputImage, inputImage->GetLargestPossibleRegion() );
   inIter.SetDirection(dim);
   
   itk::ImageLinearIteratorWithIndex<TRealOutputImageType> outIter(
      doubledImage, doubledImage->GetLargestPossibleRegion() );
   outIter.SetDirection(dim);
   
   while (!inIter.IsAtEnd())
   {
      inIter.GoToBeginOfLine();
      
      // 1. fill every second entry in result buffer with original values  
      
      unsigned int counter = 0;
      while (!inIter.IsAtEndOfLine())
      {
         resultBuffer[2*counter] = inIter.Get();
          
         ++counter;
         ++inIter;
      }
      
      // 2. fill in the rest of the entries in resultBuffer by linear 
      // interpolation
      
      for(counter = 1; counter<(sizeOfResultBuffer-1); counter=counter+2)
      {
         resultBuffer[counter] = 
            0.5f * (resultBuffer[counter-1]+resultBuffer[counter+1]);
      }
      
      // 3. store the result in the result image
      
      outIter.GoToBeginOfLine();
      counter = 0;
      while (!outIter.IsAtEndOfLine())
      {
         outIter.Set( resultBuffer[counter++] );
         ++outIter;
      }
      
      inIter.NextLine();
      outIter.NextLine();
   }
   
   delete [] resultBuffer;
   
//   for (unsigned int y=0; y<shrinkedSize[1]; ++y)
//   {
//      for (unsigned int x=0; x<shrinkedSize[0]; ++x)
//      {
//         index[0] = x;
//         index[1] = y;
//         std::cout << x << ";" << y << ":" << shrinkedImage->GetPixel(index) << 
//            " ... ";
//      }
//      std::cout << std::endl;
//   }
   
   return doubledImage;
}


template< typename TInputImageType, typename TRealOutputImageType >
typename TRealOutputImageType::Pointer
lineBasedInterpolationWithDoublingEvenCase(
   const typename TInputImageType::Pointer& inputImage,
   const unsigned int dim )
{  
   typename TInputImageType::SizeType originalSize = 
      inputImage->GetLargestPossibleRegion().GetSize();
      
//   typename TInputImageType::IndexType index;
//   for (unsigned int y=0; y<originalSize[1]; ++y)
//   {
//      for (unsigned int x=0; x<originalSize[0]; ++x)
//      {
//         index[0] = x;
//         index[1] = y;
//         std::cout << x << ";" << y << ":" << inputImage->GetPixel(index) << 
//            " ... ";
//      }
//      std::cout << std::endl;
//   }
      
      
   typename TRealOutputImageType::SizeType doubledSize = originalSize; 
   doubledSize[dim] = 2*originalSize[dim];
   
   typename TRealOutputImageType::SpacingType doubledSpacing =
      inputImage->GetSpacing();
   doubledSpacing[dim] *= 0.5;
   
   //std::cout << "intermediate original size: " << originalSize << std::endl;
   //std::cout << "intermediate doubled size: " << doubledSize << std::endl;
 
   // allocate the intermediate image
   
   typename TRealOutputImageType::RegionType doubledRegion =
      inputImage->GetLargestPossibleRegion();
   doubledRegion.SetSize( doubledSize );
   
   // create the upsampled image
   typename TRealOutputImageType::Pointer doubledImage = 
      TRealOutputImageType::New(); 
   doubledImage->SetRegions( doubledRegion );
   doubledImage->SetSpacing( doubledSpacing ); 
   doubledImage->SetOrigin( inputImage->GetOrigin() ); 
   doubledImage->Allocate();
   
   
   // create the result buffer that is used for the calculations
   
   const unsigned int sizeOfResultBuffer = doubledSize[dim];
   float *resultBuffer = new float[sizeOfResultBuffer];

   itk::ImageLinearIteratorWithIndex<TInputImageType> inIter(
      inputImage, inputImage->GetLargestPossibleRegion() );
   inIter.SetDirection(dim);
   
   itk::ImageLinearIteratorWithIndex<TRealOutputImageType> outIter(
      doubledImage, doubledImage->GetLargestPossibleRegion() );
   outIter.SetDirection(dim);
   
   while (!inIter.IsAtEnd())
   {
      inIter.GoToBeginOfLine();
      
      // 1. calculate the first value of the result buffer
      unsigned int counter = 0;
      float previousInputValue = inIter.Get(); 
      resultBuffer[counter++] = previousInputValue; 
      ++inIter;
      
      float currentInputValue = 0.0f;
      
      // 2. calculate all interior values
      for (counter=1;counter<(sizeOfResultBuffer-1); ++counter)
      {
         currentInputValue = inIter.Get();
         
         resultBuffer[counter] = (counter%2 == 1) ? 
            (0.75f*previousInputValue + 0.25f*currentInputValue) : 
            (0.25f*previousInputValue + 0.75f*currentInputValue);
         
         if (counter%2 == 0)
         {
            previousInputValue = currentInputValue;
            ++inIter;
         }
      }
      
      // 3. calculate the last entry
      resultBuffer[counter] = currentInputValue;
      
      // 4. store the result in the result image
      
      outIter.GoToBeginOfLine();
      counter = 0;
      while (!outIter.IsAtEndOfLine())
      {
         outIter.Set( resultBuffer[counter++] );
         ++outIter;
      }
      
      inIter.NextLine();
      outIter.NextLine();
   }
   
   delete [] resultBuffer;
   
//   for (unsigned int y=0; y<shrinkedSize[1]; ++y)
//   {
//      for (unsigned int x=0; x<shrinkedSize[0]; ++x)
//      {
//         index[0] = x;
//         index[1] = y;
//         std::cout << x << ";" << y << ":" << shrinkedImage->GetPixel(index) << 
//            " ... ";
//      }
//      std::cout << std::endl;
//   }
   
   return doubledImage;
}






template< typename TImageType >
typename TImageType::Pointer doubleImageByLinearInterpolationImpl( 
   const typename TImageType::Pointer image,
   const PyramidDownsampleModeNTuple& downsampleModeNTuple )
{
   typename TImageType::SizeType originalImageSize =
      image->GetLargestPossibleRegion().GetSize();

   std::cout << "originalImageSize: " << originalImageSize << std::endl;
   std::cout << "originalImageSpacing: " << image->GetSpacing() << std::endl;
   
   const unsigned int dimensions = 
      itk::GetImageDimension<TImageType>::ImageDimension;
      
   // main loop
   
   typename TImageType::Pointer doubledImage = 0;
   
   for (unsigned int dim = 0; dim < dimensions; ++dim)
   {
      // sanity checks
      if (originalImageSize[dim] % 2 == 0 && // the even case
          !(downsampleModeNTuple[dim] == EVEN_TO_EVEN || 
            downsampleModeNTuple[dim] == ODD_TO_EVEN) )
      {
         std::cout << "doubleImageByLinearInterpolation: inconsistency even!" 
            << std::endl;
      }
      if (originalImageSize[dim] % 2 == 1 && // the odd case
          !(downsampleModeNTuple[dim] == EVEN_TO_ODD || 
            downsampleModeNTuple[dim] == ODD_TO_ODD) )
      {
         std::cout << "doubleImageByLinearInterpolation: inconsistency odd!" 
            << std::endl;
      }
      
      typename TImageType::Pointer workImage = (dim==0) ? image : doubledImage; 
      
      switch (downsampleModeNTuple[dim])
      {
         case EVEN_TO_EVEN:
         {
            doubledImage = lineBasedInterpolationWithDoublingEvenCase<
               TImageType, TImageType>( workImage, dim );
            break;
         }
         case EVEN_TO_ODD:
         {
            doubledImage = lineBasedInterpolationWithDoublingEvenCase<
               TImageType, TImageType>( workImage, dim );
            break;
         }
         case ODD_TO_EVEN:
         {
            typename TImageType::SizeType doubledSize = 
               workImage->GetLargestPossibleRegion().GetSize(); 
            doubledSize[dim] = 2*originalImageSize[dim]-1;
   
            doubledImage = lineBasedInterpolationWithDoublingOddCase<
               TImageType, TImageType>( workImage, doubledSize, dim );
            break;
         }
         case ODD_TO_ODD:
         {
            typename TImageType::SizeType doubledSize = 
               workImage->GetLargestPossibleRegion().GetSize();
            doubledSize[dim] = 1+2*(originalImageSize[dim]-1);
            
            doubledImage = lineBasedInterpolationWithDoublingOddCase<
               TImageType, TImageType>( workImage, doubledSize, dim );
            break;
         }
         default:
         {
            std::cerr << "Unknown downsample mode: " << 
               downsampleModeNTuple[dim] << std::endl;
            break;
         }
      }
   }
   
   std::cout << "doubledImageSize: " << 
      doubledImage->GetLargestPossibleRegion().GetSize() << std::endl;
   
   return doubledImage;
}
 
template< typename TInputImageType, typename TOutputImageType >
typename TOutputImageType::Pointer doubleImageByLinearInterpolation( 
   const typename TInputImageType::Pointer image,
   const PyramidDownsampleModeNTuple& downsampleModeNTuple )
{
   if (SameType<TInputImageType,TOutputImageType>::value == true)
   {
      return doubleImageByLinearInterpolationImpl<TOutputImageType>(image,
         downsampleModeNTuple);
   }
   else
   {
      typename TInputImageType::Pointer doubledImage = 
         doubleImageByLinearInterpolationImpl<TInputImageType>(image,
            downsampleModeNTuple);
   
      // create final result image
      ITK_IMAGE_ALLOCATE_MACRO_TN( TOutputImageType, result,
                                   TInputImageType, doubledImage );
   
      // copy intermediate result to final result
      ITK_IMAGE_COPY_MACRO_TN( TInputImageType, doubledImage, 
                               TOutputImageType, result );

      return result;
   }
}
 
 
 
 
 
 
 
template< typename TVectorFieldType >
typename TVectorFieldType::Pointer
lineBasedVectorFieldInterpolationWithDoublingOddCase(
   const typename TVectorFieldType::Pointer& inputImage,
   const typename TVectorFieldType::SizeType& doubledSize,
   const unsigned int dim )
{  
   typename TVectorFieldType::SizeType originalSize = 
      inputImage->GetLargestPossibleRegion().GetSize();

   typename TVectorFieldType::SpacingType doubledSpacing =
      inputImage->GetSpacing();
   doubledSpacing[dim] *= 0.5;
   
   //std::cout << "intermediate original size: " << originalSize << std::endl;
   //std::cout << "intermediate doubled size: " << doubledSize << std::endl;
 
   // allocate the intermediate image
   
   typename TVectorFieldType::RegionType doubledRegion =
      inputImage->GetLargestPossibleRegion();
   doubledRegion.SetSize( doubledSize );
   
   // create the upsampled image
   typename TVectorFieldType::Pointer doubledImage = 
      TVectorFieldType::New(); 
   doubledImage->SetRegions( doubledRegion );
   doubledImage->SetSpacing( doubledSpacing ); 
   doubledImage->SetOrigin( inputImage->GetOrigin() ); 
   doubledImage->Allocate();
   
   
   // create the result buffer that is used for the calculations
   
   const unsigned int sizeOfResultBuffer = doubledSize[dim];
   typename TVectorFieldType::PixelType *resultBuffer = 
      new typename TVectorFieldType::PixelType[sizeOfResultBuffer];

   itk::ImageLinearIteratorWithIndex<TVectorFieldType> inIter(
      inputImage, inputImage->GetLargestPossibleRegion() );
   inIter.SetDirection(dim);
   
   itk::ImageLinearIteratorWithIndex<TVectorFieldType> outIter(
      doubledImage, doubledImage->GetLargestPossibleRegion() );
   outIter.SetDirection(dim);
   
   while (!inIter.IsAtEnd())
   {
      inIter.GoToBeginOfLine();
      
      // 1. fill every second entry in result buffer with original values  
      
      unsigned int counter = 0;
      while (!inIter.IsAtEndOfLine())
      {
         resultBuffer[2*counter] = inIter.Get();
          
         ++counter;
         ++inIter;
      }
      
      // 2. fill in the rest of the entries in resultBuffer by linear 
      // interpolation
      
      for(counter = 1; counter<(sizeOfResultBuffer-1); counter=counter+2)
      {
         resultBuffer[counter] = 
            0.5f * (resultBuffer[counter-1]+resultBuffer[counter+1]);
      }
      
      // 3. store the result in the result image
      
      outIter.GoToBeginOfLine();
      counter = 0;
      while (!outIter.IsAtEndOfLine())
      {
         outIter.Set( resultBuffer[counter++] );
         ++outIter;
      }
      
      inIter.NextLine();
      outIter.NextLine();
   }
   
   delete [] resultBuffer;
   
   return doubledImage;
}


template< typename TVectorFieldType >
typename TVectorFieldType::Pointer
lineBasedVectorFieldInterpolationWithDoublingEvenCase(
   const typename TVectorFieldType::Pointer& inputImage,
   const unsigned int dim )
{  
   typename TVectorFieldType::SizeType originalSize = 
      inputImage->GetLargestPossibleRegion().GetSize();
      
      
   typename TVectorFieldType::SizeType doubledSize = originalSize; 
   doubledSize[dim] = 2*originalSize[dim]; 
   
   typename TVectorFieldType::SpacingType doubledSpacing =
      inputImage->GetSpacing();
   doubledSpacing[dim] *= 0.5;
   
   //std::cout << "intermediate original size: " << originalSize << std::endl;
   //std::cout << "intermediate doubled size: " << doubledSize << std::endl;
 
   // allocate the intermediate image
   
   typename TVectorFieldType::RegionType doubledRegion =
      inputImage->GetLargestPossibleRegion();
   doubledRegion.SetSize( doubledSize );
   
   // create the upsampled image
   typename TVectorFieldType::Pointer doubledImage = 
      TVectorFieldType::New(); 
   doubledImage->SetRegions( doubledRegion );
   doubledImage->SetSpacing( doubledSpacing ); 
   doubledImage->SetOrigin( inputImage->GetOrigin() ); 
   doubledImage->Allocate();
   
   
   // create the result buffer that is used for the calculations
   
   const unsigned int sizeOfResultBuffer = doubledSize[dim];
   typename TVectorFieldType::PixelType *resultBuffer = 
      new typename TVectorFieldType::PixelType[sizeOfResultBuffer];

   itk::ImageLinearIteratorWithIndex<TVectorFieldType> inIter(
      inputImage, inputImage->GetLargestPossibleRegion() );
   inIter.SetDirection(dim);
   
   itk::ImageLinearIteratorWithIndex<TVectorFieldType> outIter(
      doubledImage, doubledImage->GetLargestPossibleRegion() );
   outIter.SetDirection(dim);
   
   while (!inIter.IsAtEnd())
   {
      inIter.GoToBeginOfLine();
      
      // 1. calculate the first value of the result buffer
      unsigned int counter = 0;
      typename TVectorFieldType::PixelType previousInputValue = inIter.Get(); 
      resultBuffer[counter++] = previousInputValue; 
      ++inIter;
      
      typename TVectorFieldType::PixelType currentInputValue;
      
      // 2. calculate all interior values
      for (counter=1;counter<(sizeOfResultBuffer-1); ++counter)
      {
         currentInputValue = inIter.Get();
         
         resultBuffer[counter] = (counter%2 == 1) ? 
            (0.75f*previousInputValue + 0.25f*currentInputValue) : 
            (0.25f*previousInputValue + 0.75f*currentInputValue); 
         
         if (counter%2 == 0)
         {
            previousInputValue = currentInputValue;
            ++inIter;
         }
      }
      
      // 3. calculate the last entry
      resultBuffer[counter] = currentInputValue;
      
      // 4. store the result in the result image
      
      outIter.GoToBeginOfLine();
      counter = 0;
      while (!outIter.IsAtEndOfLine())
      {
         outIter.Set( resultBuffer[counter++] );
         ++outIter;
      }
      
      inIter.NextLine();
      outIter.NextLine();
   }
   
   delete [] resultBuffer;
   
   return doubledImage;
}
 
template< typename TImageType >
typename TImageType::Pointer doubleVectorFieldByLinearInterpolation( 
   const typename TImageType::Pointer image,
   const PyramidDownsampleModeNTuple& downsampleModeNTuple )
{
   typename TImageType::SizeType originalImageSize =
      image->GetLargestPossibleRegion().GetSize();
   
   const unsigned int dimensions = 
      itk::GetImageDimension<TImageType>::ImageDimension;
      
   // main loop
   
   typename TImageType::Pointer doubledImage = 0;
   
   for (unsigned int dim = 0; dim < dimensions; ++dim)
   {
      // sanity checks
      if (originalImageSize[dim] % 2 == 0 && // the even case
          !(downsampleModeNTuple[dim] == EVEN_TO_EVEN || 
            downsampleModeNTuple[dim] == ODD_TO_EVEN) )
      {
         std::cout << "doubleDefFieldByLinearInterpolation: inconsistency even!"
            << std::endl;
      }
      if (originalImageSize[dim] % 2 == 1 && // the odd case
          !(downsampleModeNTuple[dim] == EVEN_TO_ODD || 
            downsampleModeNTuple[dim] == ODD_TO_ODD) )
      {
         std::cout << "doubleDefFieldByLinearInterpolation: inconsistency odd!" 
            << std::endl;
      }
      
      typename TImageType::Pointer workImage = (dim==0) ? image : doubledImage; 
      
      switch (downsampleModeNTuple[dim])
      {
         case EVEN_TO_EVEN:
         {
            doubledImage = lineBasedVectorFieldInterpolationWithDoublingEvenCase
               <TImageType>( workImage, dim );
            break;
         }
         case EVEN_TO_ODD:
         {
            doubledImage = lineBasedVectorFieldInterpolationWithDoublingEvenCase
               <TImageType>( workImage, dim );
            break;
         }
         case ODD_TO_EVEN:
         {
            typename TImageType::SizeType doubledSize = 
               workImage->GetLargestPossibleRegion().GetSize();
            doubledSize[dim] = 2*originalImageSize[dim]-1;
            
            doubledImage = lineBasedVectorFieldInterpolationWithDoublingOddCase
               <TImageType>( workImage, doubledSize, dim );
            break;
         }
         case ODD_TO_ODD:
         {
            typename TImageType::SizeType doubledSize = 
               workImage->GetLargestPossibleRegion().GetSize();
            doubledSize[dim] = 1+2*(originalImageSize[dim]-1);
            
            doubledImage = lineBasedVectorFieldInterpolationWithDoublingOddCase
               <TImageType>( workImage, doubledSize, dim );
            break;
         }
         default:
         {
            std::cerr << "Unknown downsample mode: " << 
               downsampleModeNTuple[dim] << std::endl;
            break;
         }
      }
   }
      
   return doubledImage;
}  
   

} // end of namespace

#endif // __DOWN_AND_UPSAMPLE_UTILITIES_H__