<div dir="ltr">HI,<div>I have always the same issue with my demons registration code, "Inputs do not occupy the same physical space".</div><div> the fixed image is a Ct  512x512 dicom and the second one is a Pet  256x256 dicom.</div><div>May you help me to find out what is the problem?</div><div>Thank you very much.<div><br></div><div><br></div><div><br></div><div>#include "itkImageFileReader.h"</div><div>#include "itkImageFileWriter.h"</div><div>#include <time.h></div><div>#include <fstream></div><div>#include "itkSymmetricForcesDemonsRegistrationFilter.h"</div><div>#include "itkHistogramMatchingImageFilter.h"</div><div>#include "itkCastImageFilter.h"</div><div>#include "itkWarpImageFilter.h"</div><div><br></div><div>#include "itkGDCMImageIO.h"</div><div>#include "itkShrinkImageFilter.h"</div><div><br></div><div><br></div><div><br></div><div>  class CommandIterationUpdate : public itk::Command</div><div>  {</div><div>  public:</div><div>    typedef  CommandIterationUpdate                     Self;</div><div>    typedef  itk::Command                               Superclass;</div><div>    typedef  itk::SmartPointer<CommandIterationUpdate>  Pointer;</div><div>    itkNewMacro( CommandIterationUpdate );</div><div>  protected:</div><div>    CommandIterationUpdate() {};</div><div><br></div><div>    typedef itk::Image< double, 2 >            InternalImageType;</div><div>    typedef itk::Vector< double, 2 >           VectorPixelType;</div><div>    typedef itk::Image<  VectorPixelType, 2 > DisplacementFieldType;</div><div><br></div><div>    typedef itk::SymmetricForcesDemonsRegistrationFilter<</div><div>                                InternalImageType,</div><div>                                InternalImageType,</div><div>                                DisplacementFieldType>   RegistrationFilterType;</div><div><br></div><div>  public:</div><div><br></div><div>    void Execute(itk::Object *caller, const itk::EventObject & event) ITK_OVERRIDE</div><div>      {</div><div>        Execute( (const itk::Object *)caller, event);</div><div>      }</div><div><br></div><div>    void Execute(const itk::Object * object, const itk::EventObject & event) ITK_OVERRIDE</div><div>      {</div><div>         const RegistrationFilterType * filter = static_cast< const RegistrationFilterType * >( object );</div><div>        if( !(itk::IterationEvent().CheckEvent( &event )) )</div><div>          {</div><div>          return;</div><div>          }</div><div>        std::cout << filter->GetMetric() << std::endl;</div><div>      }</div><div>  };</div><div><br></div><div><br></div><div>int main( int argc, char *argv[] )</div><div>{</div><div>  /*if( argc < 4 )</div><div>    {</div><div>    std::cerr << "Missing Parameters " << std::endl;</div><div>    std::cerr << "Usage: " << argv[0];</div><div>    std::cerr << " fixedImageFile movingImageFile ";</div><div>    std::cerr << " outputImageFile " << std::endl;</div><div>    std::cerr << " [outputDisplacementFieldFile] " << std::endl;</div><div>    return EXIT_FAILURE;</div><div>    }*/</div><div><br></div><div><br></div><div><span class="" style="white-space:pre">   </span>time_t     now = time(0);</div><div>      struct tm  tstruct;</div><div>      char       buf[80];</div><div>      tstruct = *localtime(&now);</div><div>      strftime(buf, sizeof(buf), "%Y-%m-%d %X", &tstruct);</div><div><span class="" style="white-space:pre">     </span>  </div><div>     std::ofstream resultfile;</div><div>      char Risultati[] = "D:/Images/MetricResultDeformaion.txt";</div><div>      <span class="" style="white-space:pre">   </span>   </div><div>    resultfile.open(Risultati, std::ios::app);  </div><div>    if (resultfile.is_open())</div><div>     {</div><div>        std::cout << "File Open exists\n";</div><div>  </div><div>     }</div><div>   </div><div><span class="" style="white-space:pre">  </span>resultfile<<"Data :"<<buf<<std::endl;</div><div><br></div><div><br></div><div><br></div><div>  </div><div>  const unsigned int Dimension = 2;</div><div>  typedef unsigned int PixelType;</div><div><br></div><div>  typedef itk::Image< PixelType, Dimension >  FixedImageType;</div><div>  typedef itk::Image< PixelType, Dimension >  MovingImageType;</div><div>   </div><div>  </div><div><br></div><div>   </div><div>  typedef itk::ImageFileReader< FixedImageType  > FixedImageReaderType;</div><div>  typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;</div><div>  typedef itk::GDCMImageIO GDCMType;</div><div>  GDCMType::Pointer gdcm=GDCMType::New();</div><div>  FixedImageReaderType::Pointer fixedImageReader   = FixedImageReaderType::New();</div><div>  MovingImageReaderType::Pointer movingImageReader = MovingImageReaderType::New();</div><div><br></div><div>  char* filename1="D:/Images/def00001.dcm";</div><div>  char* filename2="D:/Images/def100001.dcm";</div><div>  fixedImageReader->SetFileName( filename1 );</div><div>  movingImageReader->SetFileName( filename2 );</div><div>  fixedImageReader->SetImageIO(gdcm);</div><div>  movingImageReader->SetImageIO(gdcm);</div><div><br></div><div>  resultfile<<filename1<<std::endl;</div><div>  resultfile<<filename2<<std::endl;</div><div><br></div><div>  fixedImageReader->Update();</div><div>  movingImageReader->Update();</div><div> </div><div>  typedef itk::ShrinkImageFilter<MovingImageType,MovingImageType> ShrinkFilterType;</div><div>ShrinkFilterType::Pointer shrinkFilter = ShrinkFilterType::New();</div><div>shrinkFilter->SetShrinkFactors( 2 );</div><div>shrinkFilter->SetInput( fixedImageReader->GetOutput() );</div><div>shrinkFilter->Update();</div><div>  </div><div>  typedef double                                      InternalPixelType;</div><div>  typedef itk::Image< InternalPixelType, Dimension > InternalImageType;</div><div>  typedef itk::CastImageFilter< FixedImageType, InternalImageType >  FixedImageCasterType;</div><div>  typedef itk::CastImageFilter< MovingImageType, InternalImageType >  MovingImageCasterType;</div><div><br></div><div>  FixedImageCasterType::Pointer fixedImageCaster = FixedImageCasterType::New();</div><div>  MovingImageCasterType::Pointer movingImageCaster= MovingImageCasterType::New();</div><div><br></div><div>  fixedImageCaster->SetInput( shrinkFilter->GetOutput() );</div><div>  movingImageCaster->SetInput( fixedImageReader->GetOutput() );</div><div><br></div><div>  fixedImageCaster->Update();</div><div>  movingImageCaster->Update();</div><div>  </div><div>  typedef itk::HistogramMatchingImageFilter< InternalImageType,InternalImageType >   MatchingFilterType;</div><div>  MatchingFilterType::Pointer matcher = MatchingFilterType::New();</div><div>    </div><div>  matcher->SetInput( movingImageCaster->GetOutput() );</div><div>  </div><div>  matcher->SetReferenceImage( fixedImageCaster->GetOutput() );</div><div>  </div><div>  </div><div>  matcher->SetNumberOfHistogramLevels(1024);</div><div>  matcher->SetNumberOfMatchPoints(7  );</div><div> </div><div>  </div><div>  matcher->ThresholdAtMeanIntensityOn();</div><div>  </div><div>  typedef itk::Vector< double, Dimension >           VectorPixelType;</div><div>  typedef itk::Image<  VectorPixelType, Dimension > DisplacementFieldType;</div><div><br></div><div><br></div><div><span class="" style="white-space:pre">      </span>  </div><div><span class="" style="white-space:pre">        </span></div><div><span class="" style="white-space:pre">   </span>  typedef itk::SymmetricForcesDemonsRegistrationFilter<InternalImageType,InternalImageType, DisplacementFieldType> RegistrationFilterType;</div><div>  RegistrationFilterType::Pointer filter = RegistrationFilterType::New();</div><div>  </div><div><br></div><div><br></div><div>  CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();</div><div>  filter->AddObserver( itk::IterationEvent(), observer );</div><div>  </div><div>  </div><div>  </div><div>  filter->SetFixedImage( fixedImageCaster->GetOutput() );</div><div>  filter->SetMovingImage( matcher->GetOutput());</div><div>  filter->SetNumberOfIterations( 256 );</div><div>  filter->SetStandardDeviations( 4.0 );</div><div> </div><div><br></div><div><br></div><div>  try</div><div>{</div><div>filter->Update();</div><div>}</div><div>catch( itk::ExceptionObject & err )</div><div>{</div><div>std::cerr << "ExceptionObject caught !" << std::endl;</div><div>std::cerr << err << std::endl;</div><div>return EXIT_FAILURE;</div><div>}</div><div>  </div><div>  </div><div>  typedef itk::WarpImageFilter<  MovingImageType, MovingImageType, DisplacementFieldType  >     WarperType;</div><div>  typedef itk::LinearInterpolateImageFunction< MovingImageType,   double >  InterpolatorType;</div><div>  WarperType::Pointer warper = WarperType::New();</div><div>  InterpolatorType::Pointer interpolator = InterpolatorType::New();</div><div>  FixedImageType::Pointer fixedImage = fixedImageReader->GetOutput();</div><div><br></div><div>  warper->SetInput( movingImageReader->GetOutput() );</div><div>  warper->SetInterpolator( interpolator );</div><div>  warper->SetOutputSpacing( fixedImage->GetSpacing() );</div><div>  warper->SetOutputOrigin( fixedImage->GetOrigin() );</div><div>  warper->SetOutputDirection( fixedImage->GetDirection() );</div><div>  </div><div>  warper->SetDisplacementField( filter->GetOutput() );</div><div>  </div><div>  warper->Update();</div><div><br></div><div><br></div><div>   </div><div>  typedef  unsigned char                           OutputPixelType;</div><div>  typedef itk::Image< OutputPixelType, Dimension > OutputImageType;</div><div>  typedef itk::CastImageFilter<</div><div>                        MovingImageType,</div><div>                        OutputImageType >          CastFilterType;</div><div>  typedef itk::ImageFileWriter< OutputImageType >  WriterType;</div><div><br></div><div>  WriterType::Pointer      writer =  WriterType::New();</div><div>  CastFilterType::Pointer  caster =  CastFilterType::New();</div><div><br></div><div><br></div><div><br></div><div>  </div><div><br></div><div><br></div><div><br></div><div>  std::cout<<"last metric value "<<filter->GetMetric()<<std::endl;</div><div>  resultfile<<"last metric value "<<filter->GetMetric()<<std::endl;</div><div>  writer->SetFileName( "D:/Images/new.dcm" );</div><div><br></div><div>  caster->SetInput( warper->GetOutput() );</div><div>  caster->Update();</div><div>    writer->SetInput( caster->GetOutput()   );</div><div><span class="" style="white-space:pre">   </span>writer->SetImageIO(gdcm);</div><div>  writer->Update();</div><div><br></div><div><br></div><div>    return EXIT_SUCCESS;</div><div>}</div><div><br></div><div><br></div><div><br></div><div><br></div><h1 class="" style="font-size:30px;margin:0px 150px 0px 0px;line-height:1.1;font-weight:normal;word-wrap:break-word;color:rgb(51,51,51);font-family:Helvetica,arial,nimbussansl,liberationsans,freesans,clean,sans-serif,'Segoe UI Emoji','Segoe UI Symbol'"><br></h1></div></div>