[Insight-users] Question regarding centered versor transform

[Dan Mueller]

I've had a dig into the code and seem to have found an inconsistency...

In itkVersorTransform.txx in SetRotation(.) the versor is set and then
both ComputeMatrix(.) and ComputeOffset(.) are called.

In itkVersorTransform.txx in SetParameters(.) the verser is set and
then only ComputeMatrix(.) is called. If I add ComputeOffset(.) here,
then my code works as expected.

Should I file a bug, or have I misunderstood something?

Regards, Dan

2008/10/28 Dan Mueller <dan.muel at gmail.com>:
> 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
>