[Insight-users] Bad rescale slope/intercept with gdcm 2.0.15 via itkGDCMImageIO

[Jesús Spí­nola]

Hi,

So far I was using gdcm 2.0.14 and I was able to read any image with no
problem, but I upgraded recently to gdcm 2.0.15 and I found a problem
reading some CT and MR DICOM image files.
The way I'm reading DICOM files is using itkGDCMImageIO with
itkImageSeriesReader and then filtering the output to vtk for the image
rendering.

The problem is that some images are now (gdcm 2.0.15) loaded with a
different scalar range than before (gdcm 2.0.14), hence they are not
displayed properly with the default window/level setting (almost completely
white in most CT cases for instance). I observed that the images with
trouble were those containing rescale intercept/slope tags with some value.
The images with no rescale intercept/slope present are displayed as
expected, so my guess is there is some problem regarding the intercept/slope
rescale computation.

I've written a little program which reads a DICOM image and renders it with
vtk to reproduce the bug. Then I compiled it against gdcm 2.0.14 and 2.0.15
to see the differences between each other. The source code is at the end of
this e-mail.

I've been testing with some DICOM images from
http://www.creatis.insa-lyon.fr/~jpr/PUBLIC/gdcm/gdcmData.tar.gz. Some of
these test images are valid to reproduce the problem, in this case I'll show
you the results with
MR-MONO2-12-shoulder.dcm, CT-SIEMENS-Icone-With-PaletteColor.dcm
and 05119865-mr-siemens-avanto-syngo.dcm.

In the following table I show some of the relevant information of each image
and the scalar range I get building the program against gdcm 2.0.14 and
2.0.15. The right scalar range is the one obtained with gdcm 2.0.14. The
result is specially bad with the MR-MONO2-12-shoulder.dcm image.

*File** **WW,**WL** **Bits Allocated/**Stored** **Rescale Slope/**Intercept*
* **Scalar Range:** **2.0.14 Build           **2.0.15 Build*
MR-MONO2-12-shoulder.dcm 2000,1000 16/12 3.774114, 0.000061 0..2245
-32065..32556
CT-SIEMENS-Icone-With-PaletteColor.dcm 350,40 16/12 1, -1024 -1024..981
0..2005
05119865-mr-siemens-avanto-syngo.dcm 1538,98 16/12 2, -4096 -1186..526
1455..2311


I tested with itk 3.18 and vtk 5.6,  (the vtk version might not make any
sense in the results because the problem is just in the reading stage).
(Look at some screenshots<https://docs.google.com/uc?id=0B-5DPBfN4U4FMzUzNjI5YjQtNTU5MS00ODlkLThhNDMtNjA3NDRiZjUyMjJh&export=download&hl=en>of
the results)

//
// SOURCE CODE
//

//
// You'll need also itkImageToVTKImageFilter class, available in the
InsightApplications<http://voxel.dl.sourceforge.net/sourceforge/itk/InsightApplications-3.20.0.zip>
 package under "Auxiliary\vtk" folder
//

// Usage: main.exe DICOMFilename [WindowWidth WindowLevel]

#include <vtkImageviewer2.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkRenderWindow.h>

#include <itkImage.h>
#include <itkImageSeriesReader.h>
#include <itkGDCMImageIO.h>
#include "itkImageToVTKImageFilter.h" //Converts an ITK image into a VTK
image and plugs a itk data pipeline to a VTK datapipeline.

int main(int argc, char* argv[])
{
    // Voxel type and image dimension
    typedef signed short int VoxelType;
    typedef VoxelType ItkPixelType;
    static const unsigned int VDimension = 3;

    // Itk Image type
    typedef itk::Image<ItkPixelType, VDimension> ItkImageType;
    typedef ItkImageType::Pointer ItkImageTypePointer;

    // Series reader
    typedef itk::ImageSeriesReader<ItkImageType> SeriesReaderType;
    typedef itk::GDCMImageIO ImageIOType;

    // itk2vtk filter
    typedef itk::ImageToVTKImageFilter<ItkImageType> ItkToVtkFilterType;

    /// DICOM Reader
    SeriesReaderType::Pointer m_seriesReader = SeriesReaderType::New();
    ImageIOType::Pointer m_gdcmIO = ImageIOType::New();

    m_seriesReader->SetImageIO(m_gdcmIO);
    ItkToVtkFilterType::Pointer m_itkToVtkFilter =
ItkToVtkFilterType::New();

    // Check the number of parameters given
    if (argc >= 2)
    {
        // Create the renderer object and configure it
        vtkImageViewer2 *viewer = vtkImageViewer2::New();
        vtkRenderWindowInteractor *interactor =
vtkRenderWindowInteractor::New();
        viewer->SetupInteractor(interactor);

        // Set the first argument as the input filename to read
        // Atention: We don't check whether the file is valid or not
        m_seriesReader->SetFileName(argv[1]);
        m_itkToVtkFilter->SetInput(m_seriesReader->GetOutput());
        try
        {
            m_itkToVtkFilter->Update();
        }
        catch (itk::ExceptionObject & exception)
        {
            std::cerr << "Exception caught!" << exception.GetDescription()
<< std::endl;
        }

        // Show some information about the image
        m_itkToVtkFilter->GetOutput()->Print(std::cout);
        std::cerr << "** Scalar Range **" << std::endl;
        std::cerr << m_itkToVtkFilter->GetOutput()->GetScalarRange()[0] <<
".." << m_itkToVtkFilter->GetOutput()->GetScalarRange()[1] << std::endl;

        // Set the viewer input (DICOM File) and do the rendering stuff
        viewer->SetInput(m_itkToVtkFilter->GetOutput());

        // We can optionally set the window/level through the arguments
        if (argc >=4)
        {
            std::cout << "Applying given ww/wl: " << argv[2] << ", " <<
argv[3] << std::endl;
            viewer->SetColorWindow(atof(argv[2]));
            viewer->SetColorLevel(atof(argv[3]));
        }
        else
        {
            std::cout << "Applying default ww/wl: 256, 128" << std::endl;
            viewer->SetColorWindow(256);
            viewer->SetColorLevel(128);
        }

        viewer->Render();
        interactor->Start();
    }
    else
    {
        std::cerr << "Usage: main.exe DICOMFilename [WindowWidth
WindowLevel]" << std::endl;
    }

    return 0;
}

-- 
Jesús Spínola
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20100723/e6e8003b/attachment.htm>