[Insight-users] DICOM orientation
Bill Lorensen
bill.lorensen at gmail.com
Tue Oct 6 11:55:19 EDT 2009
John,
I've attached some code thta passes the metat data dictionary from the
input to the output. This example does a resample, bout you should do
a similar thing in your program.
Bill
On Tue, Oct 6, 2009 at 11:38 AM, John Drozd <john.drozd at gmail.com> wrote:
> I have attached a screen capture of the Slicer views. The Slicer views on
> the left are the original DICOM series (that are not flipped and fine). The
> Slicer Views on the right are the outputted volume (after reading the series
> from memory) and are flipped. Note the -ve vs +ve values on the sliders,
> indicating the radiologist vs neurologist views.
>
> John
>
> On Tue, Oct 6, 2009 at 10:21 AM, John Drozd <john.drozd at gmail.com> wrote:
>>
>> I invoke the program with:
>>
>> ./DeformableRegistration "param.file" "m000-talairach.dcm"
>> "/trumpet/downloads/DeformableRegistration_Plugin/DeformableRegistration/datasubject"
>>
>> where the last item on the list is the name of the directory containing
>> the dicom series.
>>
>> john
>>
>> On Tue, Oct 6, 2009 at 10:07 AM, John Drozd <john.drozd at gmail.com> wrote:
>>>
>>> Hi Bill,
>>>
>>> Yes, I did read in the original DICOM series into 3D Slicer and the
>>> original series was not flipped and fine. I was checking if ITK had read
>>> the series correctly into memory, so I wrote it to a volume. I was writing
>>> the volume out to see what was in memory and when I viewed the outputted
>>> volume in 3D Slicer it was flipped. In this case, I was basing my code on
>>> Examples/IO/DicomSeriesReadImageWrite2.cxx
>>> I plan to today read in the series and write it to a series (instead of a
>>> volume) and see if the outputted series will be flipped (hopefully not) when
>>> I view it in 3D Slicer.
>>> By the way, as shown in Examples/IO/DicomImageReadWrite.cxx, I did read
>>> in a dicom volume (a single file containing the talairach atlas) and wrote
>>> it out to another dicom volume and the outputted volume this time was not
>>> flipped when I viewed it in 3D Slicer.
>>>
>>> Below is my code that I used for reading a dicom series and writing to a
>>> volume:
>>> (It is part of a deformable registration program that I am writing based
>>> on Examples/Registration/DeformableRegistration1.cxx so there are some
>>> header files and code from that)
>>> The pertinent reading and writing code is in red. (about 60 lines down)
>>>
>>>
>>> #if defined(_MSC_VER)
>>> #pragma warning ( disable : 4786 )
>>> #endif
>>>
>>>
>>> #include "itkImageFileReader.h"
>>> #include "itkImageFileWriter.h"
>>>
>>> #include "itkRescaleIntensityImageFilter.h"
>>> #include "itkHistogramMatchingImageFilter.h"
>>>
>>> //Added from DicomImageReadWrite.cxx
>>> #include "itkGDCMImageIO.h"
>>>
>>> //Added from DicomSeriesReadImageWrite2.cxx
>>> #include "itkOrientedImage.h"
>>> #include "itkGDCMImageIO.h"
>>> #include "itkGDCMSeriesFileNames.h"
>>> #include "itkImageSeriesReader.h"
>>>
>>> #include "itkFEM.h"
>>> #include "itkFEMRegistrationFilter.h"
>>>
>>> // Next, we use \code{typedef}s to instantiate all necessary classes.
>>> We
>>> // define the image and element types we plan to use to solve a
>>> // two-dimensional registration problem. We define multiple element
>>> // types so that they can be used without recompiling the code.
>>> //
>>> //
>>> // Note that in order to solve a three-dimensional registration
>>> // problem, we would simply define 3D image and element types in lieu
>>> // of those above. The following declarations could be used for a 3D
>>> // problem:
>>>
>>> typedef itk::Image<short, 3> fileImage3DType;
>>> typedef itk::Image<short, 3> Image3DType;
>>>
>>> typedef itk::fem::Element3DC0LinearHexahedronMembrane Element3DType;
>>> typedef itk::fem::Element3DC0LinearTetrahedronMembrane Element3DType2;
>>>
>>> // Here, we instantiate the load types and explicitly template the
>>> // load implementation type. We also define visitors that allow the
>>> // elements and loads to communicate with one another.
>>>
>>> typedef itk::fem::FiniteDifferenceFunctionLoad<Image3DType,Image3DType>
>>> ImageLoadType;
>>> template class itk::fem::ImageMetricLoadImplementation<ImageLoadType>;
>>>
>>> typedef Element3DType::LoadImplementationFunctionPointer LoadImpFP;
>>> typedef Element3DType::LoadType
>>> ElementLoadType;
>>>
>>> typedef Element3DType2::LoadImplementationFunctionPointer LoadImpFP2;
>>> typedef Element3DType2::LoadType
>>> ElementLoadType2;
>>>
>>> typedef itk::fem::VisitorDispatcher<Element3DType,ElementLoadType,
>>> LoadImpFP>
>>>
>>> DispatcherType;
>>>
>>> typedef itk::fem::VisitorDispatcher<Element3DType2,ElementLoadType2,
>>> LoadImpFP2>
>>>
>>> DispatcherType2;
>>>
>>> // Once all the necessary components have been instantiated, we can
>>> // instantiate the \doxygen{FEMRegistrationFilter}, which depends on the
>>> // image input and output types.
>>>
>>> typedef itk::fem::FEMRegistrationFilter<Image3DType,Image3DType>
>>> RegistrationType;
>>>
>>> int main(int argc, char *argv[])
>>> {
>>> char *paramname;
>>> if ( argc < 2 )
>>> {
>>> std::cout << "Parameter file name missing" << std::endl;
>>> std::cout << "Usage: " << argv[0] << " param.file" << " atlas.file"
>>> << " subject.file" <<std::endl;
>>> return EXIT_FAILURE;
>>> }
>>> else
>>> {
>>> paramname=argv[1];
>>> }
>>>
>>> // The \doxygen{fem::ImageMetricLoad} must be registered before it
>>> // can be used correctly with a particular element type. An example
>>> // of this is shown below for ElementType. Similar
>>> // definitions are required for all other defined element types.
>>>
>>>
>>> // Register the correct load implementation with the element-typed
>>> visitor dispatcher.
>>> {
>>> // Software Guide : BeginCodeSnippet
>>> Element3DType::LoadImplementationFunctionPointer fp =
>>>
>>> &itk::fem::ImageMetricLoadImplementation<ImageLoadType>::ImplementImageMetricLoad;
>>> DispatcherType::RegisterVisitor((ImageLoadType*)0,fp);
>>> // Software Guide : EndCodeSnippet
>>> }
>>> {
>>> Element3DType2::LoadImplementationFunctionPointer fp =
>>>
>>> &itk::fem::ImageMetricLoadImplementation<ImageLoadType>::ImplementImageMetricLoad;
>>> DispatcherType2::RegisterVisitor((ImageLoadType*)0,fp);
>>> }
>>>
>>> // In order to begin the registration, we declare an instance of the
>>> // FEMRegistrationFilter. For simplicity, we will call
>>> // it \code{registrationFilter}.
>>>
>>> RegistrationType::Pointer registrationFilter = RegistrationType::New();
>>>
>>> typedef short InputPixelType;
>>> const unsigned int InputDimension = 3;
>>>
>>> typedef itk::Image< InputPixelType, InputDimension > InputImageType;
>>>
>>> typedef itk::ImageSeriesReader< InputImageType > ReaderSeriesType;
>>>
>>>
>>> ReaderSeriesType::Pointer fixedsubjectfilter = ReaderSeriesType::New();
>>>
>>> typedef itk::GDCMImageIO ImageIOType;
>>>
>>> ImageIOType::Pointer gdcmImageIO = ImageIOType::New();
>>>
>>> fixedsubjectfilter->SetImageIO( gdcmImageIO );
>>>
>>> typedef itk::GDCMSeriesFileNames NamesGeneratorType;
>>> NamesGeneratorType::Pointer nameGenerator = NamesGeneratorType::New();
>>>
>>> nameGenerator->SetUseSeriesDetails( true );
>>> nameGenerator->AddSeriesRestriction("0008|0021" );
>>>
>>> nameGenerator->SetDirectory( argv[3] );
>>>
>>> try
>>> {
>>> std::cout << std::endl << "The directory: " << std::endl;
>>> std::cout << std::endl << argv[3] << std::endl << std::endl;
>>> std::cout << "Contains the following DICOM Series: ";
>>> std::cout << std::endl << std::endl;
>>>
>>> typedef std::vector< std::string > SeriesIdContainer;
>>>
>>> const SeriesIdContainer & seriesUID = nameGenerator->GetSeriesUIDs();
>>>
>>> //std::cout << seriesUID.begin() << endl;
>>>
>>> SeriesIdContainer::const_iterator seriesItr = seriesUID.begin();
>>> SeriesIdContainer::const_iterator seriesEnd = seriesUID.end();
>>> while( seriesItr != seriesEnd )
>>> {
>>> std::cout << seriesItr->c_str() << std::endl;
>>> seriesItr++;
>>> }
>>>
>>> std::string seriesIdentifier;
>>>
>>> if( argc > 4 ) // If no optional series identifier
>>> {
>>> seriesIdentifier = argv[3];
>>> }
>>> else
>>> {
>>> seriesIdentifier = seriesUID.begin()->c_str();
>>> }
>>>
>>> std::cout << std::endl << std::endl;
>>> std::cout << "Now reading series: " << std::endl << std::endl;
>>> std::cout << seriesIdentifier << std::endl;
>>> std::cout << std::endl << std::endl;
>>>
>>> typedef std::vector< std::string > FileNamesContainer;
>>> FileNamesContainer fileNames;
>>>
>>> fileNames = nameGenerator->GetFileNames( seriesIdentifier );
>>>
>>> fixedsubjectfilter->SetFileNames( fileNames );
>>>
>>> try
>>> {
>>> fixedsubjectfilter->Update();
>>> std::cout << "Subject read successfully" << std::endl;
>>> }
>>> catch (itk::ExceptionObject &ex)
>>> {
>>> std::cout << ex << std::endl;
>>> return EXIT_FAILURE;
>>> }
>>>
>>> typedef itk::ImageFileWriter< InputImageType > WriterSubjectType;
>>>
>>> WriterSubjectType::Pointer fixedsubjectfilterwriter =
>>> WriterSubjectType::New();
>>>
>>> fixedsubjectfilterwriter->SetFileName( "subjectout.dcm" );
>>>
>>> fixedsubjectfilterwriter->SetInput( fixedsubjectfilter->GetOutput() );
>>>
>>> try
>>> {
>>> fixedsubjectfilterwriter->Update();
>>> }
>>> catch (itk::ExceptionObject &ex)
>>> {
>>> std::cout << ex << std::endl;
>>> return EXIT_FAILURE;
>>> }
>>> }
>>> catch (itk::ExceptionObject &ex)
>>> {
>>> std::cout << ex << std::endl;
>>> return EXIT_FAILURE;
>>> }
>>>
>>>
>>>
>>> return EXIT_SUCCESS;
>>> }
>>>
>>>
>>> john
>>>
>>> On Mon, Oct 5, 2009 at 7:09 PM, Bill Lorensen <bill.lorensen at gmail.com>
>>> wrote:
>>>>
>>>> 3D Slicer can read the DICOM directly. It will use the direction
>>>> information properly. Why are you reading and writing it? I am
>>>> suprised (and concerned) that 3D SLicer would flip the images.
>>>>
>>>> Bill
>>>>
>>>> On Mon, Oct 5, 2009 at 6:13 PM, John Drozd <john.drozd at gmail.com> wrote:
>>>> > Hello,
>>>> >
>>>> > I am reading a brain subject in the form of a DICOM series using ITK
>>>> > and
>>>> > writing it out from memory as a volume.
>>>> > I used the method shown in Examples/IO/DicomSeriesReadImageWrite2.cxx
>>>> > I viewed the original DICOM series in 3D Slicer, and also viewed the
>>>> > outputted volume in 3D Slicer.
>>>> > I am using 3D Slicer because I am developing a segmentation module for
>>>> > 3D
>>>> > Slicer using ITK.
>>>> > I noticed that the original DICOM series and the outputted volume are
>>>> > inverted in orientation, such that the original series is the
>>>> > radiologist
>>>> > view and the outputted volume is the neurosurgean view as per
>>>> > http://www.itk.org/pipermail/insight-users/2009-June/031128.html and
>>>> >
>>>> > http://www.itk.org/Wiki/Proposals:Orientation#DICOM_LPS_Differences_in_Visualization_presented_to_Radiologist_and_NeuroSurgeons
>>>> >
>>>> > Does this have something to do with what I read on the internet that
>>>> > vtk's
>>>> > images are flipped vertically from itk's images?
>>>> > Is there a way to manipulate the image in memory using ITK so that
>>>> > Slicer
>>>> > will display the outputted volume in the same radiologist orientation
>>>> > as the
>>>> > original DICOM series?
>>>> >
>>>> > Thank you.
>>>> >
>>>> > john
>>>> >
>>>> > --
>>>> > John Drozd
>>>> > Postdoctoral Fellow
>>>> > Imaging Research Laboratories
>>>> > Robarts Research Institute
>>>> > Room 1256
>>>> > 100 Perth Drive
>>>> > London, Ontario, Canada
>>>> > N6A 5K8
>>>> > (519) 661-2111 ext. 25306
>>>> >
>>>> > _____________________________________
>>>> > Powered by www.kitware.com
>>>> >
>>>> > Visit other Kitware open-source projects at
>>>> > http://www.kitware.com/opensource/opensource.html
>>>> >
>>>> > 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
>>>> >
>>>> >
>>>
>>>
>>>
>>> --
>>> John Drozd
>>> Postdoctoral Fellow
>>> Imaging Research Laboratories
>>> Robarts Research Institute
>>> Room 1256
>>> 100 Perth Drive
>>> London, Ontario, Canada
>>> N6A 5K8
>>> (519) 661-2111 ext. 25306
>>
>>
>>
>> --
>> John Drozd
>> Postdoctoral Fellow
>> Imaging Research Laboratories
>> Robarts Research Institute
>> Room 1256
>> 100 Perth Drive
>> London, Ontario, Canada
>> N6A 5K8
>> (519) 661-2111 ext. 25306
>
>
>
> --
> John Drozd
> Postdoctoral Fellow
> Imaging Research Laboratories
> Robarts Research Institute
> Room 1256
> 100 Perth Drive
> London, Ontario, Canada
> N6A 5K8
> (519) 661-2111 ext. 25306
>
-------------- next part --------------
/*=========================================================================
Program: Insight Segmentation & Registration Toolkit
Module: $RCSfile: DicomResample.cxx,v $
Language: C++
Date: $Date: 2005/11/19 16:31:50 $
Version: $Revision: 1.9 $
Copyright (c) Insight Software Consortium. All rights reserved.
See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE. See the above copyright notices for more information.
=========================================================================*/
#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif
// Resample a DICOM study
// Usage: DicomResample InputDirectory OutputDirectory
// xSpacing ySpacing zSpacing
//
// Example: DicomResample 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:
// DicomResample 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 "itkOrientedImage.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>
#include "gdcm/src/gdcmFile.h"
#include "gdcm/src/gdcmUtil.h"
#include <string>
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;
}
const unsigned int InputDimension = 3;
const unsigned int OutputDimension = 2;
typedef signed short PixelType;
typedef itk::OrientedImage< PixelType, InputDimension >
InputImageType;
typedef itk::OrientedImage< 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();
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] = atof(argv[3]);
outputSpacing[1] = atof(argv[4]);
outputSpacing[2] = atof(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;
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);
ResampleFilterType::Pointer 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.
std::string seriesUID = gdcm::Util::CreateUniqueUID( gdcmIO->GetUIDPrefix());
std::string frameOfReferenceUID = gdcm::Util::CreateUniqueUID( gdcmIO->GetUIDPrefix());
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);
std::string sopInstanceUID = gdcm::Util::CreateUniqueUID( gdcmIO->GetUIDPrefix());
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 (unsigned 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);
}
////////////////////////////////////////////////
// 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
// 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]);
SeriesWriterType::Pointer seriesWriter = SeriesWriterType::New();
seriesWriter->SetInput( shiftScale->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;
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;
}
}
More information about the Insight-users
mailing list