[Insight-users] Re: AffineTransform
Luis Ibanez
luis.ibanez@kitware.com
Mon, 16 Sep 2002 22:33:10 -0400
Hi Suresh,
1) The Affine transform in 3D is composed by a 3x3 matrix
and a 3 dimensional vector.
The matrix represents rotations, scaling and shearing
while the vector represents translations.
The 12 parameters returned by the affine transform in 3D
are the 9 (=3x3) matrix coefficients plus the 3 vector
components. They are ordered as:
0 1 2
Matrix = 3 4 5
6 7 8
9
Vector = 10
11
2) The piece of code:
> vnl_quaternion<double>
> quat(solution[0],solution[1],solution[2],solution[3]);
>
> vnl_matrix_fixed<double,3,3> vnlmat = quat.rotation_matrix();
>
Is only valid if "solution" is the set of parameters corresponding
to a QuaternionRigidTransform. In this case, the array of parameters
has 7 values corresponding to 4 quaternion coefficients followed by
3 vector components.
If you are replacing the QuaternionRigidTransform by an
AffineTransform, this piece of code is no longer valid.
However, the good news is that the Affine transform will provide
directly the matrix to you. Access to the rotation matrix and the
translation vector is provided by the methods:
affineTransform.GetMatrix();
affineTransform.GetOffset();
3) You may want to keep the MRI as the Fixed image and the SPECT as the
moving image. Once the registration is done, the final transform
indicates how to map points from the Fixed image space to the Moving
image space. In order to map one point of the SPECT image on the
space of the MRI image you need the inverse of the transform.
Note however that in order to map the SPECT image on top of the MRI
image you use directly the transform resulting from the registration.
This is due to the way in which mapping is done. In practice you
visit every pixel of the MRI image and compute what is the SPECT
pixel corresponding to this location. This actually requires to map
a point of the MRI image into the SPECT space.
This means, that in order to reformat the SPECT image and see how
it will look after being deformed to fit the MRI image, you can use
the ResampleImageFilter by giving to it the actual transform
resulting from the registration procedure. You don't need the
inverse of the transform for this.
4) BTW A number of 7 iterations is probably insufficient for
registration, typically 300 to 500 maybe a better choice.
Please let us know if you find any problems,
Thanks
Luis
=================================================================
suresh wrote:
> Hi Luis,
>
> Thanks a lot for help.It will be alot more helpful if you clarify a few
> things on registration with AffineTransform
>
> 1. What about the transform parametres. I observed(by calling
> transform->GetNumberOfParameters()) that there are 12 to be given.What
> are they?
> 2. Is this peice of code still valid?
> vnl_quaternion<double>
> quat(solution[0],solution[1],solution[2],solution[3]);
> vnl_matrix_fixed<double,3,3> vnlmat = quat.rotation_matrix();
> 3. I'm not going to swap the fixed and moving images(MRI, SPECT), but
> while Transformation i'll get the Inverse matrix and apply.Is that right?
>
> This is modified code. for registration.
>
> Matrix* mat = NULL;
> double result[16];
>
> typedef BufferToImageConversion<unsigned char> ConverterType;
> typedef ConverterType::ImageType UCharImage;
>
> typedef itk::AffineTransform<double> TransformType;
> typedef itk::RegularStepGradientDescentBaseOptimizer OptimizerType;
> typedef itk::LinearInterpolateImageFunction <UCharImage, double>
> InterpolatorType;
> typedef itk::ImageRegistrationMethod <UCharImage , UCharImage >
> RegistrationType1;
>
> typedef itk::MutualInformationImageToImageMetric<UCharImage ,
> UCharImage> MetricType;
> ProgressUpdate(50, "Performing MutualInformation Registration",
> MRIVolumeName);
> MetricType::Pointer metric = MetricType::New();
> TransformType::Pointer transform = TransformType::New();
> OptimizerType::Pointer optimizer = OptimizerType::New();
> InterpolatorType::Pointer interpolator = InterpolatorType::New();
> RegistrationType1::Pointer registration =
> RegistrationType1::New();
>
> ConverterType converter;
> UCharImage::Pointer fixedImage = converter.GetImage(MRI);
> UCharImage::Pointer movingImage = converter.GetImage(SPECT);
>
> RegistrationType1::ParametersType
> guess(transform->GetNumberOfParameters() );
>
> CString str;
> str.Format(" Count %d", transform->GetNumberOfParameters());
> AfxMessageBox(str);
>
> guess[0] = 0.0; guess[1] = 0.0; guess[2] = 0.0; guess[3] =
> 1.0;
> guess[4] = 0.0; guess[5] = 0.0; guess[6] = 0.0; guess[7] = 0.0;
> guess[8] = 0.0; guess[9] = 40.0; guess[10] = 20.0; guess[11] =
> 0.0;
> try{
> registration->SetInitialTransformParameters (guess);
> metric->SetMovingImageStandardDeviation( 2 );
> metric->SetFixedImageStandardDeviation( 2 );
> metric->SetNumberOfSpatialSamples( 1000 );
> typedef OptimizerType::ScalesType ScaleType;
> ScaleType scales(transform->GetNumberOfParameters());
> scales.Fill( 1.0 );
> for( unsigned j = 4; j < 12; j++ )
> {
> scales[j] = 1.0 / vnl_math_sqr(300.0);
> }
> //Scale parameters
> optimizer->SetNumberOfIterations( 7 );
>
> //set registration parameters
> registration->SetMetric(metric);
> registration->SetOptimizer(optimizer);
> registration->SetTransform(transform);
>
> registration->SetInterpolator(interpolator);
>
> registration->SetFixedImage(movingImage);
> registration->SetMovingImage(fixedImage);
> registration->SetFixedImageRegion(
> fixedImage->GetBufferedRegion() );
> // Setup the optimizer
> optimizer->SetScales(scales);
> registration->DebugOn();
> registration->StartRegistration();
>
> ProgressUpdate(60, "Performing Normalized Registration15",
> MRIVolumeName);
> // Get the results
> itk::Array<double> *arr=new itk::Array<double>;
> RegistrationType1::ParametersType solution
> =registration->GetLastTransformParameters();
>
> vnl_quaternion<double>
> quat(solution[0],solution[1],solution[2],solution[3]);
>
> vnl_matrix_fixed<double,3,3> vnlmat = quat.rotation_matrix();
>
> result[0] = vnlmat[0][0]; result[1] = vnlmat[0][1];
> result[2] = vnlmat[0][2]; result[3] = solution[4];
> result[4] = vnlmat[1][0]; result[5] = vnlmat[1][1];
> result[6] = vnlmat[1][2]; result[7] = solution[5];
> result[8] = vnlmat[2][0]; result[9] = vnlmat[2][1];
> result[10] = vnlmat[2][2]; result[11] = solution[6];
> result[12] = 0; result[13] = 0; result[14] = 0;
> result[15] = 1;
>
> With this code i could not go beyond the StartRegistration call? Please
> help in fixing this..
>
> And in Resample code i'm using the Inverse of Transform.
> TransformType::Pointer inverseTransform = transform->Inverse();
> resampleFilter ->SetTransform(inverseTransform);
>
> Thanks
> suresh
>
>