[ITK] Image Registration Initialization
Matt McCormick
matt.mccormick at kitware.com
Thu Oct 5 16:46:38 EDT 2017
Hi Brian,
Thanks for the investigation and sharing the results.
It seems that adding an option to change the smoothing filter to
itk::SmoothingRecursiveGuassianImageFilter or to the even faster
SmothingRecursiveYvvGaussianImageFilter [1] could help.
We will also be discussing performance in general [2] -- your
participation is welcome.
Thanks,
Matt
[1] https://itk.org/Doxygen/html/classitk_1_1SmoothingRecursiveYvvGaussianImageFilter.html
[2] https://discourse.itk.org/t/next-itk-october-performance-confab/268/1
On Wed, Sep 27, 2017 at 12:44 PM, Brian Eastwood
<brian at mbfbioscience.com> wrote:
> Hi Brad,
>
>
>
> Thank you for your suggestions. The performance profiling I did was rather
> rudimentary, just inserting chrono timings in ITK code. I created a small
> example program to illustrate the situation, but I’ve seen this in many
> cases of multiresolution registration, especially if the coarse levels of
> the registration require a lot of smoothing.
>
>
>
> In my example program (below), I’m performing 3D-to-3D registration using
> itk::CenteredEuler3DTransform, itk::CorrelationImageToImageMetricv4, and
> itk::RegularStepGradientDescentOptimizerv4. My fixed image is a 3D mouse
> brain image from the Allen Institute, with size [456, 320, 528] and 25 um
> isotropic spacing; my moving image is just a synthetically rotated and
> translated version of the same image. I run registration with a single
> resolution level, smoothing the image with a standard deviation of 100 um
> and shrinking by a factor of 8.
>
>
>
> Sample timings from this program show that the smoothing that occurs in
> itk::ImageRegistrationMethodv4::InitializeRegistrationAtEachLevel()
> dominates the level initialization (15 s), the moving image’s gradient image
> computation dominates the metric initialization (4 s), and the registration
> itself runs pretty speedily (1.6 s). (The chrono timing part is not in my
> example code, so I’ve indicated which method the timings have been added
> to.)
>
>
>
> Fixed image:
>
> size [456, 320, 528]
>
> spacing [25, 25, 25]
>
> Initial transform:
>
> parameters [0, 0, 0, 5695.034988478505, 3726.9738602938187,
> 7309.779123793645, 279.5341527061055, -247.51563557593, -1.3564556700512185]
>
> center [5695.03, 3726.97, 7309.78]
>
> ImageRegistrationMethodv4::GenerateData() level: 0
>
> InitializeRegistrationAtEachLevel() smoothing: 15.4068
>
> this->InitializeRegistrationAtEachLevel(): 15.4224
>
> ImageToImageMetricv4::ComputeMovingImageGradientFilterImage(): 3.9064
>
> this->m_Metric->Initialize(): 3.9064
>
> this->m_Optimizer->StartOptimization(): 1.62503
>
> Final: -0.987401 @ [0.0007390244255403159, 0.03171806411143246,
> 0.25914684071268346, 5693.978269101054, 3725.58761384107, 7309.850962050369,
> 278.47743332865645, -248.90188202867876, -1.2846174133290214]
>
>
>
> I fully understand that smoothing takes time and that it’s a critically
> important part of multiresolution registration. My example program will not
> converge with lower smoothing values. It’s just that there are situations
> where the smoothed images could be faster to compute another way or may be
> readily available from elsewhere, so I was hoping there would be a way to
> exert more control over the multiresolution component of the framework. In
> fact, switching to itk::SmoothingRecursiveGaussianImageFilter reduces the
> initial smoothing time significantly (1.4 s), but it has different features
> that may have minor knock-on effects to ImageRegistrtaionMethodv4. I may
> also take your suggestion to separately configure and run each level, but as
> I understand it there will always be some smoothing applied to the fixed and
> moving images in itk::ImageRegistrationMethodv4, no matter how those images
> have been preprocessed upstream of registration.
>
>
>
> Thank you,
>
> Brian Eastwood
>
>
>
>
>
> /**
>
> * \file TestRegistrationv4.cpp
>
> * \brief ITKv4 registration test.
>
> * \author Brian Eastwood
>
> * \date 9/26/2017
>
> *
>
> * Example using ITK registration framework. This is a place to try out
> several
>
> * different registration procedures on different data sets.
>
> */
>
> #include <iostream>
>
>
>
> #include "itkCenteredEuler3DTransform.h"
>
> #include "itkCenteredTransformInitializer.h"
>
> #include "itkCommand.h"
>
> #include "itkCorrelationImageToImageMetricv4.h"
>
> #include "itkImage.h"
>
> #include "itkImageFileReader.h"
>
> #include "itkImageRegistrationMethodv4.h"
>
> #include "itkLinearInterpolateImageFunction.h"
>
> #include "itkObjectToObjectOptimizerBase.h"
>
> #include "itkRegularStepGradientDescentOptimizerv4.h"
>
> #include "itkVersorRigid3DTransform.h"
>
>
>
>
>
> /**
>
> * 3D-to-3D volume registration.
>
> */
>
> void testVolumeRegistration(int argc, char* argv[])
>
> {
>
> // define image types
>
> const unsigned int Dimension = 3;
>
> typedef float PixelType;
>
> typedef itk::Image< PixelType, Dimension > ImageType;
>
> typedef itk::ImageFileReader< ImageType > ReaderType;
>
>
>
> // define registration types
>
> typedef itk::CenteredEuler3DTransform<double> TransformType;
>
> typedef itk::CenteredTransformInitializer< TransformType,
> ImageType, ImageType > InitializerType;
>
> typedef itk::CorrelationImageToImageMetricv4< ImageType,
> ImageType > MetricType;
>
> typedef itk::RegularStepGradientDescentOptimizerv4< double >
> OptimizerType;
>
> typedef itk::ImageRegistrationMethodv4< ImageType, ImageType
>> RegistrationType;
>
>
>
> if (argc < 3)
>
> {
>
> std::cerr << "Usage " << argv[0] << "
> fixedImage movingImage" << std::endl;
>
> exit(1);
>
> }
>
>
>
> // create image pre-processing pipeline
>
> std::string fixedIn(argv[1]);
>
> std::string movingIn(argv[2]);
>
> ReaderType::Pointer fixedReader = ReaderType::New();
>
> fixedReader->SetFileName(fixedIn);
>
> fixedReader->Update();
>
> ReaderType::Pointer movingReader = ReaderType::New();
>
> movingReader->SetFileName(movingIn);
>
> movingReader->Update();
>
> std::cout << "Fixed image:" << std::endl;
>
> std::cout << " size " <<
> fixedReader->GetOutput()->GetLargestPossibleRegion().GetSize() << std::endl;
>
> std::cout << " spacing " <<
> fixedReader->GetOutput()->GetSpacing() << std::endl;
>
>
>
> // initialize transform parameters
>
> TransformType::Pointer transform = TransformType::New();
>
> InitializerType::Pointer initializer =
> InitializerType::New();
>
> initializer->SetFixedImage(fixedReader->GetOutput());
>
> initializer->SetMovingImage(movingReader->GetOutput());
>
> initializer->SetTransform(transform);
>
> initializer->MomentsOn();
>
> initializer->InitializeTransform();
>
> std::cout << "Initial transform:" << std::endl;
>
> std::cout << " parameters " << transform->GetParameters()
> << std::endl;
>
> std::cout << " center " << transform->GetCenter() <<
> std::endl;
>
>
>
> // create registration objects
>
> OptimizerType::Pointer optimizer = OptimizerType::New();
>
> MetricType::Pointer metric = MetricType::New();
>
> RegistrationType::Pointer registration =
> RegistrationType::New();
>
> registration->SetFixedImage(fixedReader->GetOutput());
>
> registration->SetMovingImage(movingReader->GetOutput());
>
> registration->SetOptimizer(optimizer);
>
> registration->SetMetric(metric);
>
> registration->SetInitialTransform(transform);
>
> registration->InPlaceOn();
>
>
>
> // configure resolution levels
>
> itk::SizeValueType levels = 1;
>
> ImageType::SpacingType spacing =
> fixedReader->GetOutput()->GetSpacing();
>
> RegistrationType::SmoothingSigmasArrayType sigmas(levels);
>
> RegistrationType::ShrinkFactorsArrayType shrink(levels);
>
> sigmas[0] = spacing[0] * 4;
>
> shrink[0] = 8;
>
> registration->SetNumberOfLevels(levels);
>
> registration->SetSmoothingSigmasPerLevel(sigmas);
>
>
> registration->SetSmoothingSigmasAreSpecifiedInPhysicalUnits(true);
>
> registration->SetShrinkFactorsPerLevel(shrink);
>
>
>
> // configure optimizer
>
> OptimizerType::ScalesType
> scales(transform->GetNumberOfParameters());
>
> scales.fill(1e-4);
>
> scales[0] = scales[1] = scales[2] = 1;
>
> optimizer->SetScales(scales);
>
> optimizer->SetLearningRate(1);
>
> optimizer->SetMinimumStepLength(1e-3);
>
> optimizer->SetMinimumConvergenceValue(1e-4);
>
> optimizer->SetNumberOfIterations(60);
>
>
>
> // run registration
>
> try
>
> {
>
> registration->Update();
>
> }
>
> catch (itk::ExceptionObject& err)
>
> {
>
> std::cerr << "Exception occurred during
> registration:\n" << err << std::endl;
>
> }
>
>
>
> // print result of registration
>
> std::cout << "Final: " <<
> optimizer->GetCurrentMetricValue() << " @ " << transform->GetParameters() <<
> std::endl;
>
> }
>
>
>
> int main(int argc, char* argv[])
>
> {
>
> testVolumeRegistration(argc, argv);
>
> return 0;
>
> }
>
>
>
>
>
>
>
> From: Lowekamp, Bradley (NIH/NLM/LHC) [C] [mailto:blowekamp at mail.nih.gov]
> Sent: Monday, September 25, 2017 10:02 AM
> To: Brian Eastwood <brian at mbfbioscience.com>; community at itk.org
> Cc: Nick Tustison <ntustison at gmail.com>
> Subject: Re: [ITK] Image Registration Initialization
>
>
>
> Hello,
>
>
>
> Did you do any specific profiling to determine these bottlenecks? If so
> could you please share the profiling code along with the timing you got?
>
>
>
> Also you did not provide any details about the situation where you found
> these performance issues. What transform you are using? What is the size of
> the image? What optimizer are you using?
>
>
>
> With the v4 framework, in some situations the initialization, and iterations
> are slower, however for many registration setups with a large number of
> parameters good results can be achieved with fewer iterations. Additionally,
> when the framework is combined with the parameters estimators the results
> are more robust. And with the built in multi-resolution and metric sampling
> efficient results can be achieved by using all these features of the
> framework.
>
>
>
> Based on your comments, I looked through the itk::ImageRegistrationMethodv4.
> I see that it is using the itk::DiscreteGaussianImageFilter. The
> itk::SmoothingRecursiveGuassianImageFilter is generally known to be more
> efficient. I would suggest hacking the ImageRegistrationMethodv4 to use the
> itk::SmoothingRecursiveGaussianImageFilter, and determine if there is a
> performance improvement.
>
>
>
> You cannot swap the fixed and moving inputs images at each resolution scale
> iteration. However, you could take control of scale iterations yourself, by
> just running the ImageRegistrationMethod at one level, where you choose the
> appropriate scale for setting the fixed and moving input images.
>
>
>
> HTH,
>
> Brad
>
>
>
>
>
> From: Brian Eastwood <brian at mbfbioscience.com>
> Date: Friday, September 22, 2017 at 4:55 PM
> To: "community at itk.org" <community at itk.org>
> Subject: [ITK] Image Registration Initialization
>
>
>
> Hi Folks,
>
>
>
> I really like the changes implemented for ITK v4 registration. The composite
> transform makes multistage registration much easier to implement. Nice work!
>
>
>
> Is there a way to explicitly set the fixed and moving images to be used at
> each level of a multiresolution registration? How about the gradient images
> used by the image metric? I see you can set gradient filters on the metric,
> but what if you already have the gradients precomputed?
>
>
>
> I get the feeling that v4 registration runs slower than similar v3
> registration, but I think it’s mainly down to the initialization where the
> fixed and moving images are smoothed and the metric computes the image
> gradients.
>
>
>
> It seems like there could be some real advantages to using precomputed
> smoothed and gradient images. For example, the smoothing could be
> implemented with a cascade of Gaussians where the standard deviation used at
> a level k is given by s^2 = s_k^2 - s_{k-1}^2. Or the smoothing could be
> implemented on the GPU. Or we could be reusing a resolution level in a
> different registration stage, perhaps optimizing a different set of
> transform parameters.
>
>
>
> Have I missed something in the API where it is possible to explicitly set
> the images and gradients to use at each level of the registration? I could
> probably create an itk::Command that swaps the fixed and moving images
> before each registration iteration, but is there a way to disable the
> smoothing step?
>
>
>
> Regards,
>
> Brian
>
>
> _______________________________________________
> Community mailing list
> Community at itk.org
> http://public.kitware.com/mailman/listinfo/community
>
More information about the Community
mailing list