[Insight-users] Unable to write dicom tag using GDCMImageIO

Bill Lorensen bill.lorensen at gmail.com
Sun Oct 7 12:44:36 EDT 2012


Please keep replies on the list...

Why are you resampling? I thought you wanted to change the directions.

On Sun, Oct 7, 2012 at 10:15 AM, suicheng gu <gusuch at yahoo.com> wrote:

> attached, I just modified the ResampleDICOM.cxx file
>
>
>
> // 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
> //
>
> #include "itkVersion.h"
>
> #include "itkImage.h"
> #include "itkMinimumMaximumImageFilter.h"
>
> #include "itkGDCMImageIO.h"
> #include "itkGDCMSeriesFileNames.h"
> #include "itkNumericSeriesFileNames.h"
>
> #include "itkImageSeriesReader.h"
> #include "itkImageSeriesWriter.h"
>
> #include "itkResampleImageFilter.h"
> #include "itkShiftScaleImageFilter.h"
> #include "itkMetaDataObject.h"
> #include "itkIdentityTransform.h"
> #include "itkLinearInterpolateImageFunction.h"
> #include "itkOrientImageFilter.h"
> #include <itksys/SystemTools.hxx>
>
> #if ITK_VERSION_MAJOR >= 4
> #include "gdcmUIDGenerator.h"
> #else
> #include "gdcm/src/gdcmFile.h"
> #include "gdcm/src/gdcmUtil.h"
> #endif
>
> #include <string>
>
> static void CopyDictionary (itk::MetaDataDictionary &fromDict,
>                      itk::MetaDataDictionary &toDict);
>
>
> int main( int argc, char* argv[] )
> {
>   // Validate input parameters
>   if( argc < 2 )
>     {
>     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();
>   const double
> inputOrigin[3]={reader->GetOutput()->GetOrigin()[0],reader->GetOutput()->GetOrigin()[1],reader->GetOutput()->GetOrigin()[2]};
>  // const double outputOrigin[3] =
> {inputOrigin[0],inputOrigin[1],inputOrigin[2]};
>
>   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] = inputSpacing[0];
>   outputSpacing[1] = inputSpacing[2];
>   outputSpacing[2] = inputSpacing[1];
>
>   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] = inputSize[0];
>   outputSize[1] = inputSize[2];
>   outputSize[2] = inputSize[1];
>
>
>
> ////////////////////////////////////////////////
> // 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.
> #if ITK_VERSION_MAJOR >= 4
>   gdcm::UIDGenerator suid;
>   std::string seriesUID = suid.Generate();
>   gdcm::UIDGenerator fuid;
>   std::string frameOfReferenceUID = fuid.Generate();
> #else
>   std::string seriesUID = gdcm::Util::CreateUniqueUID(
> gdcmIO->GetUIDPrefix());
>   std::string frameOfReferenceUID = gdcm::Util::CreateUniqueUID(
> gdcmIO->GetUIDPrefix());
> #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);
>
>  #if ITK_VERSION_MAJOR >= 4
>   gdcm::UIDGenerator sopuid;
>   std::string sopInstanceUID = sopuid.Generate();
>  #else
>   std::string sopInstanceUID = gdcm::Util::CreateUniqueUID(
> gdcmIO->GetUIDPrefix());
>  #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 orientation
>   value.str("");
>   value<<"1\\0\\0\\0\\0\\1";
>   itk::EncapsulateMetaData<std::string>(*dict,"0020|0037", value.str());
>   std::string orientation;
>   itk::ExposeMetaData<std::string>(*dict, "0020|0037", orientation);
>   //std::cout<<orientation<<std::endl;
>   // Image Position Patient: This is calculated by computing the
>   // physical coordinate of the first pixel in each slice.
>   value.str("");
>   value << inputOrigin[0] << "\\"
> <<(inputOrigin[1]+outputSpacing[2]*f)<<"\\"<< inputOrigin[2] ;
>   //std::cout<<"position: "<<value.str();
>   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 << inputOrigin[1]+outputSpacing[2]*f;
>   itk::EncapsulateMetaData<std::string>(*dict,"0020|1041",
> value.str());
>
>   // 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
>  /*InputImageType::Pointer image = InputImageType::New();
>  image->SetSpacing( outputSpacing );
>  InputImageType::IndexType start;
>  start[0] =   0;
>  start[1] =   0;
>  start[2] =   0;
>  InputImageType::RegionType region;
>  region.SetSize( outputSize );
>  region.SetIndex( start );
>  image->SetRegions( region );
>  image->SetOrigin(inputOrigin);
>
>  image->Allocate();
>  //image->FillBuffer(0);
>
>  PixelType *ppt = image->GetBufferPointer();
>  for(int z = 0;z<outputSize[2];z++)
>  {
>   for(int y = 0;y<outputSize[1];y++)
>   {
>    for(int x = 0;x<outputSize[0];x++)
>    {
>     InputImageType::IndexType pixelIndex;
>     pixelIndex[0] = x;
>     pixelIndex[1] = z;
>     pixelIndex[2] = y;
>     PixelType pt = reader->GetOutput()->GetPixel(pixelIndex);
>     *ppt = pt;
>     ppt ++;
>    }
>   }
>  }*/
>   InputImageType::DirectionType direction;
>    direction(0,0) = 1;
>    direction(0,1) = 0;
>    direction(0,2) = 0;
>    direction(1,0) = 0;
>    direction(1,1) = 1;
>    direction(1,2) = 0;
>    direction(2,0) = 0;
>    direction(2,1) = 0;
>    direction(2,2) = 1;
>
>  //image->SetDirection(direction);
>  //image->Update();
>
>  /*
>
> itk::OrientImageFilter<InputImageType,InputImageType>::Pointer orienter =
>
> itk::OrientImageFilter<InputImageType,InputImageType>::New();
>
> orienter->UseImageDirectionOn();
> orienter->SetDesiredCoordinateOrientationToCoronal();
> orienter->SetInput(reader->GetOutput());
>
> orienter->Update();
> */
> itk::SpatialOrientation::ValidCoordinateOrientationFlags fileOrientation;
>
>
> //itk::ExposeMetaData<itk::SpatialOrientation::ValidCoordinateOrientationFlags>(reader->GetOutput()->GetMetaDataDictionary(),itk::ITK_CoordinateOrientation,fileOrientation);
>
>
> itk::OrientImageFilter<InputImageType,InputImageType>::Pointer orienter =
>
> itk::OrientImageFilter<InputImageType,InputImageType>::New();
>
> orienter->SetGivenCoordinateOrientation(itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_LPS);
> // deprecated
>
>
> orienter->SetDesiredCoordinateOrientation(itk::SpatialOrientation::ITK_COORDINATE_ORIENTATION_RIP);
>
>
> orienter->SetInput(reader->GetOutput());
>
> orienter->Update();
>
> ////////////////////////////////////////////////
> // 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]);
>
>   //itk::MetaDataDictionary & dict1 = gdcmIO->GetMetaDataDictionary();
>   //itk::EncapsulateMetaData<std::string>(dict1,"0020|0037",
> "1\\0\\0\\0\\0\\1");
>
>  SeriesWriterType::Pointer seriesWriter = SeriesWriterType::New();
>  seriesWriter->SetInput( orienter->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;
>  system("pause");
>   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;
>     }
> }
>
>    *From:* Bill Lorensen <bill.lorensen at gmail.com>
> *To:* suicheng gu <gusuch at yahoo.com>
> *Cc:* "insight-users at itk.org" <insight-users at itk.org>
> *Sent:* Sunday, October 7, 2012 8:14 PM
>
> *Subject:* Re: [Insight-users] Unable to write dicom tag using GDCMImageIO
>
> Can you post the entire program please?
>
> On Sat, Oct 6, 2012 at 9:04 PM, suicheng gu <gusuch at yahoo.com> wrote:
>
> I used the following code to obtain the coronal view, but the orientation
> tag is still 1\0\0\0\1\0
> And one more problem is that some dcm files are with size of 0.
>
> itk::OrientImageFilter<InputImageType,InputImageType>::Pointer orienter =
> itk::OrientImageFilter<InputImageType,InputImageType>::New();
> orienter->UseImageDirectionOn();
> orienter->SetDesiredCoordinateOrientationToCoronal();
> orienter->SetInput(reader->GetOutput());
> orienter->Update();
>
>    *From:* Bill Lorensen <bill.lorensen at gmail.com>
> *To:* suicheng gu <gusuch at yahoo.com>
> *Cc:* "insight-users at itk.org" <insight-users at itk.org>
> *Sent:* Sunday, October 7, 2012 4:33 AM
>
> *Subject:* Re: [Insight-users] Unable to write dicom tag using GDCMImageIO
>
> You can use itkOrientImageFilter to change the direction of the volume.
> Then when you write the IDCOM with itk, the direction cosines will be
> correct.
>
> On Sat, Oct 6, 2012 at 11:24 AM, suicheng gu <gusuch at yahoo.com> wrote:
>
> Thanks very much.
> what I want to do is to reconstruct the coronal and saggittal views and
> save them as dcm files.
> I can manually change the direction of the volume data
> If I can't change these tags, can I just create one dictionary? such as
> the example: ImageReadDicomSeriesWrite.cxx
>
>     *From:* Bill Lorensen <bill.lorensen at gmail.com>
> *To:* suicheng gu <gusuch at yahoo.com>
> *Cc:* "insight-users at itk.org" <insight-users at itk.org>
> *Sent:* Saturday, October 6, 2012 11:07 PM
> *Subject:* Re: [Insight-users] Unable to write dicom tag using GDCMImageIO
>
> If I recall, tags that are related to orientation and spacing cannot be
> changed. You must change the orientation of the data before you write the
> dicom file. The writer will output the proper tag that goes with the data.
>
> On Sat, Oct 6, 2012 at 5:21 AM, suicheng gu <gusuch at yahoo.com> wrote:
>
>  Hi,
>
> I was trying to chage the image orientation tag using
> itk::EncapsulateMetaData<std::string>(dict1,"0020|0037",
> "1\\0\\0\\0\\0\\1");
> But the image orientation was still 1\0\0\0\1\0
>
>
>
> I also tested the provided exmaple:
> InsightToolkit-4.1.0\Examples\IO\DicomImageReadChangeHeaderWrite.cxx
>
> It works for some tags, such as the instance number:
> DicomImageReadChangeHeaderWrite $infile$ $outfile$ 0020|0013 10
>
> But it failed on the image orientation: DicomImageReadChangeHeaderWrite
> $infile$ $outfile$  0020|0037 1\0\0\0\0\1
>
>
>
> Anyone can help?
> Thanks in advance.
> Suicheng
>
>
> _____________________________________
> Powered by http://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://www.itk.org/mailman/listinfo/insight-users
>
>
>
>
> --
> Unpaid intern in BillsBasement at noware dot com
>
>
>
>
>
>
> --
> Unpaid intern in BillsBasement at noware dot com
>
>
>
>
>
>
> --
> Unpaid intern in BillsBasement at noware dot com
>
>
>
>


-- 
Unpaid intern in BillsBasement at noware dot com
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20121007/8d09a1aa/attachment.htm>


More information about the Insight-users mailing list