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

[lien lee]

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 =====================================
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20121025/a7ce5f69/attachment.htm>