ITK/Examples/DICOM/ResampleDICOM

ResampleDICOM.cxx

<source lang="cpp"> // Resample a DICOM study // Usage: ResampleDICOM InputDirectory OutputDirectory // xSpacing ySpacing zSpacing // // Example: ResampleDICOM CT CTResample 0 0 1.5 // will read a series from the CT directory and create a // new series in the CTResample directory. The new series // will have the same x,y spacing as the input series, but // will have a z-spacing of 1.5. // // Description: // ResampleDICOM resamples a DICOM series with user-specified // spacing. The program outputs a new DICOM series with a series // number set to 1001. All non-private DICOM tags are moved from the input // series to the output series. The Image Position Patient is adjusted // for each slice to reflect the z-spacing. The number of slices in // the output series may be larger or smaller due to changes in the // z-spacing. To retain the spacing for a given dimension, specify 0. // // The program progresses as follows: // 1) Read the input series // 2) Resample the series according to the user specified x-y-z // spacing. // 3) Create a MetaDataDictionary for each slice. // 4) Shift data to undo the effect of a rescale intercept by the // DICOM reader // 5) Write the new DICOM series //

1. include "itkVersion.h"

1. include "itkImage.h"
2. include "itkMinimumMaximumImageFilter.h"

1. include "itkGDCMImageIO.h"
2. include "itkGDCMSeriesFileNames.h"
3. include "itkNumericSeriesFileNames.h"

1. include "itkImageSeriesReader.h"
2. include "itkImageSeriesWriter.h"

1. include "itkResampleImageFilter.h"
2. include "itkShiftScaleImageFilter.h"

1. include "itkIdentityTransform.h"
2. include "itkLinearInterpolateImageFunction.h"

1. include <itksys/SystemTools.hxx>

1. if ITK_VERSION_MAJOR >= 4
2. include "gdcmUIDGenerator.h"
3. else
4. include "gdcm/src/gdcmFile.h"
5. include "gdcm/src/gdcmUtil.h"
6. endif

1. include <string>

static void CopyDictionary (itk::MetaDataDictionary &fromDict,

                    itk::MetaDataDictionary &toDict);

int main( int argc, char* argv[] ) {

 // Validate input parameters
 if( argc < 4 )
   {
   std::cerr << "Usage: " 
             << argv[0]
             << " InputDicomDirectory OutputDicomDirectory spacing_x spacing_y spacing_z"
             << std::endl;
   return EXIT_FAILURE;
   }

 const unsigned int InputDimension = 3;
 const unsigned int OutputDimension = 2;

 typedef signed short PixelType;

 typedef itk::Image< PixelType, InputDimension >
   InputImageType;
 typedef itk::Image< PixelType, OutputDimension >
   OutputImageType;
 typedef itk::ImageSeriesReader< InputImageType >
   ReaderType;
 typedef itk::GDCMImageIO
   ImageIOType;
 typedef itk::GDCMSeriesFileNames
   InputNamesGeneratorType;
 typedef itk::NumericSeriesFileNames
   OutputNamesGeneratorType;
 typedef itk::IdentityTransform< double, InputDimension >
   TransformType;
 typedef itk::LinearInterpolateImageFunction< InputImageType, double >
   InterpolatorType;
 typedef itk::ResampleImageFilter< InputImageType, InputImageType >
   ResampleFilterType;
 typedef itk::ShiftScaleImageFilter< InputImageType, InputImageType >
   ShiftScaleType;
 typedef itk::ImageSeriesWriter< InputImageType, OutputImageType >
   SeriesWriterType;

//////////////////////////////////////////////// // 1) Read the input series

 ImageIOType::Pointer gdcmIO = ImageIOType::New();
 InputNamesGeneratorType::Pointer inputNames = InputNamesGeneratorType::New();
 inputNames->SetInputDirectory( argv[1] );

 const ReaderType::FileNamesContainer & filenames = 
                           inputNames->GetInputFileNames();

 ReaderType::Pointer reader = ReaderType::New();

 reader->SetImageIO( gdcmIO );
 reader->SetFileNames( filenames );
 try
   {
   reader->Update();
   }
 catch (itk::ExceptionObject &excp)
   {
   std::cerr << "Exception thrown while reading the series" << std::endl;
   std::cerr << excp << std::endl;
   return EXIT_FAILURE;
   }

//////////////////////////////////////////////// // 2) Resample the series

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

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

 const InputImageType::SpacingType& inputSpacing =
   reader->GetOutput()->GetSpacing();
 const InputImageType::RegionType& inputRegion =
   reader->GetOutput()->GetLargestPossibleRegion();
 const InputImageType::SizeType& inputSize =
   inputRegion.GetSize();

 std::cout << "The input series in directory " << argv[1]
           << " has " << filenames.size() << " files with spacing "
           << inputSpacing
           << std::endl;

 // Compute the size of the output. The user specifies a spacing on
 // the command line. If the spacing is 0, the input spacing will be
 // used. The size (# of pixels) in the output is recomputed using
 // the ratio of the input and output sizes.
 InputImageType::SpacingType outputSpacing;
 outputSpacing[0] = atof(argv[3]);
 outputSpacing[1] = atof(argv[4]);
 outputSpacing[2] = atof(argv[5]);

 bool changeInSpacing = false;
 for (unsigned int i = 0; i < 3; i++)
   {
   if (outputSpacing[i] == 0.0)
     {
     outputSpacing[i] = inputSpacing[i];
     }
   else
     {
     changeInSpacing = true;
     }
   }
 InputImageType::SizeType   outputSize;
 typedef InputImageType::SizeType::SizeValueType SizeValueType;
 outputSize[0] = static_cast<SizeValueType>(inputSize[0] * inputSpacing[0] / outputSpacing[0] + .5);
 outputSize[1] = static_cast<SizeValueType>(inputSize[1] * inputSpacing[1] / outputSpacing[1] + .5);
 outputSize[2] = static_cast<SizeValueType>(inputSize[2] * inputSpacing[2] / outputSpacing[2] + .5);

 ResampleFilterType::Pointer resampler = ResampleFilterType::New();
   resampler->SetInput( reader->GetOutput() );
   resampler->SetTransform( transform );
   resampler->SetInterpolator( interpolator );
   resampler->SetOutputOrigin ( reader->GetOutput()->GetOrigin());
   resampler->SetOutputSpacing ( outputSpacing );
   resampler->SetOutputDirection ( reader->GetOutput()->GetDirection());
   resampler->SetSize ( outputSize );
   resampler->Update ();

//////////////////////////////////////////////// // 3) Create a MetaDataDictionary for each slice.

 // Copy the dictionary from the first image and override slice
 // specific fields
 ReaderType::DictionaryRawPointer inputDict = (*(reader->GetMetaDataDictionaryArray()))[0];
 ReaderType::DictionaryArrayType outputArray;
   
 // To keep the new series in the same study as the original we need
 // to keep the same study UID. But we need new series and frame of
 // reference UID's.

1. if ITK_VERSION_MAJOR >= 4

 gdcm::UIDGenerator suid;
 std::string seriesUID = suid.Generate();
 gdcm::UIDGenerator fuid;
 std::string frameOfReferenceUID = fuid.Generate();

1. else

 std::string seriesUID = gdcm::Util::CreateUniqueUID( gdcmIO->GetUIDPrefix());
 std::string frameOfReferenceUID = gdcm::Util::CreateUniqueUID( gdcmIO->GetUIDPrefix());

1. endif

 std::string studyUID;
 std::string sopClassUID;
 itk::ExposeMetaData<std::string>(*inputDict, "0020|000d", studyUID);
 itk::ExposeMetaData<std::string>(*inputDict, "0008|0016", sopClassUID);
 gdcmIO->KeepOriginalUIDOn();

 for (unsigned int f = 0; f < outputSize[2]; f++)
   {
   // Create a new dictionary for this slice
   ReaderType::DictionaryRawPointer dict = new ReaderType::DictionaryType;

   // Copy the dictionary from the first slice
   CopyDictionary (*inputDict, *dict);

   // Set the UID's for the study, series, SOP  and frame of reference
   itk::EncapsulateMetaData<std::string>(*dict,"0020|000d", studyUID);
   itk::EncapsulateMetaData<std::string>(*dict,"0020|000e", seriesUID);
   itk::EncapsulateMetaData<std::string>(*dict,"0020|0052", frameOfReferenceUID);

1. if ITK_VERSION_MAJOR >= 4

   gdcm::UIDGenerator sopuid;
   std::string sopInstanceUID = sopuid.Generate();

1. else

   std::string sopInstanceUID = gdcm::Util::CreateUniqueUID( gdcmIO->GetUIDPrefix());

1. endif

   itk::EncapsulateMetaData<std::string>(*dict,"0008|0018", sopInstanceUID);
   itk::EncapsulateMetaData<std::string>(*dict,"0002|0003", sopInstanceUID);

   // Change fields that are slice specific
   itksys_ios::ostringstream value;
   value.str("");
   value << f + 1;

   // Image Number
   itk::EncapsulateMetaData<std::string>(*dict,"0020|0013", value.str());

   // Series Description - Append new description to current series
   // description
   std::string oldSeriesDesc;
   itk::ExposeMetaData<std::string>(*inputDict, "0008|103e", oldSeriesDesc);

   value.str("");
   value << oldSeriesDesc
         << ": Resampled with pixel spacing "
         << outputSpacing[0] << ", " 
         << outputSpacing[1] << ", " 
         << outputSpacing[2];
   // This is an long string and there is a 64 character limit in the 
   // standard
   unsigned lengthDesc = value.str().length();
   
   std::string seriesDesc( value.str(), 0,
                           lengthDesc > 64 ? 64
                           : lengthDesc);
   itk::EncapsulateMetaData<std::string>(*dict,"0008|103e", seriesDesc);

   // Series Number
   value.str("");
   value << 1001;
   itk::EncapsulateMetaData<std::string>(*dict,"0020|0011", value.str());

   // Derivation Description - How this image was derived
   value.str("");
   for (int i = 0; i < argc; i++)
     {
     value << argv[i] << " ";
     }
   value << ": " << ITK_SOURCE_VERSION;

   lengthDesc = value.str().length();
   std::string derivationDesc( value.str(), 0,
                               lengthDesc > 1024 ? 1024
                               : lengthDesc);
   itk::EncapsulateMetaData<std::string>(*dict,"0008|2111", derivationDesc);
   
   // Image Position Patient: This is calculated by computing the
   // physical coordinate of the first pixel in each slice.
   InputImageType::PointType position;
   InputImageType::IndexType index;
   index[0] = 0;
   index[1] = 0;
   index[2] = f;
   resampler->GetOutput()->TransformIndexToPhysicalPoint(index, position);

   value.str("");
   value << position[0] << "\\" << position[1] << "\\" << position[2];
   itk::EncapsulateMetaData<std::string>(*dict,"0020|0032", value.str());      
   // Slice Location: For now, we store the z component of the Image
   // Position Patient.
   value.str("");
   value << position[2];
   itk::EncapsulateMetaData<std::string>(*dict,"0020|1041", value.str());      

   if (changeInSpacing)
     {
     // Slice Thickness: For now, we store the z spacing
     value.str("");
     value << outputSpacing[2];
     itk::EncapsulateMetaData<std::string>(*dict,"0018|0050",
                                           value.str());
     // Spacing Between Slices
     itk::EncapsulateMetaData<std::string>(*dict,"0018|0088",
                                           value.str());
     }
     
   // Save the dictionary
   outputArray.push_back(dict);
   }
   

//////////////////////////////////////////////// // 4) Shift data to undo the effect of a rescale intercept by the // DICOM reader

 std::string interceptTag("0028|1052");
 typedef itk::MetaDataObject< std::string > MetaDataStringType;
 itk::MetaDataObjectBase::Pointer entry = (*inputDict)[interceptTag];
   
 MetaDataStringType::ConstPointer interceptValue = 
   dynamic_cast<const MetaDataStringType *>( entry.GetPointer() ) ;
   
 int interceptShift = 0;
 if( interceptValue )
   {
   std::string tagValue = interceptValue->GetMetaDataObjectValue();
   interceptShift = -atoi ( tagValue.c_str() );
   }

 ShiftScaleType::Pointer shiftScale = ShiftScaleType::New();
   shiftScale->SetInput( resampler->GetOutput());
   shiftScale->SetShift( interceptShift );

//////////////////////////////////////////////// // 5) Write the new DICOM series

 // Make the output directory and generate the file names.
 itksys::SystemTools::MakeDirectory( argv[2] );

 // Generate the file names
 OutputNamesGeneratorType::Pointer outputNames = OutputNamesGeneratorType::New();
 std::string seriesFormat(argv[2]);
 seriesFormat = seriesFormat + "/" + "IM%d.dcm";
 outputNames->SetSeriesFormat (seriesFormat.c_str());
 outputNames->SetStartIndex (1);
 outputNames->SetEndIndex (outputSize[2]);
 
 SeriesWriterType::Pointer seriesWriter = SeriesWriterType::New();
   seriesWriter->SetInput( shiftScale->GetOutput() );
   seriesWriter->SetImageIO( gdcmIO );
   seriesWriter->SetFileNames( outputNames->GetFileNames() );
   seriesWriter->SetMetaDataDictionaryArray( &outputArray );
 try
   {
   seriesWriter->Update();
   }
 catch( itk::ExceptionObject & excp )
   {
   std::cerr << "Exception thrown while writing the series " << std::endl;
   std::cerr << excp << std::endl;
   return EXIT_FAILURE;
   }
 std::cout << "The output series in directory " << argv[2]
           << " has " << outputSize[2] << " files with spacing "
           << outputSpacing
           << std::endl;
 return EXIT_SUCCESS;

}

void CopyDictionary (itk::MetaDataDictionary &fromDict, itk::MetaDataDictionary &toDict) {

 typedef itk::MetaDataDictionary DictionaryType;

 DictionaryType::ConstIterator itr = fromDict.Begin();
 DictionaryType::ConstIterator end = fromDict.End();
 typedef itk::MetaDataObject< std::string > MetaDataStringType;

 while( itr != end )
   {
   itk::MetaDataObjectBase::Pointer  entry = itr->second;

   MetaDataStringType::Pointer entryvalue = 
     dynamic_cast<MetaDataStringType *>( entry.GetPointer() ) ;
   if( entryvalue )
     {
     std::string tagkey   = itr->first;
     std::string tagvalue = entryvalue->GetMetaDataObjectValue();
     itk::EncapsulateMetaData<std::string>(toDict, tagkey, tagvalue); 
     }
   ++itr;
   }

} </source>

CMakeLists.txt

<source lang="cmake"> cmake_minimum_required(VERSION 2.6)

PROJECT(ResampleDICOM)

FIND_PACKAGE(ITK REQUIRED) INCLUDE(${ITK_USE_FILE})

ADD_EXECUTABLE(ResampleDICOM ResampleDICOM.cxx) TARGET_LINK_LIBRARIES(ResampleDICOM ${ITK_LIBRARIES}) </source>