[Insight-users] Much more computation with Itk v4 as compared to v3 ?
brian avants
stnava at gmail.com
Fri Oct 26 09:55:51 EDT 2012
just to be a bit clearer ---- there are several key points in a
reasonable comparison:
1. how many threads are you using?
2. how many points are you using to sample the mutual information?
3. what is the shape of the convergence curve showing the mutual information?
4. what version of v4 are you using?
point 2 is probably the most critical.
for instance, the v4 registration methods have the option
registration->SetMetricSamplingPercentage(0.1);
which ( with the above setting ) will sample 10% of your image -----
at each resolution.
this is very different than v3 which , in general , uses a fixed
number of samples.
v4 supports that strategy too but it is not the default.
so, in comparison, you may have set up your v3 metric with something like
metric->SetNumberOfSamples( 20000 )
which would only be about 0.1 % of the image points.
so you might first try doing
registration->SetMetricSamplingPercentage(0.01);
to see how this affects speed.
brian
On Fri, Oct 26, 2012 at 9:12 AM, lien lee <lienlee at gmail.com> wrote:
> Thanks for the informative reply, Brian. I am going to double check, and
> will let you know what I could get.
>
>
>
> 2012/10/25 brian avants <stnava at gmail.com>
>>
>> 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
>> >
>
>
More information about the Insight-users
mailing list