<html><head><meta http-equiv="Content-Type" content="text/html charset=iso-8859-1"></head><body style="word-wrap: break-word; -webkit-nbsp-mode: space; -webkit-line-break: after-white-space; ">Hello lien,<div><br></div><div>Are you comparing 4.2.1 to 3.20? Or the v3 vs v4 registration in 4.2.1?</div><div><br></div><div>Also how does v3 compare to v4 when only a single thread is used?</div><div><br></div><div>Thanks for the additional information.</div><div><br></div><div>Brad</div><div><br><div><div>On Oct 26, 2012, at 9:12 AM, lien lee <<a href="mailto:lienlee@gmail.com">lienlee@gmail.com</a>> wrote:</div><br class="Apple-interchange-newline"><blockquote type="cite">Thanks for the informative reply, Brian. I am going to double check, and will let you know what I could get.<br><br><br><div class="gmail_quote">2012/10/25 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">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 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>
<div><div class="h5"><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 (as<br>
> attached at the bottom of this message) was created and tested on a pair 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 data,<br>
> but it just took about 12s on the same machine with almost the same matching<br>
> result.<br>
><br>
> As a newbie, I am not sure whether I did the right thing and I am trying 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) 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, I am<br>
> not so sure.<br>
><br>
> Of course, I can just stick to v3, but, I am just curious whether there 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 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 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>
> 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>
><br>
> vComposite->AddTransform(const_cast<RigidTransformType*>(vRigidRegistration->GetOutput()->Get()));<br>
> return true;<br>
> }<br>
> //<br>
> //=== End of the code =====================================<br>
><br>
><br>
</div></div>> _____________________________________<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>
</blockquote></div><br>
_____________________________________<br>Powered by <a href="http://www.kitware.com">www.kitware.com</a><br><br>Visit other Kitware open-source projects at<br><a href="http://www.kitware.com/opensource/opensource.html">http://www.kitware.com/opensource/opensource.html</a><br><br>Kitware offers ITK Training Courses, for more information visit:<br>http://www.kitware.com/products/protraining.php<br><br>Please keep messages on-topic and check the ITK FAQ at:<br>http://www.itk.org/Wiki/ITK_FAQ<br><br>Follow this link to subscribe/unsubscribe:<br>http://www.itk.org/mailman/listinfo/insight-users<br></blockquote></div><br></div></body></html>