[Insight-users] Question regarding centered versor transform

[Dan Mueller]

Hi Insight Users,

I am experiencing a problem with setting parameters for centered
versor transforms.

Some background: I ultimately want to perform registration considering
only (centered) rotation. From my understanding the VersorTransform
compactly represents a centered rotation in 3D space (three rotation
parameters, no translation parameters optimized). However, when the
VersorTransformOptimizer optimizes the metric I seem to be noticing
"translation". I think this "translation" is occurring because the
center of rotation/translation/offset ? is not being set correctly
when setting the versor parameters.

To further investigate the issue I have simplified the problem,
totally removing the registration component, only considering image
resampling with a VersorTransform. The source code attached below does
the following:
1. Read the fixed image
2. Read the moving image
3. Create and setup the transform
4. Create interpolator
5. Resample moving image using the transform
6. Write resampled image

What I am observing is the following:
a. If I set the rotation via Versor.SetRotation(.) everything behaves
as I expect (i.e. the image is rotated around the center).
b. If I set the rotation via Transform.SetParameters(.) (which is what
the optimizer does) the rotation is NOT around the (desired) center.

I have made screen shots showing the results of (a) and (b) above,
only changing the #define SETBYPARAMETERS.
    https://fileshare.qut.edu.au/public/muellerd/resampled_direct.png
    https://fileshare.qut.edu.au/public/muellerd/resampled_param.png

The problem described in (b) seems to be what I am noticing with my
registration (i.e. a rotation around the wrong center, introducing
"translation").

What am I doing wrong? Have I misunderstood something? Is this a bug
with VersorTransform?

Thanks for any advice!

Regards, Dan

=== CMake ===
PROJECT(VersorTest)
CMAKE_MINIMUM_REQUIRED(VERSION 2.6)

FIND_PACKAGE(ITK REQUIRED)
INCLUDE(${ITK_USE_FILE})

ADD_EXECUTABLE(VersorTest main.cxx)
TARGET_LINK_LIBRARIES(VersorTest ITKCommon ITKIO)

=== Source Code (main.cxx) ===
#include "itkImage.h"
#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkVersorTransform.h"
#include "itkLinearInterpolateImageFunction.h"
#include "itkResampleImageFilter.h"

#define SETBYPARAMETERS 0

int main(int argc, char* argv[]) {

    try {
        // Typedefs
        const unsigned int Dimension = 3;
        typedef float PixelType;
        typedef double CoordRepType;
        typedef itk::Image< PixelType, Dimension > ImageType;
        typedef ImageType::PointType PointType;
        typedef itk::Vector< CoordRepType, Dimension > VectorType;

        // 1. Read fixed image
        std::cout << "Reading fixed image..." << std::endl;
        typedef itk::ImageFileReader< ImageType > ReaderType;
        ReaderType::Pointer readerFixed = ReaderType::New( );
        readerFixed->SetFileName( "D:/Temp/fixed.mhd" );
        readerFixed->Update( );
        ImageType::Pointer imageFixed = readerFixed->GetOutput();
        imageFixed->DisconnectPipeline( );

        // 2. Read moving image
        std::cout << "Reading moving image..." << std::endl;
        ReaderType::Pointer readerMoving = ReaderType::New( );
        readerMoving->SetFileName( "D:/Temp/moving.mhd" );
        readerMoving->Update( );
        ImageType::Pointer imageMoving = readerMoving->GetOutput();
        imageMoving->DisconnectPipeline( );

        // 3. Create and setup transform
        // NOTE: This basically does the same as CenteredTransformInitializer
        // with geometry on
        std::cout << "Creating transform..." << std::endl;
        typedef itk::VersorTransform< CoordRepType > TransformType;
        TransformType::Pointer transform = TransformType::New();
        PointType centerFixed;
        for (int i=0; i<Dimension; i++) {
            centerFixed[i] =
                imageFixed->GetOrigin()[i] +
                imageFixed->GetLargestPossibleRegion().GetSize()[0] *
                imageFixed->GetSpacing()[0] / 2.0;
        }
        PointType centerMoving;
        for (int i=0; i<Dimension; i++) {
            centerMoving[i] =
                imageMoving->GetOrigin()[i] +
                imageMoving->GetLargestPossibleRegion().GetSize()[0] *
                imageMoving->GetSpacing()[0] / 2.0;
        }
        VectorType translation = centerFixed - centerMoving;
        TransformType::VersorType versor;
        versor.Set( 0.0, 0.0, 0.9, 0.43589 );
        transform->SetCenter( centerFixed );
        transform->SetTranslation( translation );
        TransformType::ParametersType param( Dimension );
        for (int i=0; i<Dimension; i++) {
            param[i] = versor.GetRight()[i];
        }
#if SETBYPARAMETERS
        std::cout << "Setting rotation using parameters..." << param
<< std::endl;
        transform->SetParameters( param );
#else
        std::cout << "Setting rotation directly..." << versor << std::endl;
        transform->SetRotation( versor );
#endif

        // 4. Create interpolator
        typedef itk::LinearInterpolateImageFunction< ImageType,
CoordRepType > InterpolatorType;
        InterpolatorType::Pointer interpolator = InterpolatorType::New();

        // 5. Resample
        std::cout << "Resampling moving image..." << std::endl;
        typedef itk::ResampleImageFilter< ImageType, ImageType> ResampleType;
        ResampleType::Pointer resample = ResampleType::New();
        resample->SetInput( imageMoving );
        resample->SetOutputParametersFromImage( imageFixed );
        resample->SetInterpolator( interpolator );
        resample->SetTransform( transform );
        resample->Update( );

        // 6. Write resample image
        std::cout << "Writing resampled image..." << std::endl;
        typedef itk::ImageFileWriter< ImageType > WriterType;
        WriterType::Pointer writer = WriterType::New( );
        writer->SetFileName( "D:/Temp/resampled.mhd" );
        writer->SetInput( resample->GetOutput() );
        writer->Update( );

        return EXIT_SUCCESS;
    }
    catch( itk::ExceptionObject& err ) {
        std::cout << err << std::endl;
        return EXIT_FAILURE;
    }
}

=== Information ===
ITK 3.8
CMake 2.6.0
Windows XP 64-bit SP2
Visual Studio 8.0.50727.762