[Insight-users] translating image volume using affine transformation

[Bill Lorensen]

What is the origin, spacing and direction of the input volume?

On Fri, Nov 27, 2009 at 9:59 AM, michiel mentink
<michael.mentink at st-hughs.ox.ac.uk> wrote:
> Dear all,
>
> I'm trying to make a simple function that translates a dicom volume by a
> certain amount.
> It first opens all dicom slides of an mri scan.
>
> To do this, I pasted DicomSeriesReadImageWrite2.cxx and
> ResampleImageFilter.cxx and adjusted
> the latter so it accepts 3 dimensions instead of two and uses integer
> instead of unsigned char.
>
> The function compiles and a .vtk volume is written. However, there's no
> image data in the file.
> It does display all pixels with value 100, as I set by:
> filter->SetDefaultPixelValue( 100 );
>
> questions:
> 1) am I too naive to think that the ResampleImageFilter can easily be
> extended from 2 to 3 dimensions?
> 2) if not, what did I do wrong?
>
> I pasted the source code below.
>
> //////////////////////////////////////////////////////////////////////////////////////////////////
> #if defined(_MSC_VER)
> #pragma warning ( disable : 4786 )
> #endif
>
> #ifdef __BORLANDC__
> #define ITK_LEAN_AND_MEAN
> #endif
>
> #include "itkOrientedImage.h"
> #include "itkGDCMImageIO.h"
> #include "itkGDCMSeriesFileNames.h"
> #include "itkImageSeriesReader.h"
> #include "itkImageFileWriter.h"
> #include "itkResampleImageFilter.h"
>
> #include "itkAffineTransform.h"
> #include "itkNearestNeighborInterpolateImageFunction.h"
> // Software Guide : EndCodeSnippet
>
> int main( int argc, char* argv[] )
> {
>
>   if( argc < 3 )
>     {
>     std::cerr << "Usage: " << std::endl;
>     std::cerr << argv[0] << " DicomDirectory  outputFileName  [seriesName]"
>               << std::endl;
>     return EXIT_FAILURE;
>     }
>
>   typedef signed short    PixelType;
>   const unsigned int      Dimension = 3;
>
>   typedef itk::OrientedImage< PixelType, Dimension >         ImageType;
>
>   typedef itk::ImageSeriesReader< ImageType >        ReaderType;
>   ReaderType::Pointer reader = ReaderType::New();
>
>   typedef itk::GDCMImageIO       ImageIOType;
>   ImageIOType::Pointer dicomIO = ImageIOType::New();
>
>   reader->SetImageIO( dicomIO );
>
>   typedef itk::GDCMSeriesFileNames NamesGeneratorType;
>   NamesGeneratorType::Pointer nameGenerator = NamesGeneratorType::New();
>
>   nameGenerator->SetUseSeriesDetails( true );
>   nameGenerator->AddSeriesRestriction("0008|0021" );
>
>   nameGenerator->SetDirectory( argv[1] );
>
>   try
>     {
>     std::cout << std::endl << "The directory: " << std::endl;
>     std::cout << std::endl << argv[1] << std::endl << std::endl;
>     std::cout << "Contains the following DICOM Series: ";
>     std::cout << std::endl << std::endl;
>
>     typedef std::vector< std::string >    SeriesIdContainer;
>
>     const SeriesIdContainer & seriesUID = nameGenerator->GetSeriesUIDs();
>
>     SeriesIdContainer::const_iterator seriesItr = seriesUID.begin();
>     SeriesIdContainer::const_iterator seriesEnd = seriesUID.end();
>     while( seriesItr != seriesEnd )
>       {
>       std::cout << seriesItr->c_str() << std::endl;
>       seriesItr++;
>       }
>
>     std::string seriesIdentifier;
>
>     if( argc > 3 ) // If no optional series identifier
>       {
>       seriesIdentifier = argv[3];
>       }
>     else
>       {
>       seriesIdentifier = seriesUID.begin()->c_str();
>       }
>
>     std::cout << std::endl << std::endl;
>     std::cout << "Now reading series: " << std::endl << std::endl;
>     std::cout << seriesIdentifier << std::endl;
>     std::cout << std::endl << std::endl;
>
>     typedef std::vector< std::string >   FileNamesContainer;
>     FileNamesContainer fileNames;
>
>     fileNames = nameGenerator->GetFileNames( seriesIdentifier );
>     reader->SetFileNames( fileNames );
>     try
>       {
>       reader->Update();
>       }
>     catch (itk::ExceptionObject &ex)
>       {
>       std::cout << ex << std::endl;
>       return EXIT_FAILURE;
>       }
>     }
>   catch (itk::ExceptionObject &ex)
>     {
>     std::cout << ex << std::endl;
>     return EXIT_FAILURE;
>     }
> //  return EXIT_SUCCESS;
>
> //////////////////////////////////////////////////////////////////////////////////////
> // FROM HERE ADJUSTED ResampleImageFilter.cxx CODE
> /////////////////////////////////////////////////////////////////////////////////////
>
> //  const     unsigned int   Dimension = 3;
>   typedef   unsigned char  InputPixelType;
>   typedef   unsigned char  OutputPixelType;
>   typedef itk::Image< InputPixelType,  Dimension >   InputImageType;
> //  typedef itk::Image< OutputPixelType, Dimension >   OutputImageType;
>
>   typedef itk::Image< PixelType, Dimension >   OutputImageType;
>   // typedef itk::ImageFileReader< InputImageType  >  ReaderType;
>
>   typedef itk::ImageFileWriter< ImageType >  WriterType;
>   WriterType::Pointer writer = WriterType::New();
>
>
>   writer->SetFileName( argv[2] );
>
> //  typedef itk::ResampleImageFilter<ImageType,OutputImageType> FilterType;
>   typedef itk::ResampleImageFilter<ImageType,ImageType> FilterType;
>   FilterType::Pointer filter = FilterType::New();
>
>   typedef itk::AffineTransform< double, Dimension >  TransformType;
>   TransformType::Pointer transform = TransformType::New();
>   filter->SetTransform( transform );
>
>
>   typedef itk::NearestNeighborInterpolateImageFunction<ImageType, double >
> InterpolatorType;
>   InterpolatorType::Pointer interpolator = InterpolatorType::New();
>   filter->SetInterpolator( interpolator );
>
>
>   filter->SetDefaultPixelValue( 100 ); // if pixel in output image not
> defined, what value?
>
>
>   double spacing[ Dimension ];
>   spacing[0] = 0.29; // pixel spacing in millimeters along X
>   spacing[1] = 0.29; // pixel spacing in millimeters along Y
>   spacing[2] = 3.3;
>
>   filter->SetOutputSpacing( spacing );
>
>   double origin[ Dimension ];
>   origin[0] = 0.0;  // X space coordinate of origin
>   origin[1] = 0.0;  // Y space coordinate of origin
>   origin[2] = 0.0;
>
>   filter->SetOutputOrigin( origin );
>   ImageType::DirectionType direction;
>   direction.SetIdentity();
>   filter->SetOutputDirection( direction );
>   ImageType::SizeType   size;
>
>   size[0] = 560;  // number of pixels along X
>   size[1] = 560;  // number of pixels along Y
>   size[2] = 300;
>
>   filter->SetSize( size );
>   filter->SetInput( reader->GetOutput() );
> //  writer->SetInput( filter->GetOutput() );
> //  writer->Update();
>   TransformType::OutputVectorType translation;
>   translation[0] = 0;  // X translation in millimeters
>   translation[1] = 0;  // Y translation in millimeters
>   translation[2] = 0;
>
>   transform->Translate( translation );
>   // Software Guide : EndCodeSnippet
>
>    writer->SetInput( filter->GetOutput() );
>
>     std::cout  << "Writing the image as " << std::endl << std::endl;
>     std::cout  << argv[2] << std::endl << std::endl;
>     writer->Update();
>
> }
>
> _____________________________________
> Powered by www.kitware.com
>
> Visit other Kitware open-source projects at
> http://www.kitware.com/opensource/opensource.html
>
> Kitware offers ITK Training Courses, for more information visit:
> http://www.kitware.com/products/protraining.html
>
> Please keep messages on-topic and check the ITK FAQ at:
> http://www.itk.org/Wiki/ITK_FAQ
>
> Follow this link to subscribe/unsubscribe:
> http://www.itk.org/mailman/listinfo/insight-users
>
>