[Insight-users] multiresolution levelset registration bug

[Luis Ibanez]

Hi Lyubomir,

Thanks for pointing this out.

Could you please open a bug entry for this problem in the
bug tracker ?

     http://public.kitware.com/Bug/index.php


Please let us know if you find any problem creating your
account, of entering the bug report.


   Thanks


     Luis


----------------------------
Lyubomir Zagorchev wrote:
> The LevelSetMotionRegistrationFunction generates an error
> if a MultiResolutionPDEDeformableRegistration is used with
> a LevelSetMotionRegistrationFilter. This is because the
> m_MovingImageSmoothingFilter in
> 
> template <class TFixedImage, class TMovingImage,
> class TDeformationField>
> void LevelSetMotionRegistrationFunction<
> TFixedImage,TMovingImage,TDeformationField>
> ::InitializeIteration()
> 
> does not get updated correctly in the transition from
> coarse to fine resolution. As a result, the multiresolution
> registration with the levelset filter works fine at the
> coarsest resolution and crashes erratically at higher
> resolutions.
> 
> Could you please verify this with the attached example?
> I am running ITK 2.8.1 with Visual Studio 2005.
> 
> One possible fix is to recreate the smoothing filter every
> time when InitializeIteration() is invoked (once per level)
> but I am not sure if there isn't a better approach.
> 
> Could you please let me know what is the best fix?
> 
> Thank you,
> 
> Lyubomir
> 
> 
> #include "itkImageFileReader.h"
> #include "itkImageFileWriter.h"
> #include "itkLevelSetMotionRegistrationFilter.h"
> #include "itkMultiResolutionPDEDeformableRegistration.h"
> #include "itkHistogramMatchingImageFilter.h"
> #include "itkCastImageFilter.h"
> #include "itkWarpImageFilter.h"
> #include "itkLinearInterpolateImageFunction.h"
> #include "itkExceptionObject.h"
> 
> const unsigned int Dimension = 3;
> 
>  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< float, Dimension > InternalImageType;
>    typedef itk::Vector< float, Dimension >    VectorPixelType;
>    typedef itk::Image<  VectorPixelType, Dimension > DeformationFieldType;
> 
>         typedef itk::LevelSetMotionRegistrationFilter<
>                                                             
> InternalImageType,
>                                                             
> InternalImageType,
>                                                             
> DeformationFieldType> RegistrationFilterType;
>  public:
> 
>    void Execute(itk::Object *caller, const itk::EventObject & event)
>      {
>        Execute( (const itk::Object *)caller, event);
>      }
> 
>    void Execute(const itk::Object * object, const itk::EventObject & event)
>      {
>         const RegistrationFilterType * filter =
>          dynamic_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;
>    return 1;
>    }
> 
>  typedef short PixelType;
> 
>  typedef itk::Image< PixelType, Dimension >  FixedImageType;
>  typedef itk::Image< PixelType, Dimension >  MovingImageType;
>  typedef itk::Image< float, Dimension > JacobianImageType;
> 
>  typedef itk::ImageFileReader< FixedImageType  > FixedImageReaderType;
>  typedef itk::ImageFileReader< MovingImageType > MovingImageReaderType;
> 
>  FixedImageReaderType::Pointer fixedImageReader   =
>    FixedImageReaderType::New();
>  MovingImageReaderType::Pointer movingImageReader =
>    MovingImageReaderType::New();
> 
>  fixedImageReader->SetFileName( argv[1] );
>  movingImageReader->SetFileName( argv[2] );
> 
>  typedef float 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( fixedImageReader->GetOutput() );
>  movingImageCaster->SetInput( movingImageReader->GetOutput() );
> 
>  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< float, Dimension >    VectorPixelType;
>  typedef itk::Image<  VectorPixelType, Dimension > DeformationFieldType;
> 
>     typedef itk::LevelSetMotionRegistrationFilter<
>                                                         InternalImageType,
>                                                         InternalImageType,
>                                                         
> DeformationFieldType> RegistrationFilterType;
> 
>  RegistrationFilterType::Pointer filter = RegistrationFilterType::New();
> 
>  CommandIterationUpdate::Pointer observer = CommandIterationUpdate::New();
>  filter->AddObserver( itk::IterationEvent(), observer );
> 
>  typedef itk::MultiResolutionPDEDeformableRegistration<
>                                InternalImageType,
>                                InternalImageType,
>                                DeformationFieldType >
> MultiResRegistrationFilterType;
>  MultiResRegistrationFilterType::Pointer multires =
> MultiResRegistrationFilterType::New();
>  multires->SetRegistrationFilter( filter );
>  multires->SetNumberOfLevels( 3 );
>  multires->SetFixedImage( fixedImageCaster->GetOutput() );
>  multires->SetMovingImage( matcher->GetOutput() );
> 
>  unsigned int nIterations[4] = { 10, 20, 50, 50 };
>  multires->SetNumberOfIterations( nIterations );
>  multires->Update();
> 
>  typedef itk::WarpImageFilter<
>                          MovingImageType,
>                          MovingImageType,
>                          DeformationFieldType  >     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->SetDeformationField( multires->GetOutput() );
> 
>  // Write warped image out to file
>  typedef unsigned short  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();
> 
>  writer->SetFileName( argv[3] );
> 
>  caster->SetInput( warper->GetOutput() );
>  writer->SetInput( caster->GetOutput() );
>  writer->Update();
> 
>  if( argc > 4 ) // if a fourth line argument has been provided...
>    {
> 
>    typedef itk::ImageFileWriter< DeformationFieldType > FieldWriterType;
> 
>    FieldWriterType::Pointer fieldWriter = FieldWriterType::New();
>    fieldWriter->SetFileName( argv[4] );
>    fieldWriter->SetInput( multires->GetOutput() );
> 
>    try
>      {
>      fieldWriter->Update();
>      }
>    catch ( itk::ExceptionObject e )
>      {
>      e.Print( std::cerr );
>      }
>    }
> 
>  return 0;
> }
> _______________________________________________
> Insight-users mailing list
> Insight-users at itk.org
> http://www.itk.org/mailman/listinfo/insight-users
>