SubsampleVolume-alex.cxx

// Alex Chekholko 2006-01-19
// based on SubsampleVolume.cxx

// should simply interpolate the input image with the given parameters
// using the WindowedSincInterpolateImageFunction with radius 5 and Welch window 

#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif

#ifdef __BORLANDC__
#define ITK_LEAN_AND_MEAN
#endif

#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkCastImageFilter.h"
#include "itkResampleImageFilter.h"
#include "itkIdentityTransform.h"
#include "itkWindowedSincInterpolateImageFunction.h"
#include "itkConstantBoundaryCondition.h"

int main( int argc, char * argv[] )
{
  if( argc < 6 )
    {
    std::cerr << "Usage: " << std::endl;
    std::cerr << argv[0]
      << "  inputImageFile  outputImageFile factorX factorY factorZ"
      << std::endl;
    return EXIT_FAILURE;
    }

  const     unsigned int    Dimension = 3;
  typedef   unsigned char   InputPixelType;
  typedef   float           InternalPixelType;
  typedef   unsigned char   OutputPixelType;

  typedef itk::Image< InputPixelType,    Dimension >   InputImageType;
  typedef itk::Image< InternalPixelType, Dimension >   InternalImageType;
  typedef itk::Image< OutputPixelType,   Dimension >   OutputImageType;

// read in the file
  typedef itk::ImageFileReader< InputImageType  >  ReaderType;
  ReaderType::Pointer reader = ReaderType::New();
  reader->SetFileName( argv[1] );

  const double factorX = atof( argv[3] );
  const double factorY = atof( argv[4] );
  const double factorZ = atof( argv[5] );

  try
    {
    reader->Update();
    }
  catch( itk::ExceptionObject & excep )
    {
    std::cerr << "Exception caught while reading input file!" << std::endl;
    std::cerr << excep << std::endl;
    }

  InputImageType::ConstPointer inputImage = reader->GetOutput();

  const InputImageType::SpacingType& inputSpacing = inputImage->GetSpacing();

// A casting filter is instantiated in order to convert the pixel type of the
// input image into the pixel type desired for computing the resampling.

  typedef itk::CastImageFilter< InputImageType,
                                InternalImageType >   CastFilterType;

  CastFilterType::Pointer caster = CastFilterType::New();
  caster->SetInput( inputImage );

  typedef itk::ResampleImageFilter<
                  InternalImageType, OutputImageType >  ResampleFilterType;

  ResampleFilterType::Pointer resampler = ResampleFilterType::New();

// Since the resampling is performed in the same physical extent of the input
// image, we select the IdentityTransform as the one to be used by the resampling
// filter.
  typedef itk::IdentityTransform< double, Dimension >  TransformType;

  TransformType::Pointer transform = TransformType::New();
  transform->SetIdentity();
  resampler->SetTransform( transform );

  const unsigned int WindowRadius = 5;
  typedef itk::Function::WelchWindowFunction< WindowRadius > WindowFunctionType;

// by default, the constant = 0 
  typedef itk::ConstantBoundaryCondition< InternalImageType > BoundaryConditionType;

  typedef itk::WindowedSincInterpolateImageFunction<InternalImageType, WindowRadius, 
		WindowFunctionType, BoundaryConditionType, double>  InterpolatorType;

  InterpolatorType::Pointer interpolator = InterpolatorType::New();

  resampler->SetInterpolator( interpolator );
  resampler->SetDefaultPixelValue( 0 ); // value for regions without source

  OutputImageType::SpacingType spacing;

  spacing[0] = inputSpacing[0] * factorX;
  spacing[1] = inputSpacing[1] * factorY;
  spacing[2] = inputSpacing[2] * factorZ;

  resampler->SetOutputSpacing( spacing );

// The origin of the input image is preserved and passed to the output image.
  resampler->SetOutputOrigin( inputImage->GetOrigin() );

// The number of pixels to use along each direction on the grid of the
// resampled image is computed using the number of pixels in the input image
// and the sampling factors.
  InputImageType::SizeType   inputSize =
              inputImage->GetLargestPossibleRegion().GetSize();

  typedef InputImageType::SizeType::SizeValueType SizeValueType;

  InputImageType::SizeType   size;

  size[0] = static_cast< SizeValueType >( inputSize[0] / factorX );
  size[1] = static_cast< SizeValueType >( inputSize[1] / factorY );
  size[2] = static_cast< SizeValueType >( inputSize[2] / factorZ );

  resampler->SetSize( size );

//connect to the pipeline
  resampler->SetInput( caster->GetOutput() );


  typedef itk::ImageFileWriter< OutputImageType >  WriterType;

  WriterType::Pointer writer = WriterType::New();

  writer->SetInput( resampler->GetOutput() );

  writer->SetFileName( argv[2] );

  try
    {
    writer->Update();
    }
  catch( itk::ExceptionObject & excep )
    {
    std::cerr << "Exception caught!" << std::endl;
    std::cerr << excep << std::endl;
    }

//  std::cout << "Resampling Done !" << std::endl;


  return EXIT_SUCCESS;
}