[ITK-users] Symmetric Demons registration

[Matt McCormick]

Hi Stefano,

The means that the images must occupy the same spatial domain.  The
domain is determined by the Origin, Spacing, Direction, and
ImageRegion Size of the Image's.  This is explained more in the ITK
Software Guide:

  http://itk.org/ITKSoftwareGuide/html/Book1/ITKSoftwareGuide-Book1ch4.html#x40-430004

HTH,
Matt

On Sun, Sep 27, 2015 at 5:56 PM, stefano serviddio
<s.serviddio at gmail.com> wrote:
> HI,
> I have always the same issue with my demons registration code, "Inputs do
> not occupy the same physical space".
>  the fixed image is a Ct  512x512 dicom and the second one is a Pet  256x256
> dicom.
> May you help me to find out what is the problem?
> Thank you very much.
>
>
>
> #include "itkImageFileReader.h"
> #include "itkImageFileWriter.h"
> #include <time.h>
> #include <fstream>
> #include "itkSymmetricForcesDemonsRegistrationFilter.h"
> #include "itkHistogramMatchingImageFilter.h"
> #include "itkCastImageFilter.h"
> #include "itkWarpImageFilter.h"
>
> #include "itkGDCMImageIO.h"
> #include "itkShrinkImageFilter.h"
>
>
>
>   class CommandIterationUpdate : public itk::Command
>   {
>   public:
>     typedef  CommandIterationUpdate                     Self;
>     typedef  itk::Command                               Superclass;
>     typedef  itk::SmartPointer<CommandIterationUpdate>  Pointer;
>     itkNewMacro( CommandIterationUpdate );
>   protected:
>     CommandIterationUpdate() {};
>
>     typedef itk::Image< double, 2 >            InternalImageType;
>     typedef itk::Vector< double, 2 >           VectorPixelType;
>     typedef itk::Image<  VectorPixelType, 2 > DisplacementFieldType;
>
>     typedef itk::SymmetricForcesDemonsRegistrationFilter<
>                                 InternalImageType,
>                                 InternalImageType,
>                                 DisplacementFieldType>
> RegistrationFilterType;
>
>   public:
>
>     void Execute(itk::Object *caller, const itk::EventObject & event)
> ITK_OVERRIDE
>       {
>         Execute( (const itk::Object *)caller, event);
>       }
>
>     void Execute(const itk::Object * object, const itk::EventObject & event)
> ITK_OVERRIDE
>       {
>          const RegistrationFilterType * filter = static_cast< const
> RegistrationFilterType * >( object );
>         if( !(itk::IterationEvent().CheckEvent( &event )) )
>           {
>           return;
>           }
>         std::cout << filter->GetMetric() << std::endl;
>       }
>   };
>
>
> 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 " << std::endl;
>     std::cerr << " [outputDisplacementFieldFile] " << std::endl;
>     return EXIT_FAILURE;
>     }*/
>
>
> time_t     now = time(0);
>       struct tm  tstruct;
>       char       buf[80];
>       tstruct = *localtime(&now);
>       strftime(buf, sizeof(buf), "%Y-%m-%d %X", &tstruct);
>
>      std::ofstream resultfile;
>       char Risultati[] = "D:/Images/MetricResultDeformaion.txt";
>
>     resultfile.open(Risultati, std::ios::app);
>     if (resultfile.is_open())
>      {
>         std::cout << "File Open exists\n";
>
>      }
>
> resultfile<<"Data :"<<buf<<std::endl;
>
>
>
>
>   const unsigned int Dimension = 2;
>   typedef unsigned int PixelType;
>
>   typedef itk::Image< PixelType, Dimension >  FixedImageType;
>   typedef itk::Image< PixelType, Dimension >  MovingImageType;
>
>
>
>
>   typedef itk::ImageFileReader< FixedImageType  > FixedImageReaderType;
>   typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;
>   typedef itk::GDCMImageIO GDCMType;
>   GDCMType::Pointer gdcm=GDCMType::New();
>   FixedImageReaderType::Pointer fixedImageReader   =
> FixedImageReaderType::New();
>   MovingImageReaderType::Pointer movingImageReader =
> MovingImageReaderType::New();
>
>   char* filename1="D:/Images/def00001.dcm";
>   char* filename2="D:/Images/def100001.dcm";
>   fixedImageReader->SetFileName( filename1 );
>   movingImageReader->SetFileName( filename2 );
>   fixedImageReader->SetImageIO(gdcm);
>   movingImageReader->SetImageIO(gdcm);
>
>   resultfile<<filename1<<std::endl;
>   resultfile<<filename2<<std::endl;
>
>   fixedImageReader->Update();
>   movingImageReader->Update();
>
>   typedef itk::ShrinkImageFilter<MovingImageType,MovingImageType>
> ShrinkFilterType;
> ShrinkFilterType::Pointer shrinkFilter = ShrinkFilterType::New();
> shrinkFilter->SetShrinkFactors( 2 );
> shrinkFilter->SetInput( fixedImageReader->GetOutput() );
> shrinkFilter->Update();
>
>   typedef double                                      InternalPixelType;
>   typedef itk::Image< InternalPixelType, Dimension > InternalImageType;
>   typedef itk::CastImageFilter< FixedImageType, InternalImageType >
> FixedImageCasterType;
>   typedef itk::CastImageFilter< MovingImageType, InternalImageType >
> MovingImageCasterType;
>
>   FixedImageCasterType::Pointer fixedImageCaster =
> FixedImageCasterType::New();
>   MovingImageCasterType::Pointer movingImageCaster=
> MovingImageCasterType::New();
>
>   fixedImageCaster->SetInput( shrinkFilter->GetOutput() );
>   movingImageCaster->SetInput( fixedImageReader->GetOutput() );
>
>   fixedImageCaster->Update();
>   movingImageCaster->Update();
>
>   typedef itk::HistogramMatchingImageFilter<
> InternalImageType,InternalImageType >   MatchingFilterType;
>   MatchingFilterType::Pointer matcher = MatchingFilterType::New();
>
>   matcher->SetInput( movingImageCaster->GetOutput() );
>
>   matcher->SetReferenceImage( fixedImageCaster->GetOutput() );
>
>
>   matcher->SetNumberOfHistogramLevels(1024);
>   matcher->SetNumberOfMatchPoints(7  );
>
>
>   matcher->ThresholdAtMeanIntensityOn();
>
>   typedef itk::Vector< double, Dimension >           VectorPixelType;
>   typedef itk::Image<  VectorPixelType, Dimension > DisplacementFieldType;
>
>
>
>  typedef
> itk::SymmetricForcesDemonsRegistrationFilter<InternalImageType,InternalImageType,
> DisplacementFieldType> RegistrationFilterType;
>   RegistrationFilterType::Pointer filter = RegistrationFilterType::New();
>
>
>
>   CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
>   filter->AddObserver( itk::IterationEvent(), observer );
>
>
>
>   filter->SetFixedImage( fixedImageCaster->GetOutput() );
>   filter->SetMovingImage( matcher->GetOutput());
>   filter->SetNumberOfIterations( 256 );
>   filter->SetStandardDeviations( 4.0 );
>
>
>
>   try
> {
> filter->Update();
> }
> catch( itk::ExceptionObject & err )
> {
> std::cerr << "ExceptionObject caught !" << std::endl;
> std::cerr << err << std::endl;
> return EXIT_FAILURE;
> }
>
>
>   typedef itk::WarpImageFilter<  MovingImageType, MovingImageType,
> DisplacementFieldType  >     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->SetOutputDirection( fixedImage->GetDirection() );
>
>   warper->SetDisplacementField( filter->GetOutput() );
>
>   warper->Update();
>
>
>
>   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();
>
>
>
>
>
>
>
>   std::cout<<"last metric value "<<filter->GetMetric()<<std::endl;
>   resultfile<<"last metric value "<<filter->GetMetric()<<std::endl;
>   writer->SetFileName( "D:/Images/new.dcm" );
>
>   caster->SetInput( warper->GetOutput() );
>   caster->Update();
>     writer->SetInput( caster->GetOutput()   );
> writer->SetImageIO(gdcm);
>   writer->Update();
>
>
>     return EXIT_SUCCESS;
> }
>
>
>
>
>
>
> _____________________________________
> Powered by www.kitware.com
>
> Visit other Kitware open-source projects at
> http://www.kitware.com/opensource/opensource.html
>
> Kitware offers ITK Training Courses, for more information visit:
> http://www.kitware.com/products/protraining.php
>
> Please keep messages on-topic and check the ITK FAQ at:
> http://www.itk.org/Wiki/ITK_FAQ
>
> Follow this link to subscribe/unsubscribe:
> http://public.kitware.com/mailman/listinfo/insight-users
>