[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 16:46:22 EDT 2007
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