[ITK-users] Resampling CT is fine... but not PET. Any clues?
Matt McCormick
matt.mccormick at kitware.com
Wed Apr 16 10:17:15 EDT 2014
Hi Nick,
It may be possible that the image information is not being set
correctly, causing the resampled locations to be outside the input
image.
The rescale slope and intercept tags, 0028|1052, and 0028|1053,
probably cannot be just copied from the first slice -- they likely
change with slice.
Creating derived DICOM series like this is generally not recommended.
A format like MetaImage is preferred.
Hope this helps,
Matt
On Wed, Apr 16, 2014 at 8:57 AM, Dr Nick Patterson
<pattersonnp.work at gmail.com> wrote:
> I am still having further problems with resampling a DICOM image and have
> adapted the ResampleDICOM example in order to resample a PET image from
> 512x512 to 128x128.
>
> The output of DICOM files are not correct as they simply return a DICOM
> image with 0 everywhere. I seem to be absolutely stuck on this. I have had
> to forcibly set the origin and spacing as this doesn’t seem to be
> automatically extracted from the files.
>
> The same code (changing only the PixelType typedef to unsigned char) works
> perfectly for resampling CT images but not for my PET images. Can anyone
> help me with this as I am getting fairly frustrated with it all!
>
> (Note, there is some use of Qt stuff in here too)
>
> Can anyone see why this does not work for resampling PET DICOM?
>
> Regards, Nick.
>
>
> #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 "itkIdentityTransform.h"
>
> #include "itkLinearInterpolateImageFunction.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>
>
>
> #include "RaydoseManager.h"
>
> #include <iostream>
>
> using namespace std;
>
> static void CopyDictionary (itk::MetaDataDictionary &fromDict,
> itk::MetaDataDictionary &toDict);
>
>
> RaydoseManager *RDManager = RaydoseManager::GetManager();
>
>
> RaydoseResampleDICOM::RaydoseResampleDICOM(RaydoseSeries* parent, int x, int
> y)
>
> {
>
> // Calculate the new X and Y pixel width based on the supplied scale
> factors.
>
> parentSeries = parent;
>
> cout << parent->GetModality().toStdString() << endl;
>
> pixelXWidth = (parentSeries->GetXPixWidth())*x;
>
> pixelYWidth = (parentSeries->GetYPixWidth())*y;
>
> Resample();
>
> }
>
>
> int RaydoseResampleDICOM::Resample() {
>
> const unsigned int InputDimension = 3;
>
> const unsigned int OutputDimension = 2;
>
>
> //typedef unsigned short PixelType;
>
>
> typedef unsigned int 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
>
> InputImageType::Pointer image = InputImageType::New();
>
> ImageIOType::Pointer gdcmIO = ImageIOType::New();
>
> InputNamesGeneratorType::Pointer inputNames =
> InputNamesGeneratorType::New();
>
> inputNames->SetInputDirectory(
> parentSeries->GetDirectoryPath().toLocal8Bit() );
>
>
> 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;
>
> }
>
>
> image = reader->GetOutput();
>
>
> ////////////////////////////////////////////////
>
> // 2) Resample the series
>
> InterpolatorType::Pointer interpolator = InterpolatorType::New();
>
> TransformType::Pointer transform = TransformType::New();
>
> transform->SetIdentity();
>
>
> //const InputImageType::SpacingType& inputSpacing =
> reader->GetOutput()->GetSpacing();
>
> double *inputSpacing = new double[3];
>
> inputSpacing[0] = parentSeries->GetXPixWidth();
>
> inputSpacing[1] = parentSeries->GetYPixWidth();
>
> inputSpacing[2] = parentSeries->GetSliceThickness();
>
>
> const InputImageType::RegionType& inputRegion =
> reader->GetOutput()->GetLargestPossibleRegion();
>
> const InputImageType::SizeType& inputSize = inputRegion.GetSize();
>
>
> InputImageType::SpacingType outputSpacing;
>
> outputSpacing[0] = pixelXWidth;
>
> outputSpacing[1] = pixelYWidth;
>
> outputSpacing[2] = inputSpacing[2];
>
>
> 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);
>
>
> ReaderType::DictionaryRawPointer inputDict =
> (*(reader->GetMetaDataDictionaryArray()))[0];
>
> ReaderType::DictionaryArrayType outputArray;
>
>
> std::string originFirstPixel;
>
> itk::ExposeMetaData<std::string>(*inputDict, "0020|0032",
> originFirstPixel);
>
>
> QString firstPixelPostion(originFirstPixel.c_str());
>
> double *origin = new double[3];
>
> origin[0] =
> firstPixelPostion.mid(0,firstPixelPostion.indexOf("\\")).toDouble();
>
> origin[1] = firstPixelPostion.mid(firstPixelPostion.indexOf("\\")+1,
>
>
> (firstPixelPostion.lastIndexOf("\\")-firstPixelPostion.indexOf("\\"))-1).toDouble();
>
> origin[2] =
> firstPixelPostion.mid(firstPixelPostion.lastIndexOf("\\")+1,firstPixelPostion.length()).toDouble();
>
>
>
> ResampleFilterType::Pointer resampler = ResampleFilterType::New();
>
> resampler->SetInput(reader->GetOutput());
>
> resampler->SetTransform(transform);
>
> resampler->SetInterpolator(interpolator);
>
> resampler->SetOutputOrigin(origin);
>
> 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
>
>
>
> // 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];
>
>
> 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("");
>
> 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];
>
> cout << position[0] << "\\" << position[1] << "\\" << position[2]
> <<endl;
>
>
>
> 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
>
> cout << " The problem occurs here !!!" << endl;
>
> // Make the output directory and generate the file names.
>
> const char* dirID = new char[200];
>
> QString directoryID =
> RDManager->GetInputPath()+parentSeries->GetModality()+parentSeries->GetDate().toString("ddmmyyyy")+"_"+
>
> parentSeries->GetTime().toString("hhMMss");
>
> dirID = directoryID.toLocal8Bit();
>
> itksys::SystemTools::MakeDirectory( dirID );
>
> cout << "Writing out to path: " << directoryID.toStdString() << endl;
>
> // Generate the file names
>
> OutputNamesGeneratorType::Pointer outputNames =
> OutputNamesGeneratorType::New();
>
> QString seriesFormat = directoryID;
>
>
> if (parentSeries->GetModality()=="CT") { seriesFormat = seriesFormat +
> "/" + "CT%d.dcm"; }
>
> else { seriesFormat = seriesFormat + "/" + "NM%d.dcm"; }
>
> outputNames->SetSeriesFormat (seriesFormat.toLocal8Bit());
>
> outputNames->SetStartIndex(1);
>
> outputNames->SetEndIndex(outputSize[2]);
>
>
> SeriesWriterType::Pointer seriesWriter = SeriesWriterType::New();
>
> seriesWriter->SetInput( resampler->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;
>
> }
>
> 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;
>
> }
>
> }
>
>
>
>
>
>
> _____________________________________
> 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://www.itk.org/mailman/listinfo/insight-users
>
More information about the Insight-users
mailing list