[Insight-users] Rigid registration of IBSR images

[Andriy Fedorov]

Luis,

Thank you for your reply and your attention to this.

As you noted correctly, I should use the image parameters of the Fixed image.

I have discovered just now (trying to understand, why your resampler
works and mine doesn't...), that IBSR_01 and IBSR_09 have different
grid parameters: spacing is 1x1x1.5 vs 0.9375x0.9375x1.5! I
initialized the parameters from the Moving image, and that caused the
illusion of a translation error in the XY plane!

This took me too long to resolve, but I am happy, I learned a very good lesson!

Thank you very very much!

Fedorov



On 7 Apr 2007 19:25:58 -0400, Luis Ibanez <luis.ibanez at kitware.com> wrote:
>
> Hi Andriy,
>
>
> Thanks for your detailed description of the problem.
>
> (That's how a reproducible experiment should be described !)
>
>
>
>       We just rerun the registration and tried
>       your resampling code.
>
>
>             Both of them work fine.
>
>
>
> You simply get the impression that your resampled image is
> different from the image resampled by the registration
> program, because you are setting the Default pixel value to
> 0.0 while the registration program uses Default value of 100.0.
>
> If you change the default value to the same number in both
> programs, you will find that the output images are identical.
> (Well, except for the pixel type,...  you may want to unify
>   that too...).
>
>
> Note however that in your resampling program you are relying
> on the fixed and moving image to have the same origin, spacing,
> and grid parameters (number of pixels along every dimension).
> In the case of these two IBSR datasets, you can get away with
> this assumption, but in general you should use the grid parameters
> of the Fixed image. This includes:
>
>    * Origin
>    * Spacing
>    * Size
>    * Direction.
>
> Just for completeness we modified your resampling code in order
> top expect both the fixed and moving images as input.
>
> Please find the modified code attached to this email.
>
>
>
>     Regards,
>
>
>        Luis
>
>
>
> -------------------------
> Andriy Fedorov wrote:
> > Hello,
> >
> > I am using ImageRegistration8 example from the
> > Insight/Examples/Registration. I only modified it to save the final
> > transform by adding these lines:
> >
> >  itk::TransformFileWriter::Pointer trWriter =
> > itk::TransformFileWriter::New();
> >  trWriter->SetFileName("transform.par");
> >  trWriter->SetInput(finalTransform);
> >  trWriter->Update();
> >
> > The problem I am having is that if I apply this saved transform to the
> > moving image, the resulting image is different from the registered
> > image, produced by ImageRegistration8. Looks like translation is
> > incorrect, possibly incorrect translation vector direction. I did
> > check the transform as it is written and read back with
> > Print(std::cout), and it appears to be exactly the same.
> >
> > I suspect something is weird or wrong with the data I am using, or
> > with the code handling this particular kind of data, because this same
> > code and resampler work fine with the BrainWeb test images (referenced
> > in ITK guide).
> >
> > Can anybody help me what could be the problem?
> >
> > I am using IBSR "New New 18 T1 MRI" dataset, Analyze format, it can be
> > downloaded here http://www.cma.mgh.harvard.edu/ibsr/data.html.
> > Specifically, I am using image IBSR_09_ana_padded_masked.hdr as the
> > reference, and IBSR_01_ana_padded_masked.hdr as the floating image.
> >
> > I attach the resampler code that I am using.
> >
> > I tried different things to figure this out, including using
> > OrientedImage instead of just Image class, converting data from
> > Analyze to Metaheader with RPI orientation passing it through
> > OrientImage, but none seemed to work...
> >
> > I would really appreciate any help or suggestion.
> >
> > Andriy Fedorov
> >
> >
> > ------------------------------------------------------------------------
> >
> > #include "itkTransformFileReader.h"
> > #include "itkTransformFileWriter.h"
> > #include "itkLinearInterpolateImageFunction.h"
> > #include "itkNearestNeighborInterpolateImageFunction.h"
> > #include "itkImage.h"
> > #include "itkVersorRigid3DTransform.h"
> > #include "itkImageFileReader.h"
> > #include "itkImageFileWriter.h"
> > #include "itkResampleImageFilter.h"
> >
> >
> >
> > int main( int argc, char *argv[] )
> > {
> >
> >   const    unsigned int    Dimension = 3;
> >   typedef  float           PixelType;
> >
> >   typedef itk::Image< PixelType, Dimension >  ImageType;
> >   typedef itk::ImageFileReader< ImageType> ReaderType;
> >   typedef itk::ImageFileWriter< ImageType> WriterType;
> >   typedef itk::VersorRigid3DTransform< double > VersorRigid3DTransformType;
> >   typedef itk::TransformFileReader TransformReader;
> >
> >   typedef itk:: LinearInterpolateImageFunction< ImageType, double >
> >     LinearInterpolatorType;
> >   typedef itk::ResampleImageFilter<ImageType,ImageType> ResamplerType;
> >
> >   if( argc != 4 ){
> >     std::cerr << "Usage: " << argv[0];
> >     std::cerr << " movingImageFile "; // argv[1]
> >     std::cerr << " transformFile ";   // argv[2]
> >     std::cerr << " outputImagefile";  // argv[3]
> >     std::cerr << std::endl;
> >     return -1;
> >   }
> >
> >   ImageType::Pointer inputImage;
> >
> >   ReaderType::Pointer imageReader = ReaderType::New();
> >   imageReader->SetFileName(argv[1]);
> >
> >   WriterType::Pointer imageWriter = WriterType::New();
> >   imageWriter->SetFileName(argv[3]);
> >
> >   TransformReader::Pointer transformReader = TransformReader::New();
> >   transformReader->SetFileName(argv[2]);
> >
> >   ResamplerType::Pointer resampler = ResamplerType::New();
> >
> >   try{
> >     imageReader->Update();
> >     transformReader->Update();
> >   } catch(itk::ExceptionObject &e){
> >     std::cerr << "Failed to read file or transform: " << e << std::endl;
> >     return -1;
> >   }
> >
> >   inputImage = imageReader->GetOutput();
> >   resampler->SetInput(inputImage);
> >
> >   typedef itk::TransformFileReader::TransformListType* TransformListType;
> >   TransformListType transforms = transformReader->GetTransformList();
> >   itk::TransformFileReader::TransformListType::const_iterator
> >     it = transforms->begin();
> >
> >   if (!strcmp((*it)->GetNameOfClass(), "VersorRigid3DTransform")) {
> >     VersorRigid3DTransformType::Pointer  versorrigid_read =
> >         static_cast<VersorRigid3DTransformType*>((*it).GetPointer());
> >     resampler->SetTransform( versorrigid_read );
> >     std::cout << "Applying this transform: " << versorrigid_read << std::endl;
> >   } else {
> >     std::cerr << "Transform " << (*it)->GetNameOfClass()
> >       << " is not supported" << std::endl;
> >     return 1;
> >   }
> >
> >   LinearInterpolatorType::Pointer linearInterpolator =
> >       LinearInterpolatorType::New();
> >   resampler->SetInterpolator(linearInterpolator);
> >
> >   resampler->SetSize(inputImage->GetLargestPossibleRegion().GetSize());
> >   resampler->SetOutputOrigin(inputImage->GetOrigin());
> >   resampler->SetOutputSpacing(inputImage->GetSpacing());
> >   resampler->SetDefaultPixelValue(0.0);
> >
> >   imageWriter->SetInput(resampler->GetOutput());
> >
> >   try{
> >     imageWriter->Update();
> >   } catch(itk::ExceptionObject &e){
> >     std::cerr << "Failed to save the resulting image: " << e << std::endl;
> >     return -1;
> >   }
> >
> >   return 0;
> > }
> >
> >
> >
> > ------------------------------------------------------------------------
> >
> > _______________________________________________
> > Insight-users mailing list
> > Insight-users at itk.org
> > http://www.itk.org/mailman/listinfo/insight-users
>
>