[Insight-users] Question regarding centered versor transform

[Dan Mueller]

Is someone able to confirm if this is a bug with the code or with the user?

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