[ITK] [ITK-users] Symmetric Demons registration

[stefano serviddio]

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;
}
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/community/attachments/20150927/9a8d89b4/attachment-0001.html>
-------------- next part --------------
_____________________________________
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