I have experienced that Dicom pixel data are interpreted differently when read using GDCM as IO provider compared to using DCMTK as provider.
The problem arises when the rescale slope and intercept in the dicom files results in non-integer values. GDCM as provider is fully able to capture this, and provides an float image with values as expected. The same image, when interpreted with DCMTK, yields a float image with floored values, i.e. 10.7 becomes 10.0.
A bit of sample code is provided below, but my question is: Is this expected behavior? Or might there be an unfortunate casting somewhere in the code?
Thanks!
#include "itkImage.h"
#include "itkGDCMImageIO.h"
#include "itkDCMTKImageIO.h"
#include "itkNumericSeriesFileNames.h"
#include "itkImageSeriesReader.h"
#include "itkImageSeriesWriter.h"
#include <vector>
#include <string>
#include <iostream>
#include <iomanip>
#include <sstream>
#include "itksys/SystemTools.hxx"
int main( int argc, char *argv[] )
{
if( argc < 2 )
{
std::cerr << "Missing Parameters " << std::endl;
std::cerr << "Usage: " << argv[0];
std::cerr << " input_ImageFile" << std::endl;
return EXIT_FAILURE;
}
typedef itk::Image< double , 2 > InputImageType;
typedef itk::ImageSeriesReader< InputImageType > ReaderType;
ReaderType::Pointer reader = ReaderType::New();
typedef itk::GDCMImageIO GDCMImageIOType;
typedef itk::DCMTKImageIO DCMTKImageIOType;
GDCMImageIOType::Pointer gdcmDicomIO = GDCMImageIOType::New();
DCMTKImageIOType::Pointer dcmtkDicomIO = DCMTKImageIOType::New();
gdcmDicomIO->SetLoadPrivateTags(true);
std::vector<std::string> fileNames;
fileNames.push_back(std::string(argv[1]));
// DCMTK first
reader->SetImageIO( dcmtkDicomIO );
reader->SetFileNames( fileNames );
try {
reader->Update();
}
catch(...) {
return EXIT_FAILURE;
}
InputImageType::Pointer dcmtk_img=reader->GetOutput();
itk::ImageRegionIterator<InputImageType> data_it( dcmtk_img, dcmtk_img->GetLargestPossibleRegion() );
data_it.GoToBegin();
size_t n_pos=0;
double min_val=12981281.;
double max_val=-12981281.;
while(! data_it.IsAtEnd()) {
if(data_it.Get() > 0.0) {
++n_pos;
}
if(data_it.Get() > max_val) {
max_val=data_it.Get();
}
if(data_it.Get() < min_val) {
min_val=data_it.Get();
}
++data_it;
}
std::cout << "DCMTK: There are " << n_pos << " positive voxels..." << std::endl;
std::cout << "DCMTK: Min value: " << min_val << ", max value: " << max_val << std::endl;
// Then GDCM
reader->SetImageIO( gdcmDicomIO );
reader->SetFileNames( fileNames );
try {
reader->Update();
}
catch(...) {
return EXIT_FAILURE;
}
InputImageType::Pointer gdcm_img=reader->GetOutput();
itk::ImageRegionIterator<InputImageType> data_it2( gdcm_img, gdcm_img->GetLargestPossibleRegion() );
data_it2.GoToBegin();
n_pos=0;
min_val=12981281.;
max_val=-12981281.;
while(! data_it2.IsAtEnd()) {
if(data_it2.Get() > 0.0) {
++n_pos;
}
if(data_it2.Get() > max_val) {
max_val=data_it2.Get();
}
if(data_it2.Get() < min_val) {
min_val=data_it2.Get();
}
++data_it2;
}
std::cout << "GDCM: There are " << n_pos << " positive voxels..." << std::endl;
std::cout << "GDCM: Min value: " << min_val << ", max value: " << max_val << std::endl;
return 1;
}