Resample DICOM Series#
Synopsis#
Resample a DICOM series.
Results#
Note
Help Wanted Implementation of Results for sphinx examples containing this message. Reconfiguration of CMakeList.txt may be necessary. Write An Example <https://itk.org/ITKExamples/Documentation/Contribute/WriteANewExample.html>
Code#
C++#
// 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 (only for ITK < 4.6)
// 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"
#if ((ITK_VERSION_MAJOR == 4) && (ITK_VERSION_MINOR < 6))
# include "itkShiftScaleImageFilter.h"
#endif
#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 <sstream>
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;
}
constexpr unsigned int InputDimension = 3;
constexpr unsigned int OutputDimension = 2;
using PixelType = signed short;
using InputImageType = itk::Image<PixelType, InputDimension>;
using OutputImageType = itk::Image<PixelType, OutputDimension>;
using ReaderType = itk::ImageSeriesReader<InputImageType>;
using ImageIOType = itk::GDCMImageIO;
using InputNamesGeneratorType = itk::GDCMSeriesFileNames;
using OutputNamesGeneratorType = itk::NumericSeriesFileNames;
using TransformType = itk::IdentityTransform<double, InputDimension>;
using InterpolatorType = itk::LinearInterpolateImageFunction<InputImageType, double>;
using ResampleFilterType = itk::ResampleImageFilter<InputImageType, InputImageType>;
#if ((ITK_VERSION_MAJOR == 4) && (ITK_VERSION_MINOR < 6))
using ShiftScaleType = itk::ShiftScaleImageFilter<InputImageType, InputImageType>;
#endif
using SeriesWriterType = itk::ImageSeriesWriter<InputImageType, OutputImageType>;
////////////////////////////////////////////////
// 1) Read the input series
auto gdcmIO = ImageIOType::New();
auto inputNames = InputNamesGeneratorType::New();
inputNames->SetInputDirectory(argv[1]);
const ReaderType::FileNamesContainer & filenames = inputNames->GetInputFileNames();
auto reader = ReaderType::New();
reader->SetImageIO(gdcmIO);
reader->SetFileNames(filenames);
try
{
reader->Update();
}
catch (const 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
auto interpolator = InterpolatorType::New();
auto 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] = std::stod(argv[3]);
outputSpacing[1] = std::stod(argv[4]);
outputSpacing[2] = std::stod(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;
using SizeValueType = InputImageType::SizeType::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);
auto 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.
#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
auto 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
std::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);
}
#if ((ITK_VERSION_MAJOR == 4) && (ITK_VERSION_MINOR < 6))
////////////////////////////////////////////////
// 4) Shift data to undo the effect of a rescale intercept by the
// DICOM reader
std::string interceptTag("0028|1052");
using MetaDataStringType = itk::MetaDataObject<std::string>;
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());
}
auto shiftScale = ShiftScaleType::New();
shiftScale->SetInput(resampler->GetOutput());
shiftScale->SetShift(interceptShift);
#endif
////////////////////////////////////////////////
// 5) Write the new DICOM series
// Make the output directory and generate the file names.
itksys::SystemTools::MakeDirectory(argv[2]);
// Generate the file names
auto 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]);
auto seriesWriter = SeriesWriterType::New();
#if ((ITK_VERSION_MAJOR == 4) && (ITK_VERSION_MINOR < 6))
seriesWriter->SetInput(shiftScale->GetOutput());
#else
seriesWriter->SetInput(resampler->GetOutput());
#endif
seriesWriter->SetImageIO(gdcmIO);
seriesWriter->SetFileNames(outputNames->GetFileNames());
seriesWriter->SetMetaDataDictionaryArray(&outputArray);
try
{
seriesWriter->Update();
}
catch (const 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)
{
using DictionaryType = itk::MetaDataDictionary;
DictionaryType::ConstIterator itr = fromDict.Begin();
DictionaryType::ConstIterator end = fromDict.End();
using MetaDataStringType = itk::MetaDataObject<std::string>;
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;
}
}
Classes demonstrated#
-
class GDCMImageIO : public itk::ImageIOBase
ImageIO class for reading and writing DICOM V3.0 and ACR/NEMA 1&2 uncompressed images. This class is only an adaptor to the GDCM library.
GDCM can be found at: https://sourceforge.net/projects/gdcm
To learn more about the revision shipped with ITK, call
git log Modules/ThirdParty/GDCM/src/
from an ITK Git checkout.
The compressors supported include “JPEG2000” (default), and “JPEG”. The compression level parameter is not supported.
DICOM tags are represented as strings in the metadata dictionary. The string format is “XXXX,XXXX”, DICOM group number followed by element number, both are hexadecimals. The separator character is either a pipe “|” or a comma “,”.
- Warning
As the metadata dictionary uses the DICOM tag strings as keys it is possible to have multiple entries representing the same DICOM tag with different values. The last one encountered will be used when writing the image to file. For example, the patient name tag “0010,0010” and “0010|0010” may both be in the dictionary with different values due to the different separator character. Similarly, the series description tag may appear multiple times as “0008,103e”, “0008,103E”, “0008|103e”, or “0008|103E”. The strings differ in the separator character and lower or upper case letters in the hexadecimal numbers. Note that when read from file, letters in the hexadecimal numbers are always set to lower case and the separator character is a pipe. To ensure consistency, it is best to always use lower case letters and the pipe separator when explicitly adding tags to the metadata dictionary.
- Warning
There are several restrictions to this current writer:
Even though during the writing process you pass in a DICOM file as input The output file may not contains ALL DICOM field from the input file. In particular:
The SeQuence DICOM field (SQ).
Fields from Private Dictionary.
Some very long (>0xfff) binary fields are not loaded (typically 0029|0010), you need to explicitly set the maximum length of elements to load to be bigger (see Get/SetMaxSizeLoadEntry).
In DICOM some fields are stored directly using their binary representation. When loaded into the MetaDataDictionary some fields are converted to ASCII (only VR: OB/OW/OF and UN are encoded as mime64).
- ITK Sphinx Examples:
