[Insight-users] Much more computation with Itk v4 as compared to v3 ?
lien lee
lienlee at gmail.com
Fri Oct 26 17:17:30 EDT 2012
Thanks for all you guys inputs.
The v3 and v4 code base are both from ITK 4.2.0. The two version code I
created both used 8 threads.
I do found that the v3 version code used smaller number of sample points.
So I changed the v3 code (as attached at the bottom) to use the same amount
of points as those of the v4 counterpart (ie., 10% of total amount of
voxels) at different resolution levels. I did this through the observer
with something like:
if (vRegistration->GetCurrentLevel() == 0)
{
Metric->SetNumberOfSpatialSamples(18176);
} else if (vRegistration->GetCurrentLevel() == 1)
{
Metric->SetNumberOfSpatialSamples( 151815 );
} else
{
Metric->SetNumberOfSpatialSamples( 1223322 );
}
To my eyes, everything looks similar (but, like I said, I am just a newbie
and not so sure about my judgement). I ran this changed v3 code, and it
took about 45 s. It is about half of the v4 code, which is around 90 s as
reported yesterday.
As for plotting the convergence curve, I am not so familiar of how to use
it (probably I would go over the ItkSoftwareGuide again). Seems to me,
multiple curves should be plotted with the 6 parameter transform. But,
anyway, I am going to try some plots.
BTW, the matching performance of the v4 code is a little bit better than
that of the v3 code (ie., DiceValue(v4-code) = 0.821 vs. DiceValue(v3-code)
= 0.803). But, I am not sure, since this is just on one data. However, I am
hoping that the automatically estimated learning rate and parameter scaling
in v4 would help.
//=== Start of v3 code ============================================
//
void RigidTransform(AffineType::Pointer vAffine, ImageType const&
vFixImage,
ImageType const& vMovImage)
{
//- Rotation center
typedef itk::Euler3DTransform<double> RigidType;
RigidType::Pointer vRigid = RigidType::New();
typedef itk::CenteredTransformInitializer<RigidType,
ImageType,
ImageType >
InitializerType;
InitializerType::Pointer Initializer = InitializerType::New();
Initializer->SetTransform( vRigid );
Initializer->SetFixedImage( &vFixImage );
Initializer->SetMovingImage( &vMovImage );
Initializer->GeometryOn();
Initializer->InitializeTransform();
//- Rigid registration
// ...
typedef itk::LinearInterpolateImageFunction<
ImageType,
double >
InterpolatorType;
typedef itk::GradientDescentOptimizer RigidOptimizerType;
//- Instantiate major objects.
RigidOptimizerType::Pointer Optimizer = RigidOptimizerType::New();
InterpolatorType::Pointer Interpolator = InterpolatorType::New();
RegistrationType::Pointer Registration = RegistrationType::New();
Registration->SetNumberOfLevels(3);
//- Assemble to the registration pipeline
Registration->SetOptimizer( Optimizer );
Registration->SetInterpolator( Interpolator );
Registration->SetFixedImage( &vFixImage );
Registration->SetMovingImage( &vMovImage );
ImageType::RegionType FixRegion = vFixImage.GetBufferedRegion();
Registration->SetFixedImageRegion( FixRegion );
Registration->SetInitialTransformParameters( vRigid->GetParameters() );
Registration->SetTransform( vRigid );
//- Optimizer
typedef RigidOptimizerType::ScalesType OptimizerScalesType;
OptimizerScalesType vOptimizerScales( vRigid->GetNumberOfParameters() );
const double vTranslationScale = 1.0 / 1000;
vOptimizerScales[0] = 1.0;
vOptimizerScales[1] = 1.0;
vOptimizerScales[2] = 1.0;
vOptimizerScales[3] = vTranslationScale;
vOptimizerScales[4] = vTranslationScale;
vOptimizerScales[5] = vTranslationScale;
Optimizer->SetScales( vOptimizerScales );
Optimizer->SetNumberOfIterations(50);
Optimizer->SetLearningRate(0.3);
RigidOptimizerType::ParametersType vParameters =
vRigid->GetParameters();
Optimizer->SetInitialPosition(vParameters);
//- Metric
//
typedef itk::MattesMutualInformationImageToImageMetric<
ImageType,
ImageType >
MetricType;
MetricType::Pointer Metric = MetricType::New();
vRegistration->SetMetric(Metric);
Metric->SetNumberOfHistogramBins(32);
Metric->ReinitializeSeed( 76926294 );
//- Setup observers
//
RigidRegistrationInterfaceCommand::Pointer vInterface =
RigidRegistrationInterfaceCommand::New();
Registration->AddObserver(itk::IterationEvent(), vInterface);
OptimizerObserver::Pointer vOptimizerObserver =
OptimizerObserver::New();
vOptimizerObserver->mRigidMetric =
dynamic_cast<itk::MattesMutualInformationImageToImageMetric<ImageType,
ImageType>*>
(Registration->GetMetric());
Optimizer->AddObserver(itk::StartEvent(), vOptimizerObserver);
Optimizer->AddObserver(itk::IterationEvent(), vOptimizerObserver);
//- Do Translation Registration
try
{
Registration->Update();
}
catch( itk::ExceptionObject & err )
{
std::cerr << "ExceptionObject caught !" << std::endl;
std::cerr << err << std::endl;
}
std::cout << "Rigid Registration completed" << std::endl;
std::cout << std::endl;
vRigid->SetParameters(Registration->GetLastTransformParameters());
}
//
//=== End of v3 code =============================================
2012/10/26 brian avants <stnava at gmail.com>
> 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
> >> >
> >
> >
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20121026/c8e3c0eb/attachment.htm>
More information about the Insight-users
mailing list