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

[lien lee]

Yes, they are.


2012/10/26 Bill Lorensen <bill.lorensen at gmail.com>

> 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/5a5eca8e/attachment.htm>