Thanks for all you guys inputs.<br><br>The v3 and v4 code base are both from ITK 4.2.0. The two version code I created both used 8 threads.<br><br>I do found that the v3 version code used smaller number of sample points. So I changed the v3 code (as attached at the bottom) to use the same amount of points as those of the v4 counterpart (ie., 10% of total amount of voxels) at different resolution levels. I did this through the observer with something like:<br>
if (vRegistration->GetCurrentLevel() == 0)<br> {<br> Metric->SetNumberOfSpatialSamples(18176);<br> } else if (vRegistration->GetCurrentLevel() == 1)<br> {<br> Metric->SetNumberOfSpatialSamples( 151815 );<br>
} else<br> {<br> Metric->SetNumberOfSpatialSamples( 1223322 );<br> }<br><br>To my eyes, everything looks similar (but, like I said, I am just a newbie and not so sure about my judgement). I ran this changed v3 code, and it took about 45 s. It is about half of the v4 code, which is around 90 s as reported yesterday.<br>
<br>As for plotting the convergence curve, I am not so familiar of how to use it (probably I would go over the ItkSoftwareGuide again). Seems to me, multiple curves should be plotted with the 6 parameter transform. But, anyway, I am going to try some plots.<br>
<br>BTW, the matching performance of the v4 code is a little bit better than that of the v3 code (ie., DiceValue(v4-code) = 0.821 vs. DiceValue(v3-code) = 0.803). But, I am not sure, since this is just on one data. However, I am hoping that the automatically estimated learning rate and parameter scaling in v4 would help. <br>
<br><br><br>//=== Start of v3 code ============================================<br>//<br><br>void RigidTransform(AffineType::Pointer vAffine, ImageType const& vFixImage, <br> ImageType const& vMovImage)<br>{<br>
//- Rotation center<br> typedef itk::Euler3DTransform<double> RigidType;<br> RigidType::Pointer vRigid = RigidType::New();<br> typedef itk::CenteredTransformInitializer<RigidType,<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> //- Rigid registration<br> // ...<br> typedef itk::LinearInterpolateImageFunction<<br>
ImageType,<br> double > InterpolatorType;<br> typedef itk::GradientDescentOptimizer RigidOptimizerType;<br><br> //- Instantiate major objects.<br>
RigidOptimizerType::Pointer Optimizer = RigidOptimizerType::New();<br> InterpolatorType::Pointer Interpolator = InterpolatorType::New();<br> RegistrationType::Pointer Registration = RegistrationType::New();<br>
<br> Registration->SetNumberOfLevels(3);<br><br> //- Assemble to the registration pipeline<br> Registration->SetOptimizer( Optimizer );<br> Registration->SetInterpolator( Interpolator );<br>
Registration->SetFixedImage( &vFixImage );<br> Registration->SetMovingImage( &vMovImage );<br> ImageType::RegionType FixRegion = vFixImage.GetBufferedRegion();<br> Registration->SetFixedImageRegion( FixRegion );<br>
Registration->SetInitialTransformParameters( vRigid->GetParameters() );<br> Registration->SetTransform( vRigid );<br><br> //- Optimizer<br> typedef RigidOptimizerType::ScalesType OptimizerScalesType;<br>
OptimizerScalesType vOptimizerScales( vRigid->GetNumberOfParameters() );<br> const double vTranslationScale = 1.0 / 1000;<br> vOptimizerScales[0] = 1.0;<br> vOptimizerScales[1] = 1.0;<br> vOptimizerScales[2] = 1.0;<br>
vOptimizerScales[3] = vTranslationScale;<br> vOptimizerScales[4] = vTranslationScale;<br> vOptimizerScales[5] = vTranslationScale;<br> Optimizer->SetScales( vOptimizerScales );<br> Optimizer->SetNumberOfIterations(50);<br>
Optimizer->SetLearningRate(0.3);<br> RigidOptimizerType::ParametersType vParameters = vRigid->GetParameters();<br> Optimizer->SetInitialPosition(vParameters);<br><br> //- Metric<br> //<br> typedef itk::MattesMutualInformationImageToImageMetric<<br>
ImageType,<br> ImageType > MetricType;<br> MetricType::Pointer Metric = MetricType::New();<br> vRegistration->SetMetric(Metric);<br>
Metric->SetNumberOfHistogramBins(32);<br> Metric->ReinitializeSeed( 76926294 );<br><br> //- Setup observers<br> //<br> RigidRegistrationInterfaceCommand::Pointer vInterface = RigidRegistrationInterfaceCommand::New();<br>
Registration->AddObserver(itk::IterationEvent(), vInterface);<br> OptimizerObserver::Pointer vOptimizerObserver = OptimizerObserver::New();<br> vOptimizerObserver->mRigidMetric = dynamic_cast<itk::MattesMutualInformationImageToImageMetric<ImageType, ImageType>*><br>
(Registration->GetMetric());<br> Optimizer->AddObserver(itk::StartEvent(), vOptimizerObserver);<br> Optimizer->AddObserver(itk::IterationEvent(), vOptimizerObserver);<br><br> //- Do Translation Registration<br>
try<br> {<br> Registration->Update();<br> }<br> catch( itk::ExceptionObject & err )<br> {<br> std::cerr << "ExceptionObject caught !" << std::endl;<br> std::cerr << err << std::endl;<br>
}<br> std::cout << "Rigid Registration completed" << std::endl;<br> std::cout << std::endl;<br><br> vRigid->SetParameters(Registration->GetLastTransformParameters());<br><br>
}<br><br>//<br>//=== End of v3 code =============================================<br><br><br><div class="gmail_quote">2012/10/26 brian avants <span dir="ltr"><<a href="mailto:stnava@gmail.com" target="_blank">stnava@gmail.com</a>></span><br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">just to be a bit clearer ---- there are several key points in a<br>
reasonable comparison:<br>
<br>
1. how many threads are you using?<br>
<br>
2. how many points are you using to sample the mutual information?<br>
<br>
3. what is the shape of the convergence curve showing the mutual information?<br>
<br>
4. what version of v4 are you using?<br>
<br>
point 2 is probably the most critical.<br>
<br>
for instance, the v4 registration methods have the option<br>
<br>
registration->SetMetricSamplingPercentage(0.1);<br>
<br>
which ( with the above setting ) will sample 10% of your image -----<br>
at each resolution.<br>
<br>
this is very different than v3 which , in general , uses a fixed<br>
number of samples.<br>
<br>
v4 supports that strategy too but it is not the default.<br>
<br>
so, in comparison, you may have set up your v3 metric with something like<br>
<br>
metric->SetNumberOfSamples( 20000 )<br>
<br>
which would only be about 0.1 % of the image points.<br>
<br>
so you might first try doing<br>
<br>
registration->SetMetricSamplingPercentage(0.01);<br>
<br>
to see how this affects speed.<br>
<span class="HOEnZb"><font color="#888888"><br>
<br>
brian<br>
</font></span><div class="HOEnZb"><div class="h5"><br>
<br>
<br>
<br>
On Fri, Oct 26, 2012 at 9:12 AM, lien lee <<a href="mailto:lienlee@gmail.com">lienlee@gmail.com</a>> wrote:<br>
> Thanks for the informative reply, Brian. I am going to double check, and<br>
> will let you know what I could get.<br>
><br>
><br>
><br>
> 2012/10/25 brian avants <<a href="mailto:stnava@gmail.com">stnava@gmail.com</a>><br>
>><br>
>> hi lien<br>
>><br>
>> thanks for trying out v4.<br>
>><br>
>> i would not give up on it yet as there are many ways to use the v4<br>
>> framework.<br>
>><br>
>> i don't know off hand which is the most computationally efficient<br>
>> approach --- to some extent we rely on users to help us understand<br>
>> these questions as the number of active developers, right now, is<br>
>> relatively few.<br>
>><br>
>> typically, if we have questions about efficiency , we use a profiler<br>
>> to help understand sources of additional computation.<br>
>><br>
>> michael stauffer has done the majority of profiling - and he is cc'd here.<br>
>><br>
>> from what i recall, the v4 and v3 metrics are roughly similar in terms<br>
>> of speed and per iteration cost for computing the metric,<br>
>><br>
>> assuming that you set up the metric in a similar manner.<br>
>><br>
>> so , without seeing your v3 code, we cannot understand what<br>
>> differences might exist in your v4 and v3 comparison. one obvious<br>
>> question is per iteration cost of assessing the metric and how many<br>
>> points one is using in the call<br>
>><br>
>> vRigidRegistration->SetMetricSamplingPercentage(0.1);<br>
>><br>
>> relative to your v3 setup.<br>
>><br>
>> the "extra" computation from the virtual domain and in the composite<br>
>> transform should be relatively little as the identity transform does<br>
>> not require any actual computation.<br>
>><br>
>> just a few thoughts ---- perhaps others will have additional feedback.<br>
>><br>
>> brian<br>
>><br>
>><br>
>><br>
>><br>
>> On Thu, Oct 25, 2012 at 5:20 PM, lien lee <<a href="mailto:lienlee@gmail.com">lienlee@gmail.com</a>> wrote:<br>
>> > Hi all,<br>
>> ><br>
>> > As a starting point with v4, a simple rigid transform image registration<br>
>> > (as<br>
>> > attached at the bottom of this message) was created and tested on a pair<br>
>> > of<br>
>> > 256x256x187 and 256x256x229 images. It takes about 90s.<br>
>> ><br>
>> > I had an same registration based on v3, and tested against the same<br>
>> > data,<br>
>> > but it just took about 12s on the same machine with almost the same<br>
>> > matching<br>
>> > result.<br>
>> ><br>
>> > As a newbie, I am not sure whether I did the right thing and I am trying<br>
>> > to<br>
>> > 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)<br>
>> > introduces<br>
>> > one more layer which needs some extra computation.<br>
>> > 2. the transformation on a point was done through two Transformxxx()<br>
>> > operations by two transform instances in<br>
>> > ImageToImageMetricv4::m_CompositeTransform, although one of which is<br>
>> > actually an identity transform.<br>
>> > and, I am guessing maybe they are reasons for more computing time, but,<br>
>> > I am<br>
>> > not so sure.<br>
>> ><br>
>> > Of course, I can just stick to v3, but, I am just curious whether there<br>
>> > are<br>
>> > 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<br>
>> > center, so<br>
>> > 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 ><br>
>> > 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,<br>
>> > 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><br>
>> > 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<br>
>> > operation.<br>
>> > // It needs all the pieces that come together to perform the<br>
>> > registration<br>
>> > // operation.<br>
>> > //<br>
>> > typedef itk::ImageRegistrationMethodv4<ImageType, ImageType,<br>
>> > itk::Euler3DTransform<double>> RigidRegistrationType;<br>
>> > RigidRegistrationType::Pointer vRigidRegistration =<br>
>> > 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>
>> ><br>
>> ><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<br>
>> > vRigidSmoothingSigmas;<br>
>> > vRigidSmoothingSigmas.SetSize(3);<br>
>> > vRigidSmoothingSigmas.Fill(0);<br>
>> ><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: " <<<br>
>> > 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>
>> ><br>
>> ><br>
>> > vComposite->AddTransform(const_cast<RigidTransformType*>(vRigidRegistration->GetOutput()->Get()));<br>
>> > return true;<br>
>> > }<br>
>> > //<br>
>> > //=== End of the code =====================================<br>
>> ><br>
>> ><br>
>> > _____________________________________<br>
>> > Powered by <a href="http://www.kitware.com" target="_blank">www.kitware.com</a><br>
>> ><br>
>> > Visit other Kitware open-source projects at<br>
>> > <a href="http://www.kitware.com/opensource/opensource.html" target="_blank">http://www.kitware.com/opensource/opensource.html</a><br>
>> ><br>
>> > Kitware offers ITK Training Courses, for more information visit:<br>
>> > <a href="http://www.kitware.com/products/protraining.php" target="_blank">http://www.kitware.com/products/protraining.php</a><br>
>> ><br>
>> > Please keep messages on-topic and check the ITK FAQ at:<br>
>> > <a href="http://www.itk.org/Wiki/ITK_FAQ" target="_blank">http://www.itk.org/Wiki/ITK_FAQ</a><br>
>> ><br>
>> > Follow this link to subscribe/unsubscribe:<br>
>> > <a href="http://www.itk.org/mailman/listinfo/insight-users" target="_blank">http://www.itk.org/mailman/listinfo/insight-users</a><br>
>> ><br>
><br>
><br>
</div></div></blockquote></div><br>