[Insight-users] Question regarding centered versor transform

[Dan Mueller]

Hi Simon, thanks for your reply.

No, I did not called SetFixedParameters(.). Doing so makes everything
behave as I expect. A bug with the user :P

Thanks again.

Regards, Dan

2008/10/30 Simon Warfield <simon.warfield at childrens.harvard.edu>:
>
> The different types of transforms have inconsistent APIs and it is tricky to
> know exactly what should be set.
>
> There is an interaction between the 'fixed parameters' and the 'offset'.
>
> Did you set the 'Fixed Parameters' before you set the 'Parameters' ?
> http://www.itk.org/Doxygen/html/classitk_1_1MatrixOffsetTransformBase.html#bab57111fda06e69309c80fdf6fbee38
>
>
>> 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.
>
>
> --
> Simon
>
> Message: 6
> Date: Thu, 30 Oct 2008 10:44:45 +0100
> From: "Dan Mueller" <dan.muel at gmail.com>
> Subject: Re: [Insight-users] Question regarding centered versor
>        transform
> To: "Insight Users" <insight-users at itk.org>
> Message-ID:
>        <2b8d077d0810300244t5dc7ffc4w7972969dfc613fbc at mail.gmail.com>
> Content-Type: text/plain; charset=ISO-8859-1
>
> 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