[Insight-users] Much more computation with Itk v4 as compared to v3 ?

[brian avants]

hi lien

thanks for trying out v4.

i would not give up on it yet as there are many ways to use the v4 framework.

i don't know off hand which is the most computationally efficient
approach --- to some extent we rely on users to help us understand
these questions as the number of active developers, right now, is
relatively few.

typically, if we have questions about efficiency , we use a profiler
to help understand sources of additional computation.

michael stauffer has done the majority of profiling - and he is cc'd here.

from what i recall, the v4 and v3 metrics are roughly similar in terms
of speed and per iteration cost for computing the metric,

assuming that you set up the metric in a similar manner.

so , without seeing your v3 code,  we cannot understand what
differences might exist in your v4 and v3 comparison.  one obvious
question is per iteration cost of assessing the metric and how many
points one is using in the call

vRigidRegistration->SetMetricSamplingPercentage(0.1);

relative to your v3 setup.

the "extra" computation from the virtual domain and in the composite
transform should be relatively little as the identity transform does
not require any actual computation.

just a few thoughts ---- perhaps others will have additional feedback.

brian




On Thu, Oct 25, 2012 at 5:20 PM, lien lee <lienlee at gmail.com> wrote:
> Hi all,
>
> 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.
>
> 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.
>
> 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:
> 1. the new virtual domain (same as the fixed image in my case) introduces
> one more layer which needs some extra computation.
> 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.
> and, I am guessing maybe they are reasons for more computing time, but, I am
> not so sure.
>
> 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.
>
>
> //=== Start of the code ===================================
> //
> bool
> RigidTransform(itk::CompositeTransform<double,3>::Pointer vComposite,
>     ImageType const& vFixImage, ImageType const& vMovImage)
> {
>     //- The Euler transform is a rotation and translation about a center, so
> we
>     //    need to find the rotation center.
>     //
>     typedef itk::Euler3DTransform<double> RigidTransformType;
>     RigidTransformType::Pointer vRigid = RigidTransformType::New();
>     typedef itk::CenteredTransformInitializer<    RigidTransformType,
>                                                 ImageType,
>                                                 ImageType >
> InitializerType;
>     InitializerType::Pointer Initializer = InitializerType::New();
>     Initializer->SetTransform(vRigid);
>     Initializer->SetFixedImage(&vFixImage);
>     Initializer->SetMovingImage(&vMovImage);
>     Initializer->GeometryOn();
>     Initializer->InitializeTransform();
>
>     vComposite->AddTransform(vRigid);
>
>     //- Metric
>     //
>     typedef itk::MattesMutualInformationImageToImageMetricv4<ImageType,
> ImageType> MetricType;
>     MetricType::Pointer vMetric = MetricType::New();
>     vMetric->SetNumberOfHistogramBins(32);
>     vMetric->SetUseFixedImageGradientFilter(false);
>
>     //- Optimizer
>     //
>     typedef itk::GradientDescentOptimizerv4 OptimizerType;
>     OptimizerType::Pointer vOptimizer = OptimizerType::New();
>     vOptimizer->SetNumberOfIterations( 50 );
>     vOptimizer->SetDoEstimateLearningRateOnce( true );
>     vOptimizer->SetMinimumConvergenceValue( 1e-6 );
>     vOptimizer->SetConvergenceWindowSize( 5 );
>     vOptimizer->SetMaximumStepSizeInPhysicalUnits( 0.5 );
>
>     //- Scale estimator
>     //
>     itk::OptimizerParameterScalesEstimator::Pointer vScalesEstimator;
>     typedef itk::RegistrationParameterScalesFromJacobian<MetricType>
> JacobianScalesEstimatorType;
>     {
>         JacobianScalesEstimatorType::Pointer vJacobianScalesEstimator
>           = JacobianScalesEstimatorType::New();
>         vJacobianScalesEstimator->SetMetric(vMetric);
>         vJacobianScalesEstimator->SetTransformForward(true);
>         vScalesEstimator = vJacobianScalesEstimator;
>     }
>     vOptimizer->SetScalesEstimator(vScalesEstimator);
>     vOptimizer->SetDoEstimateScales(true);
>
>     //- The RegistrationMethod class coordinates the registration operation.
>     //    It needs all the pieces that come together to perform the
> registration
>     //    operation.
>     //
>     typedef itk::ImageRegistrationMethodv4<ImageType, ImageType,
> itk::Euler3DTransform<double>> RigidRegistrationType;
>     RigidRegistrationType::Pointer vRigidRegistration =
> RigidRegistrationType::New();
>     vRigidRegistration->SetOptimizer(vOptimizer);
>     vRigidRegistration->SetFixedImage(&vFixImage);
>     vRigidRegistration->SetMovingImage(&vMovImage);
>     vRigidRegistration->SetMovingInitialTransform(vRigid);
>     vRigidRegistration->SetNumberOfLevels(3);
>     vRigidRegistration->SetMetric(vMetric);
>
> vRigidRegistration->SetMetricSamplingStrategy(RigidRegistrationType::RANDOM);
>     vRigidRegistration->SetMetricSamplingPercentage(0.1);
>
>     //- Shrink the virtual domain by specified factors for each level.
>     //
>     RigidRegistrationType::ShrinkFactorsArrayType vRigidShrinkFactors;
>     vRigidShrinkFactors.SetSize( 3 );
>     vRigidShrinkFactors[0] = 4;
>     vRigidShrinkFactors[1] = 2;
>     vRigidShrinkFactors[2] = 1;
>     vRigidRegistration->SetShrinkFactorsPerLevel( vRigidShrinkFactors );
>
>     //- Smoothing sigmas array
>     //
>     RigidRegistrationType::SmoothingSigmasArrayType vRigidSmoothingSigmas;
>     vRigidSmoothingSigmas.SetSize(3);
>     vRigidSmoothingSigmas.Fill(0);
>     vRigidRegistration->SetSmoothingSigmasPerLevel(vRigidSmoothingSigmas);
>
>     //- Observer
>     //
>     typedef CommandIterationUpdate< RigidRegistrationType > CommandType;
>     CommandType::Pointer observer = CommandType::New();
>     vRigidRegistration->AddObserver( itk::InitializeEvent(), observer );
>
>     try
>     {
>         std::cout << "Starting rigid registration..." << std::endl;
>         vRigidRegistration->Update();
>         std::cout << "Rigid parameters after registration: " << std::endl
>                 << vOptimizer->GetCurrentPosition() << std::endl;
>     }
>     catch( itk::ExceptionObject &e )
>     {
>         std::cerr << "Exception caught: " << e << std::endl;
>         return false;
>     }
>
> vComposite->AddTransform(const_cast<RigidTransformType*>(vRigidRegistration->GetOutput()->Get()));
>     return true;
> }
> //
> //=== End of the code =====================================
>
>
> _____________________________________
> Powered by www.kitware.com
>
> Visit other Kitware open-source projects at
> http://www.kitware.com/opensource/opensource.html
>
> Kitware offers ITK Training Courses, for more information visit:
> http://www.kitware.com/products/protraining.php
>
> Please keep messages on-topic and check the ITK FAQ at:
> http://www.itk.org/Wiki/ITK_FAQ
>
> Follow this link to subscribe/unsubscribe:
> http://www.itk.org/mailman/listinfo/insight-users
>