[Insight-users] Re: My ITK effort so far to register two structural
studies using DCM and offset data
Luis Ibanez
luis.ibanez at kitware.com
Thu Oct 11 18:00:30 EDT 2007
Hi Frank,
Your email arrives just in time.
I just committed three tests to ITK for verifying the process of
reading Direction cosines from input files.
If you update your CVS checkout you will see new files in the
directory:
Insight/Testing/Code/IO
itkImageIODirection2DTest.cxx
itkImageIODirection3DTest.cxx
They just read a 2D or 3D file (of any ITK supported formats)
and verify the components of the Direction cosines with command
line arguments.
From your test, it would seem that Nifti is not passing the
direction cosines all the way up to the ITK layer. However,
that is very surprising, considering that Nifti is heavily
used in the BIRN community.
Is there a change that you could anonymize and upload one of
the Nifti files to the following page ?:
http://www.kitware.com/KitwareScripts/uploadfile.cgi
and, if you can also indicate what the expected values of the
Direction cosines are, that will help us to track any potential
bug.
Thanks
Luis
----------------------
Frank Ezekiel wrote:
> When I run the current code (looking to transform a sagittal T1 study
> into alignment with an axial T2 study) I get only background values in
> my output image. So clearly something in my code is still not
> working properly.
>
> It's probably something that an ITK expert would find easily ; sadly I
> am not one.
>
>
> When I use the movingImageReader->GetOutput()->Print( std::cout );
> function you suggested, I get output values for image direction :
>
>
> Direction:
> 0 0 -1
> 1 0 0
> 0 -1 0
>
> but NOT the actual floating point values from the Direction Cosine
> Matrix. This causes me to wonder if the nifti files I'm getting
> from ITK contain the correct DCM values.
>
>
>
> If I sent you some test data would you have a chance to take a look and
> see what I'm doing wrong?
>
>
>
> Thanks much,
>
>
> Frank
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
>
> At 01:46 PM 10/11/2007, Luis Ibanez wrote:
>
>> Hi Frank,
>>
>> What problem are you observing with the current code ?
>>
>> Please be more specific,
>>
>>
>> Thanks
>>
>>
>> Luis
>>
>>
>> ------------------------
>> Frank Ezekiel wrote:
>>
>>> Luis:
>>> I've included below my copy of main that I'm testing,
>>> attempting to take your advice on how to solve this problem. This
>>> is modified from ImageRegistration8.cxx as you suggested. There's
>>> probably just a missing piece that I don't know about...any
>>> suggestions greatly appreciated.
>>> Frank
>>>
>>>
>>> // **********************
>>> int main( int argc, char *argv[] )
>>> {
>>> if ( argc < 4 )
>>> {
>>> std::cerr << "Missing Parameters " << std::endl;
>>> std::cerr << "Usage: " << argv[0];
>>> std::cerr << " fixedImageFile movingImageFile ";
>>> std::cerr << " outputImagefile [differenceBeforeRegistration] ";
>>> std::cerr << " [differenceAfterRegistration] ";
>>> std::cerr << " [sliceBeforeRegistration] ";
>>> std::cerr << " [sliceAfterRegistration] "<< std::endl;
>>> return 1;
>>> }
>>>
>>> const unsigned int Dimension = 3;
>>> // typedef float PixelType;
>>> typedef unsigned char InputPixelType;
>>> // change 1: WORD for pixel type
>>> typedef unsigned char OutputPixelType;
>>> // PixelType
>>>
>>> typedef itk::OrientedImage< InputPixelType, Dimension >
>>> FixedImageType; // try OrientedImage rather
>>> than Image per Ibanez
>>> typedef itk::OrientedImage< InputPixelType, Dimension >
>>> MovingImageType; // This is a spot where
>>> "Oriented" was introduced
>>>
>>> typedef itk::VersorRigid3DTransform< double > TransformType;
>>>
>>>
>>>
>>> TransformType::Pointer transform = TransformType::New();
>>> transform->SetIdentity(); // frank's add-on
>>>
>>> typedef itk::ImageFileReader< FixedImageType > FixedImageReaderType;
>>> typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;
>>> FixedImageReaderType::Pointer fixedImageReader =
>>> FixedImageReaderType::New();
>>> MovingImageReaderType::Pointer movingImageReader =
>>> MovingImageReaderType::New();
>>> fixedImageReader->SetFileName( argv[1] );
>>> movingImageReader->SetFileName( argv[2] );
>>>
>>> fixedImageReader->Update();
>>>
>>> std::cout << "Fixed Image data\n" << std::endl;
>>> fixedImageReader->GetOutput()->Print( std::cout );
>>> FixedImageReaderType::SizeType size =
>>> fixedImageReader->GetOutput()->GetLargestPossibleRegion().GetSize();
>>> movingImageReader->Update();
>>> std::cout << "Moving Image data\n" << std::endl;
>>> movingImageReader->GetOutput()->Print( std::cout );
>>>
>>>
>>>
>>> typedef itk::CenteredTransformInitializer< TransformType,
>>> FixedImageType,
>>> MovingImageType
>>> >
>>> TransformInitializerType;
>>> TransformInitializerType::Pointer initializer =
>>>
>>> TransformInitializerType::New();
>>> initializer->SetTransform( transform );
>>> initializer->SetFixedImage( fixedImageReader->GetOutput() );
>>> initializer->SetMovingImage( movingImageReader->GetOutput() );
>>> initializer->MomentsOn();
>>> initializer->InitializeTransform();
>>>
>>> TransformType::MatrixType matrix = transform->GetRotationMatrix();
>>> TransformType::OffsetType offset = transform->GetOffset();
>>>
>>> matrix = transform->GetRotationMatrix();
>>> offset = transform->GetOffset();
>>> std::cout << "Matrix = " << std::endl << matrix << std::endl;
>>> std::cout << "Offset = " << std::endl << offset << std::endl;
>>>
>>>
>>> typedef itk::ResampleImageFilter<
>>> MovingImageType,
>>> FixedImageType > ResampleFilterType;
>>> //TransformType::Pointer finalTransform = TransformType::New();
>>> //finalTransform->SetCenter( transform->GetCenter() );
>>> //finalTransform->SetParameters( finalParameters );
>>> ResampleFilterType::Pointer resampler = ResampleFilterType::New();
>>> resampler->SetTransform( transform );
>>> // this is identity transform
>>> resampler->SetInput( movingImageReader->GetOutput() );
>>> FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();
>>> resampler->SetSize(
>>> fixedImage->GetLargestPossibleRegion().GetSize() );
>>> resampler->SetOutputOrigin( fixedImage->GetOrigin() );
>>> resampler->SetOutputSpacing( fixedImage->GetSpacing() );
>>> resampler->SetDefaultPixelValue( 100 );
>>>
>>>
>>> typedef itk::OrientedImage< OutputPixelType, Dimension >
>>> OutputImageType; // This is a spot where
>>> "Oriented" was introduced
>>>
>>> typedef itk::CastImageFilter<
>>> FixedImageType,
>>> OutputImageType > CastFilterType;
>>>
>>> typedef itk::ImageFileWriter< OutputImageType > WriterType;
>>>
>>> WriterType::Pointer writer = WriterType::New();
>>> CastFilterType::Pointer caster = CastFilterType::New();
>>>
>>> writer->SetFileName( argv[3] );
>>> caster->SetInput( resampler->GetOutput() );
>>> writer->SetInput( caster->GetOutput() );
>>> writer->Update();
>>> return 0;
>>>
>>> }
>>>
>>> // ***************************************
>>> At 07:06 AM 10/10/2007, you wrote:
>>>
>>>> Hi Frank,
>>>>
>>>> One easy way to see if the DCM direction is being read correctly
>>>> from your Nifti files is to print out the image information just
>>>> after reading the image.
>>>>
>>>> Please add to your code:
>>>>
>>>> reader->GetOuptut()->Print( std::cout );
>>>>
>>>> after you call
>>>>
>>>> reader->Update();
>>>>
>>>> and in the console output look for the section on the image
>>>> direction. You can then verify whether it matches the DCM direction
>>>> that you were expecting from the Nifti files.
>>>>
>>>> The itk::Image *will not* take Direction Cosines into account.
>>>> Only the itk::OrientedImage will do so.
>>>>
>>>>
>>>>
>>>> There is an open bug related to the use of itk::OrientedImages
>>>> in the ITK registration framework, but it refers to the computation
>>>> of Metric derivatives. This will not affect what you are doing,
>>>> since you should simply be running the itk::ResampleFilter.
>>>>
>>>> You *DO NOT* need to perform a registration, just a resampling.
>>>>
>>>> The itk::ResampleImageFilter, when instantiated over
>>>> itk::OrientedImages will take image orientation into account.
>>>>
>>>> The transform that you need to use is the itk::IdentityTransform.
>>>> This is an explicit ITK class, and it doesn't have any parameters,
>>>> so there is no need for you to worry about transform initialization.
>>>>
>>>>
>>>> Regards,
>>>>
>>>>
>>>> Luis
>>>>
>>>>
>>>>
>>>> ---------------------
>>>> Frank Ezekiel wrote:
>>>>
>>>>> Luis:
>>>>> I'm using 3.30 ITK, and seeing some strange behavior.
>>>>> Using ImageRegistration8 (as you suggested) as my example, I'm able
>>>>> to read in a "fixed" and "moving" image, then store the moving
>>>>> image at the resolution of the fixed.
>>>>> However, when I try to use the exact same dataset in nifti
>>>>> format files (which contain the DCM positional data) the same code
>>>>> fails. (The resulting image is sampled entirely outside the FOV
>>>>> of the source image)
>>>>> Is ITK using the DCM and transformation data from the nifti
>>>>> format version of these data even when I use only itk::Image as the
>>>>> type? I thought that functionality was limited to
>>>>> itk::OrientedImage.
>>>>> Any thoughts?
>>>>> Also, I've still not found the correct code to initialize a
>>>>> transformation that will align the MovingImage with the FixedImage
>>>>> (resolution, dimensions, AND spatial orientation)
>>>>>
>>>>> Finally, I've seen reports of problems with
>>>>> itk::OrientedImage. Are there any code samples of how to
>>>>> manually create a correctly orienting tranformation by reading the
>>>>> DCM values from the two images, then doing the matrix inversion and
>>>>> multiplication manually? The notes I saw on itk::OrientedImage
>>>>> indicated that it wasn't yet fully fixed even in version 3.40.
>>>>>
>>>>> Thanks again,
>>>>> Frank
>>>>>
>>>>> At 05:47 PM 10/8/2007, you wrote:
>>>>>
>>>>>> Hi Frank,
>>>>>>
>>>>>> Since your two images are relating to the same scanner,
>>>>>> you should be ok by simply using an IdentityTransform:
>>>>>>
>>>>>> http://www.itk.org/Insight/Doxygen/html/classitk_1_1IdentityTransform.html
>>>>>>
>>>>>>
>>>>>>
>>>>>> Note that the itk::IdentityTransform class doesn't take any
>>>>>> parameters.
>>>>>> Simply create one instance of this transform connect it to the
>>>>>> itk::ResampleImageFilter.
>>>>>>
>>>>>>
>>>>>>
>>>>>> Regards,
>>>>>>
>>>>>>
>>>>>>
>>>>>> Luis
>>>>>>
>>>>>>
>>>>>>
>>>>>> ------------------------
>>>>>> Frank Ezekiel wrote:
>>>>>>
>>>>>>> Any hints on how to initialize the transform that will be used?
>>>>>>> finalTransform->SetParameters( finalParameters );
>>>>>>>
>>>>>>> How do I set up the transform object to properly do the
>>>>>>> "difference" between the orientation and translation parameters?
>>>>>>> Thanks again,
>>>>>>> Frank
>>>>>>>
>>>>>>> At 08:15 AM 10/6/2007, you wrote:
>>>>>>>
>>>>>>>> Hi Frank,
>>>>>>>>
>>>>>>>> From your description it seems that you only need to resample
>>>>>>>> image B using the grid of image A.
>>>>>>>>
>>>>>>>> That is, you are relying on the fact that both images have been
>>>>>>>> acquired from the same scanner, and that the image origing and
>>>>>>>> direction cosines are correct.
>>>>>>>>
>>>>>>>> That should work fine.
>>>>>>>>
>>>>>>>> Just remember to use the itk::OrientedImage<> instead of the
>>>>>>>> itk::Image.
>>>>>>>>
>>>>>>>> The itk::OrientedImage will take the Direction cosines into
>>>>>>>> account when it is used by the interpolator of the Resample
>>>>>>>> ImageFilter.
>>>>>>>>
>>>>>>>>
>>>>>>>> Please read the ITK Software Guide,
>>>>>>>>
>>>>>>>> http://www.itk.org/ItkSoftwareGuide.pdf
>>>>>>>>
>>>>>>>> in particular Section 6.9.4 "Resample Image Filter" in
>>>>>>>> pdf-pages 254-448.
>>>>>>>>
>>>>>>>> This section describes in detail how to use the
>>>>>>>> ResampleImageFilter.
>>>>>>>>
>>>>>>>>
>>>>>>>> Regards,
>>>>>>>>
>>>>>>>>
>>>>>>>> Luis
>>>>>>>>
>>>>>>>>
>>>>>>>>
>>>>>>>> -----------------------
>>>>>>>> Frank Ezekiel wrote:
>>>>>>>>
>>>>>>>>> Luis:
>>>>>>>>> Thank you for the initial response. At this point,
>>>>>>>>> my goal has changed slightly. I'm looking to use ITK to write
>>>>>>>>> a function that will align MovingImage to FixedImage WITHOUT
>>>>>>>>> any intensity based registration.
>>>>>>>>> That is, I want to use ONLY the DCM (direction cosine
>>>>>>>>> matrix) and coordinate offset data included in a dicom header.
>>>>>>>>> (The reason is that my FixedImage space is from single
>>>>>>>>> voxel spectroscopy, so there is inadequate intensity
>>>>>>>>> information to use a cost function + resampling for alignment.)
>>>>>>>>> Can ITK do this? It basically just requires
>>>>>>>>> computing and applying the "difference transformation" between
>>>>>>>>> two images acquired in the same session.
>>>>>>>>>
>>>>>>>>> Thanks again for the help,
>>>>>>>>> Frank
>>>>>>>>>
>>>>>>>>>
>>>>>>>>> At 07:34 AM 8/5/2007, Luis Ibanez wrote:
>>>>>>>>>
>>>>>>>>>> Hi Frank,
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> You may want to start with the example:
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> ImageRegistration8.cxx
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> in the directory:
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> Insight/Examples/Registration
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> This example is described in detail in the
>>>>>>>>>> ITK Software Guide
>>>>>>>>>>
>>>>>>>>>> http://www.itk.org/ItkSoftwareGuide.pdf
>>>>>>>>>>
>>>>>>>>>> in the "Image Registration" chapter.
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> We *STRONGLY* recommend you to read the entire chapter
>>>>>>>>>> before you start playing with the registration. You
>>>>>>>>>> should also read the Section on image resampling since
>>>>>>>>>> they are very closely related.
>>>>>>>>>>
>>>>>>>>>> Since your images seems to be from the same modality
>>>>>>>>>> and same patient, you can probably solve this by using
>>>>>>>>>> the following components:
>>>>>>>>>>
>>>>>>>>>> 1) Mean Squares Metric
>>>>>>>>>> 2) VersorRigid3D Transform
>>>>>>>>>> 3) Linear Interpolator
>>>>>>>>>> 4) VersorRigid3DTransformOptimizer
>>>>>>>>>>
>>>>>>>>>> Note that, since you are reading this from DICOM,
>>>>>>>>>> you want to replace the "itk::Image" type in the
>>>>>>>>>> code with "itk::OrientedImage" type, in order to
>>>>>>>>>> make sure that the direction cosines of your images
>>>>>>>>>> are taken into account.
>>>>>>>>>>
>>>>>>>>>> ---
>>>>>>>>>>
>>>>>>>>>> Examples on how to read DICOM series are available
>>>>>>>>>> in the directory:
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> Insight/Examples/IO
>>>>>>>>>>
>>>>>>>>>> and are described as well in the ITK Software Guide
>>>>>>>>>> in the Chapter "Reading and Writing Images".
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> Please let us know if you find any problems
>>>>>>>>>> while setting up your registration code.
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> Regards,
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> Luis
>>>>>>>>>>
>>>>>>>>>>
>>>>>>>>>> --------------------
>>>>>>>>>> Frank Ezekiel wrote:
>>>>>>>>>>
>>>>>>>>>>> Hi:
>>>>>>>>>>> I have multiple MRI acquisitions acquired on Siemens
>>>>>>>>>>> Medspec 4T and stored in Dicom. I'm looking for the
>>>>>>>>>>> quickest and easiest way to compute the rotations and
>>>>>>>>>>> translations needed to align one dataset with the other
>>>>>>>>>>> (assuming no movement).
>>>>>>>>>>> Could someone please enumerate my options within ITK for
>>>>>>>>>>> doing this? I really just want a robust readout of angles
>>>>>>>>>>> and mm translations given Dicom 3.0 files with tags
>>>>>>>>>>> (0020,0032) and (0020,0037) filled in.
>>>>>>>>>>> Suggestions?
>>>>>>>>>>> Thanks mucho,
>>>>>>>>>>> Frank
>>>>>>>>>>> _______________________________________________
>>>>>>>>>>> Insight-users mailing list
>>>>>>>>>>> Insight-users at itk.org
>>>>>>>>>>> http://www.itk.org/mailman/listinfo/insight-users
>
>
More information about the Insight-users
mailing list