[Insight-users] Reading Deformation Field File into ITK ImageFileReader Filter

Emma Ryan eryanvtk at yahoo.com
Thu Feb 21 10:00:26 EST 2008


Hi,

  I forgot to put an appropriate heading, so wasn't sure if people were gonne read my email.
Hence the repeat.

Thanks,
Emma


Hi,
 
  I am trying to read in a 3D deformationField.vtk file and apply it to
 a 3D image. To read the deformation field, I use ImageFileReader
 <deformafieldType>  and then connect the output of this to the warper.

The code for this is available below.  It seems to take forever to read
 the deformation field file.  Even after 2 hours, the reading is not
 complete. What is going on ?  Any clues ?

Thanks,
Emma


#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif

#include "itkImageFileReader.h" 
#include "itkImageFileWriter.h" 

#include "itkImageRegionIterator.h"


#include "itkCastImageFilter.h"
#include "itkWarpImageFilter.h"
#include "itkLinearInterpolateImageFunction.h"

#include "itkSquaredDifferenceImageFilter.h"
#include "itkCheckerBoardImageFilter.h"
#include "itkTimeProbe.h"       

#include "itkImageRegionIteratorWithIndex.h"
#include "itkPointSet.h"
#include <fstream>



int main( int argc, char *argv[] )
{
    if( argc < 4 )
    {
        std::cerr << "Missing Parameters " << std::endl;
        std::cerr << "Usage: " << argv[0];
        std::cerr << " fixedImagefile movingImageFile deformField.vtk
 ";
        std::cerr << " outputImageFile " << std::endl;
        return 1;
    }

    
    const unsigned int Dimension = 3;
    typedef unsigned char PixelType;


    typedef itk::Image< PixelType, Dimension >  FixedImageType;
    typedef itk::ImageFileReader< FixedImageType >
 FixedImageReaderType;
    FixedImageReaderType::Pointer fixedImageReader =
 FixedImageReaderType::New();
    fixedImageReader->SetFileName( argv[1] );

    
    typedef itk::Image< PixelType, Dimension >  MovingImageType;
    typedef itk::ImageFileReader< MovingImageType >
 MovingImageReaderType;
    MovingImageReaderType::Pointer movingImageReader =
 MovingImageReaderType::New();
    movingImageReader->SetFileName( argv[2] );

        typedef itk::Vector< float, Dimension >    VectorPixelType;
    typedef itk::Image<  VectorPixelType, Dimension >
 DeformationFieldType;

    typedef itk::ImageFileReader< DeformationFieldType >
  FieldReaderType;
    FieldReaderType::Pointer fieldReader = FieldReaderType::New();

    fieldReader->SetFileName( argv[3] );

   
   try {

    fieldReader->Update();
   }

   catch ( itk::ExceptionObject & err ) 
   {
    std::cerr << "ExceptionObject caught !" << std::endl; 
    std::cerr << err << std::endl; 
    return -1;
    } 


    
    typedef itk::WarpImageFilter<MovingImageType,
 MovingImageType,DeformationFieldType  >     WarperType;
    typedef itk::LinearInterpolateImageFunction<MovingImageType, double
  >  InterpolatorType;
    WarperType::Pointer warper = WarperType::New();
    InterpolatorType::Pointer interpolator = InterpolatorType::New();
    FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();

   
    warper->SetInput( movingImageReader->GetOutput() );
    warper->SetInterpolator( interpolator );
    warper->SetOutputSpacing( fixedImage->GetSpacing() );
    warper->SetOutputOrigin( fixedImage->GetOrigin() );

    warper->SetDeformationField( fieldReader->GetOutput() ); 

  
    typedef  unsigned char  OutputPixelType;
    typedef itk::Image< OutputPixelType, Dimension > OutputImageType;
    typedef itk::CastImageFilter< 
        MovingImageType,
        OutputImageType > CastFilterType;
    typedef itk::ImageFileWriter< OutputImageType >  WriterType;

  
    WriterType::Pointer      writer =  WriterType::New();
    CastFilterType::Pointer  caster =  CastFilterType::New();

    writer->SetFileName( argv[4] );

    
    caster->SetInput( warper->GetOutput() );
    writer->SetInput( caster->GetOutput()   );
    writer->Update();

   
    return 0;
}






    
  ____________________________________________________________________________________
Be a better friend, newshound, and 
know-it-all with Yahoo! Mobile.  Try it now.
  http://mobile.yahoo.com/;_ylt=Ahu06i62sR8HDtDypao8Wcj9tAcJ 
-------------- next part --------------
An HTML attachment was scrubbed...
URL:
 http://public.kitware.com/pipermail/insight-users/attachments/20080221/ae444f83/attachment-0001.html

------------------------------

Message: 3
Date: Thu, 21 Feb 2008 00:41:31 -0800 (PST)
From: Emma Ryan <eryanvtk at yahoo.com>
Subject: [Insight-users] Is it silly to attempt this ?
To: insight-users at itk.org
Cc: insight-users at itk.org
Message-ID: <625726.500.qm at web57908.mail.re3.yahoo.com>
Content-Type: text/plain; charset="us-ascii"

Hi,
 
  I am trying to read in a 3D deformationField.vtk file and apply it to
 a 3D image. To read the deformation field, I use ImageFileReader
 <deformafieldType>  and then connect the output of this to the warper.

The code for this is available below.  It seems to take forever to read
 the deformation field file.  Even after 2 hours, the reading is not
 complete. What is going on ?  Any clues ?

Thanks,
Emma


#if defined(_MSC_VER)
#pragma warning ( disable : 4786 )
#endif

#include "itkImageFileReader.h" 
#include "itkImageFileWriter.h" 

#include "itkImageRegionIterator.h"


#include "itkCastImageFilter.h"
#include "itkWarpImageFilter.h"
#include "itkLinearInterpolateImageFunction.h"

#include "itkSquaredDifferenceImageFilter.h"
#include "itkCheckerBoardImageFilter.h"
#include "itkTimeProbe.h"       

#include "itkImageRegionIteratorWithIndex.h"
#include "itkPointSet.h"
#include <fstream>



int main( int argc, char *argv[] )
{
    if( argc < 4 )
    {
        std::cerr << "Missing Parameters " << std::endl;
        std::cerr << "Usage: " << argv[0];
        std::cerr << " fixedImagefile movingImageFile deformField.vtk
 ";
        std::cerr << " outputImageFile " << std::endl;
        return 1;
    }

    
    const unsigned int Dimension = 3;
    typedef unsigned char PixelType;


    typedef itk::Image< PixelType, Dimension >  FixedImageType;
    typedef itk::ImageFileReader< FixedImageType >
 FixedImageReaderType;
    FixedImageReaderType::Pointer fixedImageReader =
 FixedImageReaderType::New();
    fixedImageReader->SetFileName( argv[1] );

    
    typedef itk::Image< PixelType, Dimension >  MovingImageType;
    typedef itk::ImageFileReader< MovingImageType >
 MovingImageReaderType;
    MovingImageReaderType::Pointer movingImageReader =
 MovingImageReaderType::New();
    movingImageReader->SetFileName( argv[2] );

        typedef itk::Vector< float, Dimension >    VectorPixelType;
    typedef itk::Image<  VectorPixelType, Dimension >
 DeformationFieldType;

    typedef itk::ImageFileReader< DeformationFieldType >
  FieldReaderType;
    FieldReaderType::Pointer fieldReader = FieldReaderType::New();

    fieldReader->SetFileName( argv[3] );

   
   try {

    fieldReader->Update();
   }

   catch ( itk::ExceptionObject & err ) 
   {
    std::cerr << "ExceptionObject caught !" << std::endl; 
    std::cerr << err << std::endl; 
    return -1;
    } 


    
    typedef itk::WarpImageFilter<MovingImageType,
 MovingImageType,DeformationFieldType  >     WarperType;
    typedef itk::LinearInterpolateImageFunction<MovingImageType, double
  >  InterpolatorType;
    WarperType::Pointer warper = WarperType::New();
    InterpolatorType::Pointer interpolator = InterpolatorType::New();
    FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();

   
    warper->SetInput( movingImageReader->GetOutput() );
    warper->SetInterpolator( interpolator );
    warper->SetOutputSpacing( fixedImage->GetSpacing() );
    warper->SetOutputOrigin( fixedImage->GetOrigin() );

    warper->SetDeformationField( fieldReader->GetOutput() ); 

  
    typedef  unsigned char  OutputPixelType;
    typedef itk::Image< OutputPixelType, Dimension > OutputImageType;
    typedef itk::CastImageFilter< 
        MovingImageType,
        OutputImageType > CastFilterType;
    typedef itk::ImageFileWriter< OutputImageType >  WriterType;

  
    WriterType::Pointer      writer =  WriterType::New();
    CastFilterType::Pointer  caster =  CastFilterType::New();

    writer->SetFileName( argv[4] );

    
    caster->SetInput( warper->GetOutput() );
    writer->SetInput( caster->GetOutput()   );
    writer->Update();

   
    return 0;
}


Luis Ibanez <luis.ibanez at kitware.com> wrote: 
Hi Emma,

    The origin coordinate do not have to be 0,0,0.

In fact, the reason for storing the origin in the DICOM
header is that the values are different from (0,0,0).

In general the values are reported relative to some origin
of coordinates in the scanner itself. Therefore they can
be positive or negative and they can easily reach several
centimeters in magnitude.

Note that the origin is measured in physical units.
Typically millimeters in a CT scanner. Not in pixels.


     Regards


        Luis


------------------------
Emma Ryan wrote:
> Hi Luis,
> 
>   Thank you again for your reply. With regard to " Yes, the Patient 
> position is the source of the data for the
> image origin values" in your reply, how does -255, -255, -30 make sense 
> ? I'm not sure I understand the logic behind the "Image Patient 
> Position" tag in a dicom file (0020 - 0032 bytes) .
> The IPP in my dicom file is listed as -255, -255, -33. Unless the images 
> are centered, should these values be (0, 0, 0) for the first slice ?
> 
> Emma
> 
> 
> Could you point me to some references if there is a lengthy explanation 
> for this ?
> 
> 
> 
> Question 1:
> 
> When the Similarity3DTransform gets a negative value in the Scale
> parameter, that performs a reflection of the space on the origin
> of coordinates.... not a good thing at all, particularly in medical
> images.
> 
> The scale shouldn't approach a zero value either...
> 
> This however, is possible when you use an optimizer that assumes
> that the parameters of the transform are in a Vector space, as
> most of the GradientDescent optimizers do.
> 
> 
> In practice you could have prevented the scale from deforming that
> far by setting proper values in the parameter scaling array that
> you pass to the optimizer. The purpose of this array is to let
> the optimizer know that we don't expect big changes to occur in
> the Scale parameter of Transform.
> 
> 
> 
> Question 2:
> 
> You should only need a Similarity transform is you are registering
> images from different patients, that will naturally have slightly
> different organ sizes. Note however that the size variability of
> human specimens is not as large as to need scales smaller than 0.5
> or larger than 2.0.
> 
> If you are registering two different patients, then, in addition
> to having the origin & spacing taken into account, you may also
> need this type of scaling in space.
> 
> 
> Question 3:
> 
> Resampling will take the image spacing into account. You just
> need to make sure that the spacing is set correctly in the images.
> Note that the resampling also use the transform anyways.
> 
> 
> Question 4:
> 
> No, you don't need to manually compensate for the resolution.
> ITK already takes the pixel spacing into account. All the
> process of image registration is done in physical coordinates.
> 
> 
> Question 5:
> 
> Yes, the Patient position is the source of the data for the
> image origin values.
> 
> 
> 
> 
>     Regards,
> 
> 
>        Luis
> 
> 
> 
> -----------------
> Emma Ryan wrote:
>  > Hi,
>  >
>  >   I am using the itk Similarity3DTranform to obtain a registration
>  > between two volumes.
>  > Dataset1 : 512 x 512 x 30 (Resolution = 0.684 x 0.684 x 1.0) :  Image
>  > Patient Position : -175, -175, 125 (from dicom header)
>  > Dataset2: 512 x 512 x 30 (Resolution = 1.0 x 1.0 x 1.0 ) : Image Patient
>  > Position: -255, -255, 101 (from dicom header)
>  >
>  > Both are of the same modality (CT).
>  >
>  > I have a few questions.
>  >
>  > 1. How does one interpret the scale value of  (-0.00008). What does
>  > negative mean ?  Division ?
>  >
>  > 2. Prior to registration I mention the size, origin and spacing for both
>  > the fixed and moving volumes. At resampling, the moving volume are
>  > resampled to match the spacing and origin of the fixed volume. So should
>  > I care to include the scale transform in my registration/optimization
>  > parameters ?
>  >
>  > 3. Does resampling take care of the resolution difference ?
>  >
>  > 4. Is it necessary that I set the initial scaling to 1.428 to compensate
>  > for the resolution ?
>  >
>  > 5. Dicom images have Image Patient Postion tag. Is this the information
>  > that goes into the  "Origin" variable ?  If not, where do I get the
>  > image origin information ?
>  >
>  > I have read a lot of literature on resampling and all the pdfs and ppts
>  > on ITK, but I fail to understand these.
>  >
>  > Please help !
>  >
>  > Emma
>  >
>  > ------------------------------------------------------------------------
>  > Get the Yahoo! toolbar and be alerted to new email
>  > 
> wherever 
>  
> 
>  > you're surfing.
>  >
>  >
>  > ------------------------------------------------------------------------
>  >
>  > _______________________________________________
>  > Insight-users mailing list
>  > Insight-users at itk.org
>  > http://www.itk.org/mailman/listinfo/insight-users
> 
> 
> ------------------------------------------------------------------------
> Got a little couch potato?
> Check out fun summer activities for kids. 
> 


       
---------------------------------
Be a better friend, newshound, and know-it-all with Yahoo! Mobile.  Try it now.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://public.kitware.com/pipermail/insight-users/attachments/20080221/7ca5f323/attachment-0001.html


More information about the Insight-users mailing list