ITK  6.0.0
Insight Toolkit
SphinxExamples/src/IO/GDCM/ResampleDICOMSeries/Code.cxx
/*=========================================================================
*
* Copyright NumFOCUS
*
* Licensed under the Apache License, Version 2.0 (the "License");
* you may not use this file except in compliance with the License.
* You may obtain a copy of the License at
*
* https://www.apache.org/licenses/LICENSE-2.0.txt
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*
*=========================================================================*/
// 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 "itkGDCMImageIO.h"
#if ((ITK_VERSION_MAJOR == 4) && (ITK_VERSION_MINOR < 6))
#endif
#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 ImageIOType = itk::GDCMImageIO;
using InputNamesGeneratorType = itk::GDCMSeriesFileNames;
using OutputNamesGeneratorType = itk::NumericSeriesFileNames;
#if ((ITK_VERSION_MAJOR == 4) && (ITK_VERSION_MINOR < 6))
#endif
// 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;
}
}
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.
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;
}
}
Pointer
SmartPointer< Self > Pointer
Definition: itkAddImageFilter.h:93
itkImageSeriesWriter.h
itk::MetaDataObject
Allows arbitrary data types to be stored as MetaDataObjectBase types, and to be stored in a MetaDataD...
Definition: itkMetaDataObject.h:68
ConstPointer
SmartPointer< const Self > ConstPointer
Definition: itkAddImageFilter.h:94
itk::IdentityTransform
Implementation of an Identity Transform.
Definition: itkIdentityTransform.h:50
itkLinearInterpolateImageFunction.h
itk::GTest::TypedefsAndConstructors::Dimension2::PointType
ImageBaseType::PointType PointType
Definition: itkGTestTypedefsAndConstructors.h:51
itk::ImageSeriesWriter
Writes image data to a series of data files.
Definition: itkImageSeriesWriter.h:85
itk::ImageSeriesReader
Data source that reads image data from a series of disk files.
Definition: itkImageSeriesReader.h:45
itk::GTest::TypedefsAndConstructors::Dimension2::SizeType
ImageBaseType::SizeType SizeType
Definition: itkGTestTypedefsAndConstructors.h:49
itkImage.h
itk::SmartPointer< Self >
itkGDCMImageIO.h
itkShiftScaleImageFilter.h
itk::GTest::TypedefsAndConstructors::Dimension2::IndexType
ImageBaseType::IndexType IndexType
Definition: itkGTestTypedefsAndConstructors.h:50
itkImageSeriesReader.h
ITK_SOURCE_VERSION
#define ITK_SOURCE_VERSION
Definition: itkVersion.h:39
itk::LinearInterpolateImageFunction
Linearly interpolate an image at specified positions.
Definition: itkLinearInterpolateImageFunction.h:51
itk::MetaDataDictionary
Provides a mechanism for storing a collection of arbitrary data types.
Definition: itkMetaDataDictionary.h:54
itkNumericSeriesFileNames.h
itk::GTest::TypedefsAndConstructors::Dimension2::RegionType
ImageBaseType::RegionType RegionType
Definition: itkGTestTypedefsAndConstructors.h:54
itk::MetaDataDictionary::Begin
Iterator Begin()
itk::SmartPointer::GetPointer
ObjectType * GetPointer() const noexcept
Definition: itkSmartPointer.h:132
itkIdentityTransform.h
itk::GDCMImageIO
ImageIO class for reading and writing DICOM V3.0 and ACR/NEMA 1&2 uncompressed images....
Definition: itkGDCMImageIO.h:103
itkMinimumMaximumImageFilter.h
itk::GDCMSeriesFileNames
Generate a sequence of filenames from a DICOM series.
Definition: itkGDCMSeriesFileNames.h:63
itk::MetaDataDictionary::End
Iterator End()
itk::ShiftScaleImageFilter
Shift and scale the pixels in an image.
Definition: itkShiftScaleImageFilter.h:40
itkVersion.h
itk::NumericSeriesFileNames
Generate an ordered sequence of filenames.
Definition: itkNumericSeriesFileNames.h:54
itkGDCMSeriesFileNames.h
itk::ResampleImageFilter
Resample an image via a coordinate transform.
Definition: itkResampleImageFilter.h:90
itk::Image
Templated n-dimensional image class.
Definition: itkImage.h:88
New
static Pointer New()
itkResampleImageFilter.h
itk::SizeValueType
unsigned long SizeValueType
Definition: itkIntTypes.h:86
itk::Size::GetSize
const SizeValueType * GetSize() const
Definition: itkSize.h:169