[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