[Insight-users] multiresolution levelset registration bug

[Lyubomir Zagorchev]

Hi Luis,

I just opened the entry.

Thanks,

Lyubomir

On 10/28/06, Luis Ibanez <luis.ibanez at kitware.com> wrote:
> 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
> >
>