[Insight-users] DICOM orientation

[John Drozd]

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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20091006/4691f33e/attachment-0001.htm>