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

[Bill Lorensen]

Are both versions built Release?

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
>
>


-- 
Unpaid intern in BillsBasement at noware dot com
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20121026/9872acbb/attachment.htm>