itkBSplineResampleImageFilterBase.txx

Go to the documentation of this file.

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

  Program:   Insight Segmentation & Registration Toolkit
  Module:    $RCSfile: itkBSplineResampleImageFilterBase.txx,v $
  Language:  C++
  Date:      $Date: 2009-05-11 21:42:52 $
  Version:   $Revision: 1.12 $

  Copyright (c) Insight Software Consortium. All rights reserved.
  See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.

  Portions of this code are covered under the VTK copyright.
  See VTKCopyright.txt or http://www.kitware.com/VTKCopyright.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 __itkBSplineResampleImageFilterBase_txx
#define __itkBSplineResampleImageFilterBase_txx

#include "itkBSplineResampleImageFilterBase.h"
#include "itkProgressReporter.h"

namespace itk
{

template <class TInputImage, class TOutputImage>
BSplineResampleImageFilterBase<TInputImage, TOutputImage>
::BSplineResampleImageFilterBase()
{
  m_SplineOrder = -1;
  int SplineOrder = 0;
  // Because of inheritance the user must explicitly set this for m_SplineOrder != 0.
  this->SetSplineOrder(SplineOrder);
}

template <class TInputImage, class TOutputImage>
void
BSplineResampleImageFilterBase<TInputImage, TOutputImage>
::PrintSelf(
  std::ostream& os, 
  Indent indent) const
{
  Superclass::PrintSelf( os, indent );
  os << indent << "Spline Order: " << m_SplineOrder << std::endl;
}

template <class TInputImage, class TOutputImage>
void BSplineResampleImageFilterBase<TInputImage, TOutputImage>
::InitializePyramidSplineFilter(int SplineOrder)
{
  switch (SplineOrder) 
    {
  
    case 0 :
      m_GSize = 1; 
      m_HSize = 1;
      break;
      
    case 1 :
      m_GSize = 9;
      m_HSize = 2;
      m_G.resize(m_GSize);
      m_H.resize(m_HSize);
      m_G[0]  =  0.707107; 
      m_G[1]  =  0.292893; 
      m_G[2]  = -0.12132; 
      m_G[3]  = -0.0502525;
      m_G[4]  =  0.0208153; 
      m_G[5]  =  0.00862197; 
      m_G[6]  = -0.00357134;
      m_G[7]  = -0.0014793; 
      m_G[8]  =  0.000612745;
      m_H[0]  = 1.; 
      m_H[1]  = 0.5;
      break;
    case 2 :
      m_GSize = 16;
      m_HSize = 10;
      m_G.resize(m_GSize);
      m_H.resize(m_HSize);
      m_G[0]  =  0.617317; 
      m_G[1]  =  0.310754; 
      m_G[2]  = -0.0949641; 
      m_G[3]  = -0.0858654;
      m_G[4]  =  0.0529153; 
      m_G[5]  =  0.0362437; 
      m_G[6]  = -0.0240408;
      m_G[7]  = -0.0160987; 
      m_G[8]  =  0.0107498; 
      m_G[9]  =  0.00718418;
      m_G[10] = -0.00480004; 
      m_G[11] = -0.00320734; 
      m_G[12] =  0.00214306;
      m_G[13] =  0.00143195; 
      m_G[14] = -0.0009568; 
      m_G[15] = -0.000639312;
      m_H[0]  =  1.; 
      m_H[1]  =  0.585786; 
      m_H[2]  =  0; 
      m_H[3]  = -0.100505; 
      m_H[4]  =  0;
      m_H[5]  =  0.0172439; 
      m_H[6]  =  0; 
      m_H[7]  = -0.00295859; 
      m_H[8]  =  0;
      m_H[9]  =  0.000507614;
      break;
    case 3 :
      m_GSize = 20;
      m_HSize = 12;
      m_G.resize(m_GSize);
      m_H.resize(m_HSize);
      m_G[0]  =  0.596797; 
      m_G[1]  =  0.313287; 
      m_G[2]  = -0.0827691; 
      m_G[3]  = -0.0921993;
      m_G[4]  =  0.0540288; 
      m_G[5]  =  0.0436996; 
      m_G[6]  = -0.0302508;
      m_G[7]  = -0.0225552; 
      m_G[8]  =  0.0162251; 
      m_G[9]  =  0.0118738;
      m_G[10] = -0.00861788; 
      m_G[11] = -0.00627964; 
      m_G[12] =  0.00456713;
      m_G[13] =  0.00332464; 
      m_G[14] = -0.00241916; 
      m_G[15] = -0.00176059;
      m_G[16] =  0.00128128; 
      m_G[17] =  0.000932349; 
      m_G[18] = -0.000678643;
      m_G[19] = -0.000493682;
      m_H[0]  =  1.; 
      m_H[1]  =  0.600481; 
      m_H[2]  =  0; 
      m_H[3]  = -0.127405; 
      m_H[4]  =  0;
      m_H[5]  =  0.034138; 
      m_H[6]  =  0; 
      m_H[7]  = -0.00914725; 
      m_H[8]  =  0;
      m_H[9]  =  0.002451; 
      m_H[10] =  0; 
      m_H[11] = -0.000656743;
      break;
    default :
      // Throw an exception
      ExceptionObject err(__FILE__, __LINE__);
      err.SetLocation( ITK_LOCATION );
      err.SetDescription( "SplineOrder for l2 pyramid filter must be between 1 and 3. Requested spline order has not been implemented." );
      throw err;
      break;
    }
}


template < class TInputImage, class TOutputImage>
void BSplineResampleImageFilterBase<TInputImage, TOutputImage>
::SetSplineOrder(int splineOrder)
{
  if (splineOrder == m_SplineOrder)
    {
    return;
    }
  m_SplineOrder = splineOrder;

  this->InitializePyramidSplineFilter(m_SplineOrder);
  this->Modified();
  
}

template <class TInputImage, class TOutputImage>
void BSplineResampleImageFilterBase<TInputImage, TOutputImage>
::Reduce1DImage( const std::vector<double> & in,   OutputImageIterator & out, 
                 unsigned int inTraverseSize, ProgressReporter &progress )
{
  
  int i1, i2;

  unsigned int outK, inK;
  unsigned int outTraverseSize = inTraverseSize/2;
  inTraverseSize = outTraverseSize*2;  // ensures that an even number is used.
  unsigned int inModK; // number for modulus math of in
  inModK = inTraverseSize - 1;


  double outVal;


  //TODO:  m_GSize < 2 has not been tested.
  if (m_GSize < 2)
    {   
    for (outK = 0; outK < outTraverseSize; outK++)
      {
      inK = 2 * outK;
      i2 = inK + 1;
      if (i2 > (int) inModK )
        {
        //Original was 
        //i2=inModK-i2;
        //I don't think this is correct since this would be negative
        i2 = inModK - (i2 % inModK);  // Should I use this always instead of the if statement?
        }
      out.Set( static_cast<OutputImagePixelType> ( ( in[inK] + in[i2] ) / 2.0 ) );
      ++out;
      progress.CompletedPixel();
      }
    }

  else
    {
    for (outK = 0; outK < outTraverseSize; outK++) 
      {
      inK = 2L * outK;

      outVal = in[inK] * m_G[0];
      
      for (int i = 1; i < m_GSize; i++) 
        {
        // Calculate indicies for left and right of symmetrical filter.
        i1 = inK - i;
        i2 = inK + i;
        // reflect at boundaries if necessary
        if (i1 < 0) 
          {
          i1 = (-i1) % inModK;
          // Removed because i1 can never be greater than inModK, right?
          //if (i1 > inModK) 
          //i1=inModK-i1;  //TODO: I don't think this is correct.
          }
        if (i2 > (int) inModK) 
          {
          i2 = i2 % inModK;
          // Removed because i1 can never be greater than inModK, right?
          //if (i2 > inModK) 
          //i2=inModK-i2;  //TODO: I don't think this is correct.
          }
        outVal = outVal + m_G[i]*(in[i1] + in[i2]);
        }
      out.Set( static_cast<OutputImagePixelType> (outVal) );
      ++out;
      progress.CompletedPixel();
      }
    }
}

template <class TInputImage, class TOutputImage>
void BSplineResampleImageFilterBase<TInputImage, TOutputImage>
::Expand1DImage( const std::vector<double> & in, OutputImageIterator & out, 
                 unsigned int inTraverseSize, ProgressReporter &progress )
{
  int i1, i2;

  int outK;
  unsigned int inK;
  unsigned int outTraverseSize = inTraverseSize * 2;
  //inTraverseSize = outTraverseSize/2;  // ensures that an even number is used.
  int inModK; // number for modulus math of in
  inModK = inTraverseSize - 1;


  double outVal;


  //TODO:  m_GSize < 2 has not been tested.
  if (m_HSize < 2)
    {   
    for (inK = 0; inK < inTraverseSize; inK++)
      {
      out.Set( static_cast<OutputImagePixelType> (in[inK]) );
      ++out;
      out.Set( static_cast<OutputImagePixelType> (in[inK]) );
      ++out;
      }
    progress.CompletedPixel();
    }

  else
    {
    for (outK = 0; outK < (int) outTraverseSize; outK++) 
      {
      outVal = 0.0;
      for (int k = (outK % 2); k < (int)  m_HSize; k += 2)
        {
        i1 = ( outK - k ) / 2;
        if ( i1 < 0 )
          {
          i1 = (-i1) % inModK;
          // The following could never happen therefore removed.
          //if (i1 > inModK)
          //i1 - inModK - i1;
          }
        outVal = outVal + m_H[k] * in[i1];
        }
      for( int k = 2 - ( outK % 2); k < (int) m_HSize; k += 2 )
        {
        i2 = (outK + k)/ 2;
        if ( i2 > inModK )
          {
          i2 = i2 % inModK;
          i2 = inModK - i2;
          // The following could never happen therefore removed.
          //if (i2 > inModK)
          //i2 - inModK - i2;
          }
        outVal += m_H[k] * in[i2];
        }
    
      out.Set( static_cast<OutputImagePixelType> (outVal) );
      ++out;
      progress.CompletedPixel();
      }
    }

}


template < class TInputImage, class TOutputImage>
void BSplineResampleImageFilterBase<TInputImage, TOutputImage>
::ReduceNDImage(OutputImageIterator &outItr)
{
  // Set up variables for waking the image regions.
  RegionType validRegion;
  SizeType startSize;
  SizeType currentSize;
  // Does not support streaming
  typename Superclass::InputImagePointer  inputPtr = const_cast< TInputImage * > ( this->GetInput() );
  startSize = inputPtr->GetBufferedRegion().GetSize();
  
  // Initilize scratchImage space and allocate memory
  InitializeScratch(startSize);
  typename TOutputImage::Pointer scratchImage;
  scratchImage =  TOutputImage::New();
  scratchImage->CopyInformation( inputPtr );
  RegionType scratchRegion;
  scratchRegion = inputPtr->GetBufferedRegion();
  currentSize = startSize;
  // scratchImage only needs the 1/2 the space of the original
  // image for the first dimension.
  // TODO: Is dividing by 2 correct or do I need something more complicated for handling odd dimensioned
  //       images?  i.e. is the rounding handled correctly?
  currentSize[0] = currentSize[0]/2;  
  scratchRegion.SetSize( currentSize );
  scratchImage->SetRegions( scratchRegion );
  scratchImage->Allocate();  
  
 
  currentSize = startSize;
  validRegion.SetSize( currentSize );
  validRegion.SetIndex ( inputPtr->GetBufferedRegion().GetIndex() );

  // The first time through the loop our input image is inputPtr
  typename TInputImage::ConstPointer workingImage;
  workingImage = inputPtr;

  unsigned int count = scratchRegion.GetNumberOfPixels() * ImageDimension;
  ProgressReporter progress(this,0,count,10);
  for (unsigned int n=0; n < ImageDimension; n++)
    {
    // Setup iterators for input image.
    ConstInputImageIterator inIterator1( workingImage, validRegion );
    ConstOutputImageIterator inIterator2( scratchImage, scratchRegion );

    if (n==0)
      {
      // First time through the loop we use the InputImage
      inIterator1.GoToBegin();
      inIterator1.SetDirection( n);
      }
    else
      {
      // After first time through the loop we use the scratch image which is of 
      // the output type.
      inIterator2.GoToBegin();
      inIterator2.SetDirection( n);
      }

    
    // Setup iterators and bounds for output image.
    currentSize[n] = currentSize[n]/2;  // reduce by a factor of 2
    validRegion.SetSize(currentSize); 
    // TODO:  Is there a way to put this in the else statement below? 
    OutputImageIterator outIterator( scratchImage, validRegion ); 
    if (n == ( ImageDimension - 1) )
      {
      // Last time through the loop write directly to the ouput
      outIterator = outItr; 
      }
       
    outIterator.GoToBegin();
    outIterator.SetDirection(n);
    
    if (n==0)
      {
      while (!inIterator1.IsAtEnd() )
        {
        // Copies one line of input to m_Scratch
        this->CopyInputLineToScratch( inIterator1 );

        this->Reduce1DImage(  m_Scratch, outIterator,  startSize[n], progress );
        inIterator1.NextLine();
        outIterator.NextLine();
        }
      }
    else
      {
      while (!inIterator2.IsAtEnd() )
        {
        // Copies one line of input to m_Scratch
        this->CopyOutputLineToScratch( inIterator2 );

        this->Reduce1DImage(  m_Scratch, outIterator,  startSize[n], progress );
        inIterator2.NextLine();
        outIterator.NextLine();
        }
      }
    
    // After first loop the input image is scratchImage
    //workingImage = scratchImage;
    
    }

}

template < class TInputImage, class TOutputImage>
void BSplineResampleImageFilterBase<TInputImage, TOutputImage>
::ExpandNDImage(OutputImageIterator &outItr)
{
  // Set up variables for waking the image regions.
  RegionType validRegion;
  SizeType startSize;
  SizeType currentSize;

  // Does not support streaming
  typename Superclass::InputImagePointer  inputPtr = const_cast< TInputImage * > ( this->GetInput() );
  startSize = inputPtr->GetBufferedRegion().GetSize();
  

  // Initilize scratchImage space and allocate memory
  InitializeScratch(startSize);
  typename TOutputImage::Pointer scratchImage;
  scratchImage =  TOutputImage::New();
  scratchImage->CopyInformation( inputPtr );
  RegionType scratchRegion;
  scratchRegion = inputPtr->GetBufferedRegion();
  currentSize = startSize;

  // scratchImage 2 times the space of the original image .

  for (unsigned int n = 0; n < ImageDimension; n++)
    {
    currentSize[n] = currentSize[n] * 2;  
    }
  scratchRegion.SetSize( currentSize );
  scratchImage->SetBufferedRegion( scratchRegion );
  scratchImage->Allocate();  
  
 
  currentSize = startSize;
  validRegion.SetSize( currentSize );
  validRegion.SetIndex ( inputPtr->GetBufferedRegion().GetIndex() );

  // The first time through the loop our input image is m_Image
  typename TInputImage::ConstPointer workingImage;
  workingImage = inputPtr;

  RegionType workingRegion = validRegion;

  unsigned int count = scratchRegion.GetNumberOfPixels() * ImageDimension;
  ProgressReporter progress(this,0,count,10);
  for (unsigned int n=0; n < ImageDimension; n++)
    {
    // Setup iterators for input image.
    ConstInputImageIterator inIterator1( workingImage, workingRegion);
    ConstOutputImageIterator inIterator2( scratchImage, validRegion);
    if (n==0)
      {
      // First time through the loop we use the InputImage
      inIterator1.GoToBegin();
      inIterator1.SetDirection( n);
      }
    else
      {
      // After first time through the loop we use the scratch image which is of 
      // the output type.
      inIterator2.GoToBegin();
      inIterator2.SetDirection( n);
      }
    
    // Setup iterators and bounds for output image.
    currentSize[n] = currentSize[n] * 2;  // expand by a factor of 2
    validRegion.SetSize(currentSize); 

    OutputImageIterator outIterator( scratchImage, validRegion );
    if (n == ( ImageDimension - 1) )
      {
      // Last time through the loop write directly to the ouput
      outIterator = outItr; 
      }


    outIterator.GoToBegin();
    outIterator.SetDirection(n);
    
    if (n==0)
      {
      while (!inIterator1.IsAtEnd() )
        {
        // Copies one line of input to m_Scratch
        this->CopyInputLineToScratch( inIterator1 );

        this->Expand1DImage(  m_Scratch, outIterator,  startSize[n], progress );
        inIterator1.NextLine();
        outIterator.NextLine();
        }
      }
    else
      {
      while (!inIterator2.IsAtEnd() )
        {
        // Copies one line of input to m_Scratch
        this->CopyOutputLineToScratch( inIterator2 );

        this->Expand1DImage(  m_Scratch, outIterator,  startSize[n], progress );
        inIterator2.NextLine();
        outIterator.NextLine();
        }
      }
      
    
    }
}

// Allocate scratch space
template <class TInputImage, class TOutputImage>
void BSplineResampleImageFilterBase<TInputImage, TOutputImage>
::InitializeScratch(SizeType DataLength)
{
  unsigned int maxLength = 0;
  for ( unsigned int n = 0; n < ImageDimension; n++ )
    {
    if ( DataLength[n] > maxLength )
      {
      maxLength = DataLength[n];
      }
    }
  m_Scratch.resize( maxLength );
}

template <class TInputImage, class TOutputImage>
void BSplineResampleImageFilterBase<TInputImage, TOutputImage>
::CopyLineToScratch(ConstInputImageIterator & Iter)
{
  unsigned int j = 0;
  while ( !Iter.IsAtEndOfLine() )
    {
    m_Scratch[j] = static_cast<double>( Iter.Get() );
    ++Iter;
    ++j;
    }
}


template <class TInputImage, class TOutputImage>
void BSplineResampleImageFilterBase<TInputImage, TOutputImage>
::CopyInputLineToScratch(ConstInputImageIterator & Iter)
{
  unsigned int j = 0;
  while ( !Iter.IsAtEndOfLine() )
    {
    m_Scratch[j] = static_cast<double>( Iter.Get() );
    ++Iter;
    ++j;
    }
}

template <class TInputImage, class TOutputImage>
void BSplineResampleImageFilterBase<TInputImage, TOutputImage>
::CopyOutputLineToScratch(ConstOutputImageIterator & Iter)
{
  unsigned int j = 0;
  while ( !Iter.IsAtEndOfLine() )
    {
    m_Scratch[j] = static_cast<double>( Iter.Get() );
    ++Iter;
    ++j;
    }
}


} // namespace itk

#endif