[ITK] [ITK-users] ResampleDICOM example and DCMTK

Nick Patterson pattersonnp.work at gmail.com
Wed Aug 27 12:28:19 EDT 2014


Hi,
 
I have been trying to test the ResampledDICOM ITK example using DCMTK instead of GDCM. I have attached the code below. The code compiles and runs, but does not output the resultant files - which is strange, as all of the expected parameters are good (i.e. number of files, origin of image, image dimensions etc).

Can anyone possibly comment on any possible reason why the code execute without error, but does not result in any output files?

Regards, Nick



#include "itkVersion.h"
#include "itkImage.h"
#include "itkMinimumMaximumImageFilter.h"
#include "itkDCMTKImageIO.h"
#include "itkDCMTKSeriesFileNames.h"
#include "itkNumericSeriesFileNames.h"
#include "itkImageSeriesReader.h"
#include "itkImageSeriesWriter.h"
#include "itkResampleImageFilter.h"
#include "itkShiftScaleImageFilter.h"
#include "itkIdentityTransform.h"
#include "itkLinearInterpolateImageFunction.h"
#include <itksys/SystemTools.hxx>
#include "dcmdata/dcuid.h"
#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::DCMTKImageIO ImageIOType;
  typedef itk::DCMTKSeriesFileNames 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.

  char seriesUID[100];
  dcmGenerateUniqueIdentifier(seriesUID);

  char frameOfReferenceUID[100];
  dcmGenerateUniqueIdentifier(frameOfReferenceUID);

  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);

    char sopInstanceUID[100];
    dcmGenerateUniqueIdentifier(sopInstanceUID);

    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]);
  std::cout << "Output Size: " << outputSize[2] << std::endl;
  std::cout << outputNames->GetFileNames().at(0) << std::endl;

  SeriesWriterType::Pointer seriesWriter = SeriesWriterType::New();
  seriesWriter->SetInput( shiftScale->GetOutput() );
  seriesWriter->SetImageIO( gdcmIO );
  seriesWriter->SetFileNames( outputNames->GetFileNames() );
  seriesWriter->SetMetaDataDictionaryArray( &outputArray );

  std::cout << seriesWriter->GetImageIO()->GetOrigin(0) << " " << seriesWriter->GetImageIO()->GetOrigin(1) <<
               " " << seriesWriter->GetImageIO()->GetOrigin(2) << std::endl;

  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;
    }
}

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/community/attachments/20140827/350a059f/attachment-0001.html>
-------------- next part --------------
_____________________________________
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.php

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://public.kitware.com/mailman/listinfo/insight-users


More information about the Community mailing list