[ITK] [ITK-users] Resampling CT is fine... but not PET. Any clues?
Dr Nick Patterson
pattersonnp.work at gmail.com
Wed Apr 16 08:57:03 EDT 2014
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;
}
}
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/community/attachments/20140416/81b68cae/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://www.itk.org/mailman/listinfo/insight-users
More information about the Community
mailing list