[ITK] Image Registration Initialization
Brian Eastwood
brian at mbfbioscience.com
Wed Sep 27 12:44:06 EDT 2017
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<mailto:brian at mbfbioscience.com>>
Date: Friday, September 22, 2017 at 4:55 PM
To: "community at itk.org<mailto:community at itk.org>" <community at itk.org<mailto: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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://public.kitware.com/pipermail/community/attachments/20170927/b3f8ccfb/attachment-0001.html>
More information about the Community
mailing list