[Insight-users] Rigid registration of IBSR images

[Luis Ibanez]

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
-------------- next part --------------
A non-text attachment was scrubbed...
Name: Resampler1.cxx
Type: text/x-c++src
Size: 3457 bytes
Desc: not available
Url : http://public.kitware.com/pipermail/insight-users/attachments/20070407/2cd53d54/Resampler1.cxx