[Insight-users] Write Rescale Intercept and Rescale slope to Dicom file
LiviaBarazzetti
barazzetti.livia at yahoo.it
Fri May 3 12:38:23 EDT 2013
Hi,
I am trying to read a Dicom series, apply the RegionOfInterestImageFilter to
the image, copy the dicom metadata from the original series to the new one,
and save the result again as a Dicom series.
I would like that the resulting images have the same pixel values and the
same Rescale Intercept (0028,1052) and Rescale slope (0028,1052) as the
original ones.
If I declare the pixeltype of my image as "signed short", the resulting
Dicom files always have Rescale slope =1 and Rescale intercept = 0.
It seems that the rescale is considered only when the pixel type is FLOAT32
or FLOAT34 (see lines 890 and 927 of itkGDCMImageIO.cxx)
So, I declared my pixeltype as double. It works well with a input Dicom file
which has Intercept = -1000 and slope=1: the pixel values and the tags
Intercept and Slope are maintained (but the smallest image pixel value
0028,0106, the largest image pixel value 0028,0107 and the pixel data
7FE0,0010 are different)
With another Dicom file which has Rescale Intercept= -1.000000e+03 and
Rescale Slope = 6.221728e-01 the program gives Access violation (gdcm
Rescaler.cxx line 107). For information, in this image the pixel values are
from -3915 to 13544.
Here below is my code, which was originally taken from
http://www.itk.org/Wiki/ITK/Examples/DICOM/ResampleDICOM (but I had to
remove the shiftscaleimagefilter)
Thanks,
Best regards
Livia
#include "itkVersion.h"
#include "itkImage.h"
#include "itkGDCMImageIO.h"
#include "itkGDCMSeriesFileNames.h"
#include "itkImageSeriesReader.h"
#include "itkImageSeriesWriter.h"
#include <itksys/SystemTools.hxx>
#include "itkShiftScaleImageFilter.h"
#if ITK_VERSION_MAJOR >= 4
#include "gdcmUIDGenerator.h"
#else
#include "gdcm/src/gdcmFile.h"
#include "gdcm/src/gdcmUtil.h"
#endif
#include "itkStatisticsImageFilter.h"
#include <string>
static void CopyDictionary (itk::MetaDataDictionary &fromDict,
itk::MetaDataDictionary &toDict);
int main( int argc, char* argv[] )
{
const unsigned int InputDimension = 3;
const unsigned int OutputDimension = 2;
typedef double PixelType;
//typedef signed short 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
NamesGeneratorType;
typedef itk::ImageSeriesWriter< InputImageType, OutputImageType >
SeriesWriterType;
typedef InputImageType::SizeType::SizeValueType
SizeValueType;
typedef itk::ShiftScaleImageFilter< InputImageType, InputImageType >
ShiftScaleType;
// Validate input parameters
if( argc < 3 ) //9
{
std::cerr << "Usage: "
<< argv[0]
<< " InputDicomDirectory "
<< " OutputDicomDirectory "
<< std::endl;
return EXIT_FAILURE;
}
std::string myInputDirectory;
std::string myOutputDirectory;
int myCropLimits[6];
try
{
myInputDirectory = argv[1];
myOutputDirectory = argv[2];
}
catch (itk::ExceptionObject &excp)
{
std::cerr << "Exception thrown while reading the input parameters" <<
std::endl;
std::cerr << excp << std::endl;
return EXIT_FAILURE;
}
////////////////////////////////////////////////
// Read the input series
InputImageType::Pointer image = InputImageType::New();
ImageIOType::Pointer gdcmIO = ImageIOType::New();
NamesGeneratorType::Pointer namesGenerator = NamesGeneratorType::New();
namesGenerator->SetInputDirectory( myInputDirectory.c_str() );
ReaderType::Pointer reader = ReaderType::New();
reader->SetImageIO( gdcmIO );
reader->SetFileNames( namesGenerator->GetInputFileNames());
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;
}
const InputImageType::RegionType& inputRegion =
reader->GetOutput()->GetLargestPossibleRegion();
const InputImageType::SizeType& inputSize =
inputRegion.GetSize();
image = reader->GetOutput();
InputImageType::SizeType outputSize;
std::cout<< "Image has been read" << std::endl;
///////////////////////////////////////////////
// Apply the ROIFilter
//....
///////////////////////////////////////////////
// Create a MetaDataDictionary for each slice.
std::cout << "Create a MetaDataDictionary for each slice" << std::endl;
// 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 < inputSize[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());
// Derivation Description - How this image was derived
value.str("");
for (int i = 0; i < argc; i++)
{
value << argv[i] << " ";
}
value << ": " << ITK_SOURCE_VERSION;
unsigned 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;
image->TransformIndexToPhysicalPoint(index, position);
value.str("");
value << position[0] << "\\" << position[1] << "\\" << position[2];
itk::EncapsulateMetaData<std::string>(*dict,"0020|0032", value.str());
value.str("");
value << position[2];
itk::EncapsulateMetaData<std::string>(*dict,"0020|1041", value.str());
// Save the dictionary
outputArray.push_back(dict);
}
////////////////////////////////////////////////
// Shift data to undo the effect of a rescale intercept by the
// DICOM reader
/*
std::string intercept;
itk::ExposeMetaData<std::string>(*inputDict, "0020|000d", intercept);
double interceptShift = -atof( intercept.c_str() );
std::string slope;
itk::ExposeMetaData<std::string>(*inputDict, "0020|000d", slope);
double slopeShift = 1/atof( slope.c_str() );
ShiftScaleType::Pointer shiftIntercept= ShiftScaleType::New();
shiftIntercept->SetInput(image);
shiftIntercept->SetShift( interceptShift );
ShiftScaleType::Pointer shiftSlope= ShiftScaleType::New();
shiftSlope->SetInput(shiftIntercept->GetOutput());
shiftSlope->SetScale( slopeShift);
shiftSlope->Update();
image=shiftSlope->GetOutput();
*/
//////////////////////////////////////////////
// Write the new DICOM series
std::cout << "Writing the new DICOM series " << std::endl;
itksys::SystemTools::MakeDirectory( myOutputDirectory.c_str() );
namesGenerator->SetOutputDirectory( myOutputDirectory.c_str() );
SeriesWriterType::Pointer seriesWriter = SeriesWriterType::New();
seriesWriter->SetInput(image);
seriesWriter->SetImageIO( gdcmIO );
seriesWriter->SetFileNames(namesGenerator->GetOutputFileNames()) ;
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;
}
}
--
View this message in context: http://itk-insight-users.2283740.n2.nabble.com/Write-Rescale-Intercept-and-Rescale-slope-to-Dicom-file-tp7582942.html
Sent from the ITK Insight Users mailing list archive at Nabble.com.
More information about the Insight-users
mailing list