Hi all,<br><br>As a starting point with v4, a simple rigid transform image registration (as attached at the bottom of this message) was created and tested on a pair of 256x256x187 and 256x256x229 images. It takes about 90s.<br>
<br>I had an same registration based on v3, and tested against the same data, but it just took about 12s on the same machine with almost the same matching result.<br><br>As a newbie, I am not sure whether I did the right thing and I am trying to understand more about v4. By debugging into the v4 code, I noticed that:<br>
1. the new virtual domain (same as the fixed image in my case) introduces one more layer which needs some extra computation.<br>2. the transformation on a point was done through two Transformxxx() operations by two transform instances in ImageToImageMetricv4::m_CompositeTransform, although one of which is actually an identity transform.<br>
and, I am guessing maybe they are reasons for more computing time, but, I am not so sure.<br><br>Of course, I can just stick to v3, but, I am just curious whether there are some ways that I can avoid those extra computations with v4.<br>
<br><br>//=== Start of the code ===================================<br>//<br>bool<br>RigidTransform(itk::CompositeTransform<double,3>::Pointer vComposite, <br> ImageType const& vFixImage, ImageType const& vMovImage)<br>
{<br> //- The Euler transform is a rotation and translation about a center, so we<br> // need to find the rotation center.<br> //<br> typedef itk::Euler3DTransform<double> RigidTransformType;<br> RigidTransformType::Pointer vRigid = RigidTransformType::New();<br>
typedef itk::CenteredTransformInitializer< RigidTransformType,<br> ImageType, <br> ImageType > InitializerType;<br>
InitializerType::Pointer Initializer = InitializerType::New();<br> Initializer->SetTransform(vRigid);<br> Initializer->SetFixedImage(&vFixImage);<br> Initializer->SetMovingImage(&vMovImage);<br>
Initializer->GeometryOn();<br> Initializer->InitializeTransform();<br><br> vComposite->AddTransform(vRigid);<br><br> //- Metric<br> //<br> typedef itk::MattesMutualInformationImageToImageMetricv4<ImageType, ImageType> MetricType;<br>
MetricType::Pointer vMetric = MetricType::New();<br> vMetric->SetNumberOfHistogramBins(32);<br> vMetric->SetUseFixedImageGradientFilter(false);<br><br> //- Optimizer<br> //<br> typedef itk::GradientDescentOptimizerv4 OptimizerType;<br>
OptimizerType::Pointer vOptimizer = OptimizerType::New();<br> vOptimizer->SetNumberOfIterations( 50 );<br> vOptimizer->SetDoEstimateLearningRateOnce( true );<br> vOptimizer->SetMinimumConvergenceValue( 1e-6 );<br>
vOptimizer->SetConvergenceWindowSize( 5 );<br> vOptimizer->SetMaximumStepSizeInPhysicalUnits( 0.5 );<br><br> //- Scale estimator<br> //<br> itk::OptimizerParameterScalesEstimator::Pointer vScalesEstimator;<br>
typedef itk::RegistrationParameterScalesFromJacobian<MetricType> JacobianScalesEstimatorType;<br> {<br> JacobianScalesEstimatorType::Pointer vJacobianScalesEstimator<br> = JacobianScalesEstimatorType::New();<br>
vJacobianScalesEstimator->SetMetric(vMetric);<br> vJacobianScalesEstimator->SetTransformForward(true);<br> vScalesEstimator = vJacobianScalesEstimator;<br> }<br> vOptimizer->SetScalesEstimator(vScalesEstimator);<br>
vOptimizer->SetDoEstimateScales(true);<br><br> //- The RegistrationMethod class coordinates the registration operation.<br> // It needs all the pieces that come together to perform the registration<br> // operation.<br>
//<br> typedef itk::ImageRegistrationMethodv4<ImageType, ImageType, itk::Euler3DTransform<double>> RigidRegistrationType;<br> RigidRegistrationType::Pointer vRigidRegistration = RigidRegistrationType::New();<br>
vRigidRegistration->SetOptimizer(vOptimizer);<br> vRigidRegistration->SetFixedImage(&vFixImage);<br> vRigidRegistration->SetMovingImage(&vMovImage);<br> vRigidRegistration->SetMovingInitialTransform(vRigid);<br>
vRigidRegistration->SetNumberOfLevels(3);<br> vRigidRegistration->SetMetric(vMetric);<br> vRigidRegistration->SetMetricSamplingStrategy(RigidRegistrationType::RANDOM);<br> vRigidRegistration->SetMetricSamplingPercentage(0.1);<br>
<br> //- Shrink the virtual domain by specified factors for each level. <br> //<br> RigidRegistrationType::ShrinkFactorsArrayType vRigidShrinkFactors;<br> vRigidShrinkFactors.SetSize( 3 );<br> vRigidShrinkFactors[0] = 4;<br>
vRigidShrinkFactors[1] = 2;<br> vRigidShrinkFactors[2] = 1;<br> vRigidRegistration->SetShrinkFactorsPerLevel( vRigidShrinkFactors );<br><br> //- Smoothing sigmas array<br> //<br> RigidRegistrationType::SmoothingSigmasArrayType vRigidSmoothingSigmas;<br>
vRigidSmoothingSigmas.SetSize(3);<br> vRigidSmoothingSigmas.Fill(0);<br> vRigidRegistration->SetSmoothingSigmasPerLevel(vRigidSmoothingSigmas);<br><br> //- Observer<br> //<br> typedef CommandIterationUpdate< RigidRegistrationType > CommandType;<br>
CommandType::Pointer observer = CommandType::New();<br> vRigidRegistration->AddObserver( itk::InitializeEvent(), observer );<br><br> try<br> {<br> std::cout << "Starting rigid registration..." << std::endl;<br>
vRigidRegistration->Update();<br> std::cout << "Rigid parameters after registration: " << std::endl<br> << vOptimizer->GetCurrentPosition() << std::endl;<br>
}<br> catch( itk::ExceptionObject &e )<br> {<br> std::cerr << "Exception caught: " << e << std::endl;<br> return false;<br> }<br> vComposite->AddTransform(const_cast<RigidTransformType*>(vRigidRegistration->GetOutput()->Get()));<br>
return true;<br>}<br>//<br>//=== End of the code =====================================<br><br>