[Insight-users] 3d multi-modality registration
Atwood, Robert C
r.atwood at imperial.ac.uk
Wed May 31 15:01:16 EDT 2006
Hi Grace:
Here is what I have done to gain a useful interpretation of the
transform matrix -- use the features already provided by the underlying
numeric library to get some decomposition of the matrix. I was not very
familiar with this aspect of linear algebra a couple of weeks ago! I
think I may have learned it in 2nd year from Prof. Greub, who was
apparently an expert in the field
http://www.amazon.com/gp/product/0387901108/102-0261825-5847335?v=glance
&n=283155
but that was 1982. (in Toronto)
So perhaps there is a better way and I would be glad to hear about it.
Anyways I have found the Singular Value Decomposition most useful in
interpreting what the transform is 'like' , the W and V seem to
represent a symmetric shearing and the U the additional rotation. I have
also read that the Polar Decomposition is a useful interpretation as
well, also that of the 3 eigenvalues of the original matrix, one will be
real, and the corresponding eigenvector will be the axis of rotation.
Here I had some code for getting all these, later I decided that the
S.V.D. was most useful for my purpose, which was to find the direction
of the maximum stretch due to the affine registration.
--Robert
#include <vnl/algo/vnl_svd.h>
//object that stores a matrix
typedef vnl_matrix<double> MyMatrix;
//object that creates the singular value decomposition
typedef vnl_svd<double> SVDGetter;
// object that creates the eigenvalue decomposition (complex
numbers)
typedef vnl_real_eigensystem EigenGetter;
... Do some registration .... Where finalTransform contains the final
transform of the registration...
... I was using ScalableAffineTransform ...
// provide a container for the underlying matrix
MyMatrix matrix; //the transformation matrix
MyMatrix myU, myW, myV; //SVD decomposition
MyMatrix myPQ, myPP; //Polar decomposition
EigenGetter * eigen; //Eigenvalue decompostion (complex)
matrix = finalTransform->GetMatrix().GetVnlMatrix();
// create the eigenvalue decomposition
eigen = new EigenGetter(matrix);
std::cout << "matrix: "<<std::endl << matrix << std::endl;
std::cout << "eigenvectors: "<<std::endl <<std::endl <<
eigen->V << std::endl;
std::cout << "eigenvalues: " << eigen->D << std::endl;
// create the singular value decomposition
SVDGetter mysvd( matrix);
myU=mysvd.U();
myV=mysvd.V();
myW=mysvd.W();
std::cout <<"U:" <<std::endl<< myU << std::endl;
std::cout <<"V:" <<std::endl<< myV << std::endl;
std::cout <<"W:" <<std::endl<< myW << std::endl;
#ifdef POLARDECOMP
// compute the polar decompostion from the SVD
myPP = myV * myW * myV.transpose();
myPQ = myU * myV.transpose();
std::cout <<std::endl<<"PP:" << std::endl << myPP <<
std::endl;
std::cout << std::endl <<"PQ:" << std::endl<< myPQ <<
std::endl;
#endif /* POLARDECOMP */
-----Original Message-----
From: insight-users-bounces+r.atwood=imperial.ac.uk at itk.org
[mailto:insight-users-bounces+r.atwood=imperial.ac.uk at itk.org] On Behalf
Of Grace Chen
Sent: 31 May 2006 18:11
To: Karthik Krishnan
Cc: insight-users at public.kitware.com
Subject: [Insight-users] 3d multi-modality registration
Hi Karthik,
I just tried that...it gives me the same result as before! To make sure
I
understand what you said correctly, I include my code below.
Thanx a lot!!
Grace
--------------------------------------------------------------------
void vtkItkMultiCenteredImReg::Execute()
{
movingImportCaster->SetInput( movingItkImporter->GetOutput() );
fixedImportCaster->SetInput( fixedItkImporter->GetOutput() );
//movingImage = movingImportCaster->GetOutput();
fixedImage = fixedImportCaster->GetOutput();
fixedImage->Update();
//movingImage->Update();
transform->SetIdentity();
resample->SetTransform( transform );
resample->SetInput( movingImportCaster->GetOutput() );
resample->SetSize( fixedImage->GetLargestPossibleRegion().GetSize() );
resample->SetOutputOrigin( fixedImage->GetOrigin() );
resample->SetOutputSpacing( fixedImage->GetSpacing() );
movingImage = resample->GetOutput();
movingImage->Update();
registration->SetMetric( metric );
registration->SetOptimizer( optimizer );
registration->SetInterpolator( interpolator );
registration->SetTransform( transform );
registration->SetFixedImage(fixedImage);
registration->SetMovingImage(movingImage);
registration->SetFixedImageRegion(
fixedImage->GetBufferedRegion() );
const FixedImageType::SpacingType&
fixedSpacing = fixedImage->GetSpacing();
const FixedImageType::PointType&
fixedOrigin = fixedImage->GetOrigin();
FixedImageType::SizeType fixedSize =
fixedImage->GetLargestPossibleRegion().GetSize();
TransformType::InputPointType centerFixed;
centerFixed[0] = fixedOrigin[0] + fixedSpacing[0] * fixedSize[0] / 2.0;
centerFixed[1] = fixedOrigin[1] + fixedSpacing[1] * fixedSize[1] / 2.0;
centerFixed[2] = fixedOrigin[2] + fixedSpacing[2] * fixedSize[2] / 2.0;
const MovingImageType::SpacingType&
movingSpacing = movingImage->GetSpacing();
const MovingImageType::PointType&
movingOrigin = movingImage->GetOrigin();
MovingImageType::SizeType movingSize =
movingImage->GetLargestPossibleRegion().GetSize();
TransformType::InputPointType centerMoving;
centerMoving[0] = movingOrigin[0] + movingSpacing[0] * movingSize[0] /
2.0;
centerMoving[1] = movingOrigin[1] + movingSpacing[1] * movingSize[1] /
2.0;
centerMoving[2] = movingOrigin[2] + movingSpacing[2] * movingSize[2] /
2.0;
transform->SetTranslation(centerMoving-centerFixed);
/*
initializer->SetTransform(transform);
initializer->SetFixedImage(fixedImage);
initializer->SetMovingImage(movingImage);
initializer->GeometryOn();
//initializer->MomentsOn();
initializer->InitializeTransform();
VersorType rotation;
VectorType axis;
axis[0] = 0.0;
axis[1] = 0.0;
axis[2] = 1.0;
const double angle = 0;
rotation.Set(axis, angle);
transform->SetRotation(rotation);
*/
registration->SetInitialTransformParameters( transform->GetParameters()
);
OptimizerScalesType optimizerScales( transform->GetNumberOfParameters()
);
const double translationScale = 1.0 / 1000.0;
optimizerScales[0] = 1.0;
optimizerScales[1] = 1.0;
optimizerScales[2] = 1.0;
optimizerScales[3] = translationScale;
optimizerScales[4] = translationScale;
optimizerScales[5] = translationScale;
optimizer->SetScales( optimizerScales );
//optimizer->MinimizeOn();
//optimizer->SetMaximumStepLength(0.2);
//optimizer->SetMinimumStepLength(0.0001);
optimizer->SetNumberOfIterations( 50 );
metric->SetNumberOfHistogramBins(20);
metric->SetNumberOfSpatialSamples(10000);
optimizer->AddObserver(itk::IterationEvent(), observer);
registration->AddObserver(itk::IterationEvent(), command);
registration->SetNumberOfLevels(3);
try
{
registration->StartRegistration();
}
catch( itk::ExceptionObject & err )
{
std::cout << "ExceptionObject caught !" << std::endl;
std::cout << err << std::endl;
return ;
}
finalParameters = registration->GetLastTransformParameters();
finalTransform->SetParameters( finalParameters );
return ;
}
------------------------------------------------------------------------
----
---------------
----- Original Message -----
From: "Karthik Krishnan" <Karthik.Krishnan at kitware.com>
To: "Grace Chen" <Grace.Chen at swri.ca>; "itk"
<insight-users at public.kitware.com>
Sent: Wednesday, May 31, 2006 12:27 PM
Subject: [Insight-users] Re: 3d multi-modality registration
> No. In the end, you should not need to do this . If you use the
> ResampleImageFilter to resample your moving image as in ,
>
> ResampleFilterType::Pointer resample = ResampleFilterType::New();
> resample->SetTransform( finalTransform ); // after registration
> resample->SetInput( movingImageReader->GetOutput() );
> resample->SetSize(
fixedImage->GetLargestPossibleRegion().GetSize() );
> resample->SetOutputOrigin( fixedImage->GetOrigin() );
> resample->SetOutputSpacing( fixedImage->GetSpacing() );
>
> Your resample filter will ensure that the moving image is resampled to
> have the same meta-data as the fixed image. Does this not work for you
?
>
>
> Try to take things one step at a time. Just resample the moving image
> first (without registration using an identity transform). Then use a
> translation transform that overlays them to a reasonable extent. Then
> register with this translation transfrom or a transform really close
to
> it as the initial transform and see if your solution converges to the
> right translation transform.
>
> After you've gone through this exercise, you can dive into the other
> details of registration.
>
> HTH
> karthik
>
>
>
> Grace Chen wrote:
>
> >Hi Karthik,
> >
> >Do you mean I should apply scalling (as in the following line) after
> >applying the registration transformation to make corresponding slices
of
> >both images matched??
> >
>
>currTransform->Scale((double)fSpacing[0]/mSpacing[0],(double)fSpacing[1
]/mS
p
> >acing[1], (double)fSpacing[2]/mSpacing[2]);
> >
> >Thank you so much!
> >
> >Grace
> >
> >
> >----- Original Message -----
> >From: "Karthik Krishnan" <Karthik.Krishnan at kitware.com>
> >To: "Grace Chen" <Grace.Chen at swri.ca>
> >Cc: <insight-users at itk.org>
> >Sent: Wednesday, May 31, 2006 10:29 AM
> >Subject: Re: 3d multi-modality registration
> >
> >
> >
> >
> >>Grace Chen wrote:
> >>
> >>
> >>
> >>>Hi there,
> >>>
> >>>My moving image looks smaller than the fixed image on screen.
(These
two
> >>>input images have different spacing.) And because it's 3D
multi-modality
> >>>registration, so only translation and rotation are involved in the
> >>>registration process. Is this why the registered moving image
still
look
> >>>smaller than the fixed image and the slices of the moving images
don't
> >>>
> >>>
> >match
> >
> >
> >>>that of the fixed image?
> >>>
> >>>
> >>>
> >>>
> >>This happens often in multi-modality registration, where CT, MR PET
> >>datasets have very different resolutions. If you use the
registration
> >>framework in ITK, the moving image is resampled to the grid of the
fixed
> >>image. In otherwords your resampled/transformed moving image should,
> >>after registration have the same meta-data (spacing, origin) as the
> >>fixed image. You definitely want to do this so you evaluate the
> >>registatration using an overlay or checkerboard.
> >>
> >>
> >>
> >>>Thanx a lot!
> >>>
> >>>Grace
> >>>
> >>>
> >>>----- Original Message -----
> >>>From: "Grace Chen" <Grace.Chen at swri.ca>
> >>>To: <Karthik.Krishnan at kitware.com>
> >>>Cc: <insight-users at itk.org>
> >>>Sent: Friday, May 26, 2006 6:00 PM
> >>>Subject: Re: [Insight-users] 3d multi-modality registration
> >>>
> >>>
> >>>
> >>>
> >>>
> >>>
> >>>>Hi Karthik,
> >>>>
> >>>>Thanx for your help!
> >>>>
> >>>>However, my program needs this process to be made automatic....
So, I
> >>>>extracted the volume of interest from the moving image, making
sure
that
> >>>>
> >>>>
> >>>>
> >>>>
> >>>the
> >>>
> >>>
> >>>
> >>>
> >>>>moving image covers the whole fixed image but yet is not too big.
The
> >>>>registration program can correct the movement in the x and y
direction.
> >>>>
> >>>>
> >>>>
> >>>>
> >>>And
> >>>
> >>>
> >>>
> >>>
> >>>>I can see the program translates the slices (in z direction) too.
> >>>>
> >>>>
> >>>>
> >>>>
> >>>However,
> >>>
> >>>
> >>>
> >>>
> >>>>because the moving image has smaller spacing, so within the same
volume
> >>>>
> >>>>
> >of
> >
> >
> >>>>interest, the moving image has more slices. The program doesn't
seem
to
> >>>>interpolate the moving image well so the subsequent slices of the
> >>>>
> >>>>
> >>>>
> >>>>
> >>>registered
> >>>
> >>>
> >>>
> >>>
> >>>>moving image matches the corresponding slices of the fixed image.
> >>>>
> >>>>Any idea why? or is there a bug in my understanding??
> >>>>
> >>>>Thanx again!!
> >>>>
> >>>>Grace
> >>>>
> >>>>
> >>>>
> >>>>----- Original Message -----
> >>>>From: "Karthik Krishnan" <Karthik.Krishnan at kitware.com>
> >>>>To: "Grace Chen" <Grace.Chen at swri.ca>
> >>>>Cc: <insight-users at itk.org>
> >>>>Sent: Friday, May 26, 2006 3:11 PM
> >>>>Subject: Re: [Insight-users] 3d multi-modality registration
> >>>>
> >>>>
> >>>>
> >>>>
> >>>>
> >>>>
> >>>>>Hi Grace,
> >>>>>
> >>>>>The correct use of image registration is to bring it close to
final
> >>>>>solution and expect image registration to take it from there. I
you
> >>>>>provide two volumes with vastly differnt extents and a poor
> >>>>>initialization, it wouldn't be surprising if registration failed
to
> >>>>>achieve results.
> >>>>>
> >>>>>In a clinincal workflow, I don't think any radiologist/clinician
would
> >>>>>perform registration using command line tools and proceed to the
next
> >>>>>step without validating the quality of the registration. This is
why
> >>>>>
> >>>>>
> >its
> >
> >
> >>>>>usually done using GUI tools. For instance
> >>>>>InsightApplications/LandmarkInitializedMutualInformation (have
you
> >>>>>
> >>>>>
> >tried
> >
> >
> >>>>>this application for your datasets ?) allows you to place
landmarks
on
> >>>>>the source and target image, so you can roughly overlay the
images
> >>>>>(usually within 0-5 mm of each other) and then allow registration
to
> >>>>>fine tune it. If the overlay looks good, you are happy.
> >>>>>
> >>>>>Please give this application a try first :
> >>>>>
> >>>>>1. Specify a landmark on the fixed and moving image. Pick
anatomical
> >>>>>correspondances. For an MRI of the brain, you could pick the tips
of
> >>>>>
> >>>>>
> >the
> >
> >
> >>>>>splenium of the corpuscallosum on a sagittal/coronal acquistion
(you
> >>>>>
> >>>>>
> >can
> >
> >
> >>>>>see it very clearly in
> >>>>>
> >>>>>
> >Insight/Examples/Data/BrainMidSagittalSlice.png).
> >
> >
> >>>>>For an axial MRI, you could pick the tips of the ventricles.
> >>>>>
> >>>>>2. Initialize using landmarks. (it should be one of the options
on
the
> >>>>>initialization menu).
> >>>>>
> >>>>>3. Register
> >>>>>
> >>>>>Does this work for you ?
> >>>>>
> >>>>>HTH
> >>>>>karthik
> >>>>>
> >>>>>On Fri, 2006-05-26 at 12:26 -0400, Grace Chen wrote:
> >>>>>
> >>>>>
> >>>>>
> >>>>>
> >>>>>>Hi Luis, Thanx a lot for your inputs!
> >>>>>>
> >>>>>>I've been struggling with this one for a longest time.... The
problem
> >>>>>>
> >>>>>>
> >>>>>>
> >>>>>>
> >>>>is
> >>>>
> >>>>
> >>>>
> >>>>
> >>>>>>that I have two 3D brain images and one brain image has greater
> >>>>>>
> >>>>>>
> >>>>>>
> >>>>>>
> >>>physical
> >>>
> >>>
> >>>
> >>>
> >>>>>>extents in z. The information of the two volumes are as
follows:
> >>>>>>
> >>>>>>Fixed image:
> >>>>>> origin = [-120, -135.805, -29.7683]
> >>>>>> spacing = [0.9375, 0.9375, 7]
> >>>>>> extent = [256,256, 7]
> >>>>>>
> >>>>>>Moving image:
> >>>>>> origin = [-142.5, -170.488, -84.8781]
> >>>>>> spacing = [1.17188, 1.17188, 5.5]
> >>>>>> extent = [256,256, 28]
> >>>>>>
> >>>>>>For these two brain volumes, the fixed image matches a section
in
the
> >>>>>>
> >>>>>>
> >>>>>>
> >>>>>>
> >>>>moving
> >>>>
> >>>>
> >>>>
> >>>>
> >>>>>>image...I tried registered them using the whole volume, the
first
> >>>>>>
> >>>>>>
> >>>>>>
> >>>>>>
> >>>slice
> >>>
> >>>
> >>>
> >>>
> >>>>of
> >>>>
> >>>>
> >>>>
> >>>>
> >>>>>>the registered moving image doesn't looked like that of the
fixed
> >>>>>>
> >>>>>>
> >>>>>>
> >>>>>>
> >>>image
> >>>
> >>>
> >>>
> >>>
> >>>>at
> >>>>
> >>>>
> >>>>
> >>>>
> >>>>>>all...Then, I tried extracting a section from the moving slice
and
> >>>>>>registered them together....but the middle slices of the the
> >>>>>>
> >>>>>>
> >>>>>>
> >>>>>>
> >>>registered
> >>>
> >>>
> >>>
> >>>
> >>>>>>image does not match the corresponding slices in the fixed
image.
> >>>>>>
> >>>>>>
> >>>>>>
> >>>>>>
> >>>>What's
> >>>>
> >>>>
> >>>>
> >>>>
> >>>>>>going on?? Is there prerequisite on the input data for
performing
> >>>>>>mullti-modality registration using ITK?
> >>>>>>
> >>>>>>Please help!! I am in deperate need to really nail this this
time!!
> >>>>>>
> >>>>>>
> >>>>>>
> >>>>>>
> >>>>Thanx
> >>>>
> >>>>
> >>>>
> >>>>
> >>>>>>a million!!
> >>>>>>
> >>>>>>Grace
> >>>>>>
> >>>>>>----- Original Message -----
> >>>>>>From: "Luis Ibanez" <luis.ibanez at kitware.com>
> >>>>>>To: "Grace Chen" <Grace.Chen at swri.ca>
> >>>>>>Cc: <insight-users at itk.org>
> >>>>>>Sent: Thursday, May 25, 2006 9:58 PM
> >>>>>>Subject: Re: [Insight-users] 3d multi-modality registration
> >>>>>>
> >>>>>>
> >>>>>>
> >>>>>>
> >>>>>>
> >>>>>>
> >>>>>>>Hi Grace,
> >>>>>>>
> >>>>>>>In principle, that's true,
> >>>>>>>
> >>>>>>>but in practice it may not happen if the spatial extension
> >>>>>>>of the moving image doesn't fully overlap with the fixed
> >>>>>>>image.
> >>>>>>>
> >>>>>>>Have you checked the other slices ?
> >>>>>>>
> >>>>>>>Do they overlap well ?
> >>>>>>>
> >>>>>>>
> >>>>>>> Luis
> >>>>>>>
> >>>>>>>
> >>>>>>>===================
> >>>>>>>Grace Chen wrote:
> >>>>>>>
> >>>>>>>
> >>>>>>>
> >>>>>>>
> >>>>>>>>Hi there,
> >>>>>>>>
> >>>>>>>>My program performs the 3D multi-modality registration for two
> >>>>>>>>
> >>>>>>>>
> >>>>>>>>
> >>>>>>>>
> >>>>volumes.
> >>>>
> >>>>
> >>>>
> >>>>
> >>>>>>>>After the registration has been performed, is it true that the
> >>>>>>>>
> >>>>>>>>
> >>>>>>>>
> >>>>>>>>
> >>>first
> >>>
> >>>
> >>>
> >>>
> >>>>>>>>slice of the registered moving image should look like the
first
> >>>>>>>>
> >>>>>>>>
> >>>>>>>>
> >>>>>>>>
> >>>>slice of
> >>>>
> >>>>
> >>>>
> >>>>
> >>>>>>>>the fixed image?
> >>>>>>>>
> >>>>>>>>Grace
> >>>>>>>>
> >>>>>>>>
> >>>>>>>>
> >>>>>>>>
> >>>>>>>>
> >>>>>>>>
>
>>>>--------------------------------------------------------------------
----
> >>>>
> >>>>
> >>>>>
> >>>>>
> >>>>>>>>_______________________________________________
> >>>>>>>>Insight-users mailing list
> >>>>>>>>Insight-users at itk.org
> >>>>>>>>http://www.itk.org/mailman/listinfo/insight-users
> >>>>>>>>
> >>>>>>>>
> >>>>>>>>
> >>>>>>>>
> >>>>>>>
> >>>>>>>
> >>>>>>_______________________________________________
> >>>>>>Insight-users mailing list
> >>>>>>Insight-users at itk.org
> >>>>>>http://www.itk.org/mailman/listinfo/insight-users
> >>>>>>
> >>>>>>
> >>>>>>
> >>>>>>
> >>>>>_______________________________________________
> >>>>>Insight-users mailing list
> >>>>>Insight-users at itk.org
> >>>>>http://www.itk.org/mailman/listinfo/insight-users
> >>>>>
> >>>>>
> >>>>>
> >>>>>
> >>>>>
> >>>>_______________________________________________
> >>>>Insight-users mailing list
> >>>>Insight-users at itk.org
> >>>>http://www.itk.org/mailman/listinfo/insight-users
> >>>>
> >>>>
> >>>>
> >>>>
> >>>>
> >>>
> >>>
> >>>
> >>>
> >
> >
> >
> >
> _______________________________________________
> Insight-users mailing list
> Insight-users at itk.org
> http://www.itk.org/mailman/listinfo/insight-users
>
_______________________________________________
Insight-users mailing list
Insight-users at itk.org
http://www.itk.org/mailman/listinfo/insight-users
More information about the Insight-users
mailing list