[Insight-users] Transforming a Volume
Luis Ibanez
luis.ibanez@kitware.com
Sat, 14 Sep 2002 23:14:46 -0400
HI Suresh,
Your registration code looks right for the most part.
Here are however some things that you could improve.
(and also answers to some of your questions)
1) The "guess" array is the initial set of parameters for the
transformation. Here a QuaternionRigidTransform is used.
This transform uses a UnitQuaternion and a Vector to represent
a rigid transform. The UnitQuaternion (Whose correct name is
"Versor", as Hamilton named in the 1800's) represents a rotation
in 3D while the vector represents the translation. The Quaternion
is defined by four numbers and the vector is defined by three number.
This makes a total of 7 parameters. The "guess" array is the sequence
of these 7 parameters. First the four coefficients of the Quaternion,
then, the 3 components of the vector. Note that if you want to play
with the values the quaternion coefficients you have to make sure
that their squared values sum 1.0.
2) I couldn't find the definition of the FixedImageRegion in your code.
The registration is performed over a user defined region.
You especify this region with the method: SetFixedImageRegion();
A typical setup could be:
registration->SetFixedImageRegion( fixed->GetBufferedRegion() );
3) The "learning rate" is a very critical parameter in this registration
method. It defines how big the step of the optimizer are. This value
is used to multiply the gradient of the metric on the parametric
space and compute a vector that moves the current set of parameters
to the new set of transform parameters to be tested in the next
iteration of the optimization. You will need some time for
fine-tunning this parameter to your data. You may want ot start
with small values (which will result in smaller step but more stable
behavior) and progresively try larger values.
4) Once the registration is doen, you can map volumes using the
Resample Image filter. However instead of ussing an affine
transform it will probably be better for you to use the actual
QuaternionRigidTransform that resulted from the registration
process. The ResampleImageFilter expect the transform class
to be polymorphic and just overload the Transform(point) method.
Please Keep in mind that the transform resulting from the
registration process is the one that map the Fixed image into
the Moving image. In your code it seems that the Fixe = MRI
and the moving=SPECT. That means that the final transform
is not the ones that "expands" SPECT to the MRI size but rather
the one that "shrinks" MRI to the SPECT size.
As it doesn't makes much sense to shrink the MRI to the resolution
of SPECT, you may want to compute the inverse of the Quaternion
Rigid Transform, and plug this inverse into the Resample image
transform.
Oopss... BTW, the RIGID transform will no change the scale of
the image (Rigid is only rotation and translation). If you also
need scale, you may have to switch to the "SimilarityTranform".
5) You may also want to change the optimizer and use the
RegularStepGradientDescentOptimizer which doesn't have the
"learning rate" parameter and is more stable in performing
optimization.
Please let us know if you have any other questions,
Thanks
Luis
===========================================
suresh wrote:
> Hi All,
>
> Tahnsks Luis.Thanks Sayan.. With your help i could complete the
> segmentation using NeighborhoodConnectedImageFilter.
>
> Now i'm facing aproblem in tarnsforming. My problem is ..
>
> I have an MRI and SPECT scan of different sizes. I'm registering SPECT
> to MRI using the code listed at the end.
> Then i'm usingthe output matrix to tarnsform the SPECT . What i'm
> excpecting is the SPECT to GROW to the MRI volume size.
> I tried with ResampleImageFilter but i could not get it to work.
>
> Please help me with the following issues..
> 1. Is my Registration code correct.?Or do i need change any parameters.
> 2. How do i apply a transformation matrix to a vaolume.?
>
> ____________________________Registration Code
> ___________________________________________________
> //******This code is written based on itk MIRegistartion example.***/
>
> typedef itk::QuaternionRigidTransform<double> TransformType;
> typedef itk::QuaternionRigidTransformGradientDescentOptimizer
> 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; // Image source
> UCharImage::Pointer fixedImage = converter.GetImage(MRI);//
> converter is my Image Source
> UCharImage::Pointer movingImage = converter.GetImage(SPECT);
>
> RegistrationType1::ParametersType
> guess(transform->GetNumberOfParameters() );
>
> //******What is this guess for????//
> guess[0] = 0.0; guess[1] = 0.0; guess[2] = 0.0; guess[3] =
> 1.0;
> guess[4] = 20.0; guess[5] = 40.0; guess[6] = 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 < 7; j++ )
> {
> scales[j] = 1.0 / vnl_math_sqr(300.0);
> }
> optimizer->SetNumberOfIterations( 7 );
>
> //******What is this learning rate??
> optimizer->SetLearningRate( 0.0000001 );
>
> registration->SetMetric(metric);
> registration->SetOptimizer(optimizer);
> registration->SetTransform(transform);
> registration->SetInterpolator(interpolator);
> registration->SetFixedImage(fixedImage);
> registration->SetMovingImage(movingImage);
> optimizer->SetScales(scales);
> registration->StartRegistration();
> ProgressUpdate(60, "Performing Normalized Registration15",
> MRIVolumeName);
> 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;
> _______________________________________________________________________________________________
>
> ***************Transform Code*********************
> typedef itk::ResampleImageFilter<UCharImage, UCharImage>
> ResampleFilter;
> typedef ResampleFilter::TransformType TransformType ;
> typedef UCharImage::SizeType ImageSizeType;
> typedef TransformType::MatrixType MatrixType;
> typedef TransformType::OffsetType VectorType;
> ConverterType converter ;
> ConverterType ::ImagePointer image = converter.GetImage(pVol);
> ImageSizeType size;
> size[0]=pTemp->width, size[1]= pTemp->height, size[2] = pTemp->depth;
>
> ResampleFilter::Pointer resampleFilter = ResampleFilter::New();
>
> TransformType::Pointer transform = resampleFilter->GetTransform();
> VectorType vector2 = transform->GetOffset();
> MatrixType matrix2 = transform->GetMatrix();
>
> // data is my Registration output matrix of size 4X4
> matrix2[0][0] = data[0][0]; matrix2[0][1] = matrix2[0][1] ;
> matrix2[0][2] = matrix2[0][2] ;
> matrix2[1][0] = data[1][0]; matrix2[1][1] = matrix2[1][1] ;
> matrix2[1][2] = matrix2[1][2] ;
> matrix2[2][0] = data[2][0]; matrix2[2][1] = matrix2[2][1] ;
> matrix2[2][2] = matrix2[2][2] ;
> transform->SetOffset(vector2);
> transform->SetMatrix(matrix2);
> UCharImage::Pointer outImage = NULL;
> resampleFilter ->SetTransform(transform);
> resampleFilter->SetInput(image);
> try{
> resampleFilter->Update();
> }catch(itk::ExceptionObject &Eo){
> Eo.GetDescription();
> }
> outImage = resampleFilter->GetOutput();
> _______________________________________________________________________________________________
>
>
> thanks,
>
> suresh
>
>
>
>
>
> _______________________________________________
> Insight-users mailing list
> Insight-users@public.kitware.com
> http://public.kitware.com/mailman/listinfo/insight-users
>