[Insight-users] Demons 3D deformable registration problem

[Luis Ibanez]

Hi 姜华,


1) It is very normal for the Demons registration class to produce
    empty regions on the borders when you use a low value of the
    standard deviation.

    The reason is that the standard deviation values are used for
    defining a Gaussian smoothing that is applied to the deformation
    field after every iteration of the filter.

    When you use large values of standard deviation, you apply
    a stronger smoothing and therefore you obtain a deformation
    field that is more uniform across the image and that contains
    displacement vectors of smaller magnitude.

    The dark regions represent the location for which the deformation
    field can't find corresponding pixels in the moving image during
    the process of resampling.

2)  Therefore, this is not really a "problem", but a normal behavior
      of the deformable registration filter. It is up to you to fine tune
      the parameters in order to generate deformations that are
      plausible and useful.

      Note that you will benefit from visualizing the deformation field
      superimposed to the fixed image.

      You can easily produce this sort of visualization by using any
      of the following applications:

      4D Slicer:
      http://www.creatis.insa-lyon.fr/rio/vv

       Paraview:
       http://www.paraview.org


    Regards,


          Luis


------------------------------------------------------
On Wed, Dec 30, 2009 at 9:55 AM, 姜华 <huajiang0518 at gmail.com> wrote:
> Hi,luis,
>    I have several questions about ITK deformable registration.
>   First, I use the Demons for 3D deformable registration, when I set the
> parameter "SetStandardDeviations" smaller(for example 1.0),
>   the image whose size is 256*176*16 will have black regions in each slice,
> but when I set this paremeter bigger (for example 10.0),
>   only the first and the last slice have  black regions, so I do not know
> why this phenomena will happen?
>   second, since the phenomena happened, so I use the "MirrorPadImageFilter"
> to pad the image in the slice direction in order to
>   make its size 256*176*18, and set the "SetStandardDeviations"  parameter
> 10.0, then after registration,
>   I use the WarpImageFilter to warp the movingImage according to the
> deformable field, at last I use the "RegionOfInterestImageFilter" extract
> the image from second slice to the seventeenth slice. Thus I can get better
> result, but I don't know this method is whether correct.
>     In order to explain these problem, I copy the part programs in the
> email, and the result images in the inclosure.
>  Thank you for your answers as soon as possible.
>
>  With best regards!
>
>                  HuaJiang,
>
>
> The Program:(this is the part of the my program, the read and write parts
> arenot attached)
> FixedImageCasterType::Pointer fixedImageCaster   =
> FixedImageCasterType::New();
>     MovingImageCasterType::Pointer movingImageCaster =
> MovingImageCasterType::New();
>  fixedImageCaster->SetInput( fixedImage );
>     movingImageCaster->SetInput( movingImage );
> /*
>  HistogramMatchFilterType::Pointer matcher =
> HistogramMatchFilterType::New();
>  matcher->SetInput( movingImageCaster->GetOutput() );
>     matcher->SetReferenceImage( fixedImageCaster->GetOutput() );
>
>  matcher->SetNumberOfHistogramLevels( 1024 );
>     matcher->SetNumberOfMatchPoints( 7 );
>     matcher->ThresholdAtMeanIntensityOn();
>     */
>
>
>  Demons3RegistrationFilterType::Pointer filter =
> Demons3RegistrationFilterType::New();
>
>  filter->SetStandardDeviations( 10.0 );
>  MultiResRegistrationFilterType:: Pointer multires =
> MultiResRegistrationFilterType::New();
>
>  multires->SetRegistrationFilter( filter );
>  multires->SetNumberOfLevels( 3 );
>
>  unsigned int nIterations[4] = {32, 16, 4 };
>     multires->SetNumberOfIterations( nIterations );
>  DemonsCommandIterationUpdate::Pointer observer =
> DemonsCommandIterationUpdate::New();
>  filter->AddObserver(itk::IterationEvent(), observer);
>
>  CommandResolutionLevelUpdate::Pointer levelobserver =
> CommandResolutionLevelUpdate::New();
>     multires->AddObserver( itk::IterationEvent(), levelobserver );
>
>  // Pad the image.
>     unsigned long upFactor[3] = {0,0,1};
>  unsigned long lowFactor[3] = {0,0,1};
>     InternalPixelType defaultPixel = 0;
>  // define warp filter.
>  WarperType::Pointer warper = WarperType::New();
>  typedef itk::LinearInterpolateImageFunction< InternalImageType, double >
> LinearInterpolatorType;
>  LinearInterpolatorType::Pointer linearInterpolator =
> LinearInterpolatorType::New();
>
>  warper->SetInterpolator( linearInterpolator );
>     ConstantPadImageFilterType::Pointer constantPadFixedFilter;
>  ConstantPadImageFilterType::Pointer constantPadMovingFilter;
>  MirrorPadImageFilterType::Pointer mirrorPadFixedFilter;
>  MirrorPadImageFilterType::Pointer mirrorPadMovingFilter;
>  if (m_padMode == ConstantPad)
>  {
>      constantPadFixedFilter = ConstantPadImageFilterType::New();
>      constantPadMovingFilter = ConstantPadImageFilterType::New();
>
>   constantPadFixedFilter->SetInput(fixedImageCaster->GetOutput());
>   constantPadMovingFilter->SetInput(movingImageCaster->GetOutput());
>
>   constantPadFixedFilter->SetPadLowerBound(lowFactor);
>   constantPadFixedFilter->SetPadUpperBound(upFactor);
>         constantPadFixedFilter->SetConstant(defaultPixel);
>         constantPadMovingFilter->SetPadLowerBound(lowFactor);
>   constantPadMovingFilter->SetPadUpperBound(upFactor);
>         constantPadMovingFilter->SetConstant(defaultPixel);
>   multires->SetFixedImage( constantPadFixedFilter->GetOutput() );
>      multires->SetMovingImage( constantPadMovingFilter->GetOutput() );
>
>   try
>   {
>    multires->Update();
>   }
>   catch( itk::ExceptionObject & excp )
>   {
>    std::cerr << excp << std::endl;
>    return ;
>     }
>         warper->SetInput(constantPadMovingFilter->GetOutput());
>         warper->SetOutputSpacing(
> constantPadFixedFilter->GetOutput()->GetSpacing() );
>   warper->SetOutputOrigin( constantPadFixedFilter->GetOutput()->GetOrigin()
> );
>
>   warper->SetDeformationField( multires->GetOutput() );
>         warper->Update();
>  }
>  else if (m_padMode == MirrorPad)
>  {
>      mirrorPadFixedFilter = MirrorPadImageFilterType::New();
>   mirrorPadMovingFilter = MirrorPadImageFilterType::New();
>
>   mirrorPadFixedFilter->SetInput(fixedImageCaster->GetOutput());
>   mirrorPadMovingFilter->SetInput(movingImageCaster->GetOutput());
>
>   mirrorPadFixedFilter->SetPadLowerBound(lowFactor);
>   mirrorPadFixedFilter->SetPadUpperBound(upFactor);
>
>         mirrorPadMovingFilter->SetPadLowerBound(lowFactor);
>   mirrorPadMovingFilter->SetPadUpperBound(upFactor);
>
>   multires->SetFixedImage( mirrorPadFixedFilter->GetOutput() );
>   multires->SetMovingImage( mirrorPadMovingFilter->GetOutput() );
>
>   try
>   {
>    multires->Update();
>   }
>   catch( itk::ExceptionObject & excp )
>   {
>    std::cerr << excp << std::endl;
>    return ;
>   }
>
>   warper->SetInput(mirrorPadMovingFilter->GetOutput());
>         warper->SetOutputSpacing(
> mirrorPadFixedFilter->GetOutput()->GetSpacing() );
>   warper->SetOutputOrigin( mirrorPadFixedFilter->GetOutput()->GetOrigin() );
>
>   warper->SetDeformationField( multires->GetOutput() );
>         warper->Update();
>  }
>
>     // extract the image.
>     RegionOfInterestImageFilterType::Pointer regionExtractFitler
>   = RegionOfInterestImageFilterType::New();
>
>     regionExtractFitler->SetInput(warper->GetOutput());
>  InternalImageType::RegionType desiredRegion ;
>  InternalImageType::SizeType size =
> fixedImage->GetLargestPossibleRegion().GetSize();
>  InternalImageType::IndexType index;
>  index[0] = 0;
>  index[1] = 0;
>  index[2] = 0;
>     desiredRegion.SetSize(size);
>  desiredRegion.SetIndex(index);
>  regionExtractFitler->SetRegionOfInterest(desiredRegion);
>  InternalImageCasterType::Pointer internalCaster =
> InternalImageCasterType::New();
>  internalCaster->SetInput(regionExtractFitler->GetOutput());
>     internalCaster->Update();
>
> _____________________________________
> 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.html
>
> 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://www.itk.org/mailman/listinfo/insight-users
>
>