[Insight-users] Registration - Mutual Information + Affine - Fine tuning of parameters for optimization
Sharath Venkatesha
sharath20284 at yahoo.com
Sat Jul 4 02:31:11 EDT 2009
Hello,
I am using MutualInformationHistogramImageToImageMetric + RegularStepGradientDescentOptimizer combination, and have a few doubts on it.
I understand that the value of the metric should be maximum at a location of optimal registration. But I also expected that the value of the gradient magnitude to be minimum at this location, which is not so (The gradient here is computed as per GetDerivative() function line 186 in itkHistogramImageToImageMetric.txx). I expect this behavior, as once the metric value is maximized (I ensured that this is not a local maxima) the value of the gradient should be least, and there should also be a change in direction. I do see the conical shape in the metric values, but the optimizer always overshoots it and ends at a location where metric value is lower. I have observed it in all my test cases.
For example in the log that I have provided below, around iteration 45, the metric value has peaked (starting from 0.1749 in iteration 31, peak of 0.1808 in iteration 45, and then the value going down afterwards), but the magnitude of the gradient (as defined in line 205, itkRegularStepGradientDescentBaseOptimizer.cxx), I get at this iteration is 5.909, which is not minimum. I also see a pattern of the value of the gradient magnitude to be gradually gradually decreasing.
I do not understand this behavior. Can you please throw some light?
*** I see that the values of m_DerivativeStepLengthScales in function GetDerivative() of itkHistogramImageToImageMetric.txx file, all the values are set to 1.0 ,using
m_DerivativeStepLengthScales.Fill(1.0);
This means that at each step, when the gradient is computed (i.e the parameters are perturbed by m_DerivativeStepLength), all the parameters are assumed to have uniform scaling. How can this be assumed here?
Shouldn't these be scaled appropriately, like the optimizer scale
values used in
RegularStepGradientDescentBaseOptimizer::AdvanceOneStep() function?
Also the value of m_DerivativeStepLength is always set to 1.0 . Shouldn't this a user defined parameter, or based on the optimizer ?
--- LOG ----
Below is the log I am printing to analyze the value of the metric and how the optimizer behaves
MaximumStepLength : 0.1
MinimumStepLength : 0.0001
......
Format:
Gradient_Magnitude Tolerance_value Current_step_Length
Iter_Number Metric_Value [Affine Parameters] Angle_in_Degrees Change_in_Metric_value Change_in_angle
........
.....
GradMag = 7.494 Tol =0.0001 StepLen = 0.1
31 0.1749 [0.5827, -0.06824, 0.06624, 0.6027, -856.8, -931.7] Angle:6.473 Metr:0.0005967 ADiff:-0.02686
GradMag = 8.051 Tol =0.0001 StepLen = 0.1
32 0.1757 [0.5827, -0.06827, 0.06578, 0.6035, -856.7, -931.7] Angle:6.448 Metr:0.0007982 ADiff:-0.02441
GradMag = 7.914 Tol =0.0001 StepLen = 0.1
33 0.1766 [0.5826, -0.06834, 0.06529, 0.6043, -856.6, -931.7] Angle:6.423 Metr:0.0009113 ADiff:-0.02468
GradMag = 7.315 Tol =0.0001 StepLen = 0.1
34 0.1773 [0.5826, -0.06841, 0.06477, 0.6052, -856.5, -931.7] Angle:6.398 Metr:0.0006155 ADiff:-0.02567
GradMag = 7.99 Tol =0.0001 StepLen = 0.1
35 0.1778 [0.5827, -0.06849, 0.06428, 0.606, -856.4, -931.6] Angle:6.374 Metr:0.0004911 ADiff:-0.02425
GradMag = 7.949 Tol =0.0001 StepLen = 0.1
36 0.1783 [0.5827, -0.06857, 0.06378, 0.6068, -856.4, -931.6] Angle:6.349 Metr:0.000521 ADiff:-0.02433
GradMag = 7.78 Tol =0.0001 StepLen = 0.1
37 0.1787 [0.5827, -0.06868, 0.06327, 0.6076, -856.3, -931.5] Angle:6.326 Metr:0.0004448 ADiff:-0.02362
GradMag = 7.599 Tol =0.0001 StepLen = 0.1
38 0.179 [0.5828, -0.06881, 0.06276, 0.6085, -856.2, -931.5] Angle:6.302 Metr:0.0003036 ADiff:-0.02324
GradMag = 7.277 Tol =0.0001 StepLen = 0.1
39 0.1795 [0.5829, -0.06896, 0.0622, 0.6093, -856.1, -931.5] Angle:6.278 Metr:0.0004738 ADiff:-0.02439
GradMag = 8.089 Tol =0.0001 StepLen = 0.1
40 0.1797 [0.5829, -0.0691, 0.06171, 0.6101, -856, -931.4] Angle:6.257 Metr:0.0002499 ADiff:-0.02088
GradMag = 8.046 Tol =0.0001 StepLen = 0.1
41 0.18 [0.583, -0.06926, 0.06122, 0.6108, -855.9, -931.4] Angle:6.237 Metr:0.0002045 ADiff:-0.01961
GradMag = 8.101 Tol =0.0001 StepLen = 0.1
42 0.1802 [0.5831, -0.06943, 0.06075, 0.6116, -855.8, -931.4] Angle:6.219 Metr:0.0002648 ADiff:-0.01865
GradMag = 7.21 Tol =0.0001 StepLen = 0.1
43 0.1803 [0.5831, -0.06964, 0.06021, 0.6124, -855.7, -931.3] Angle:6.199 Metr:6.091e-005 ADiff:-0.02007
GradMag = 8.228 Tol =0.0001 StepLen = 0.1
44 0.1805 [0.5832, -0.06986, 0.05974, 0.6131, -855.6, -931.3] Angle:6.183 Metr:0.0002261 ADiff:-0.01615
GradMag = 7.485 Tol =0.0001 StepLen = 0.1
45 0.1808 [0.5833, -0.07012, 0.05922, 0.6139, -855.5, -931.2] Angle:6.166 Metr:0.0002553 ADiff:-0.01641
GradMag = 5.909 Tol =0.0001 StepLen = 0.1
46 0.1806 [0.5834, -0.07048, 0.05857, 0.6148, -855.4, -931.2] Angle:6.147 Metr:-0.0002007 ADiff:-0.01916
GradMag = 6.526 Tol =0.0001 StepLen = 0.1
47 0.1803 [0.5835, -0.07083, 0.058, 0.6156, -855.3, -931.1] Angle:6.132 Metr:-0.000293 ADiff:-0.01502
GradMag = 7.555 Tol =0.0001 StepLen = 0.1
48 0.1803 [0.5836, -0.07113, 0.05752, 0.6163, -855.3, -931.1] Angle:6.12 Metr:3.284e-006 ADiff:-0.01187
GradMag = 6.259 Tol =0.0001 StepLen = 0.1
49 0.1803 [0.5837, -0.0715, 0.05695, 0.617, -855.2, -931.1] Angle:6.106 Metr:-3.978e-006 ADiff:-0.01394
GradMag = 6.12 Tol =0.0001 StepLen = 0.1
50 0.1803 [0.5837, -0.07192, 0.05638, 0.6177, -855.1, -931] Angle:6.095 Metr:6.792e-005 ADiff:-0.01128
GradMag = 5.317 Tol =0.0001 StepLen = 0.1
51 0.1803 [0.5839, -0.07242, 0.05577, 0.6185, -855, -931] Angle:6.086 Metr:-2.758e-005 ADiff:-0.009191
GradMag = 5.665 Tol =0.0001 StepLen = 0.1
52 0.1802 [0.584, -0.07289, 0.05525, 0.6191, -854.9, -930.9] Angle:6.079 Metr:-0.0001544 ADiff:-0.006347
GradMag = 5.083 Tol =0.0001 StepLen = 0.1
53 0.1802 [0.5842, -0.07338, 0.05469, 0.6198, -854.8, -930.9] Angle:6.072 Metr:3.066e-006 ADiff:-0.007603
GradMag = 4.135 Tol =0.0001 StepLen = 0.1
54 0.1802 [0.5844, -0.07398, 0.05404, 0.6206, -854.7, -930.8] Angle:6.064 Metr:2.267e-005 ADiff:-0.007948
GradMag = 3.486 Tol =0.0001 StepLen = 0.1
55 0.1797 [0.5848, -0.07466, 0.05324, 0.6215, -854.7, -930.8] Angle:6.052 Metr:-0.0005252 ADiff:-0.01135
GradMag = 4.202 Tol =0.0001 StepLen = 0.1
56 0.1792 [0.585, -0.07522, 0.05258, 0.6222, -854.6, -930.7] Angle:6.043 Metr:-0.0004813 ADiff:-0.009771
GradMag = 3.866 Tol =0.0001 StepLen = 0.1
57 0.1789 [0.5853, -0.07577, 0.05185, 0.6229, -854.5, -930.7] Angle:6.029 Metr:-0.0002625 ADiff:-0.01337
GradMag = 3.972 Tol =0.0001 StepLen = 0.1
58 0.1786 [0.5856, -0.07628, 0.05115, 0.6235, -854.4, -930.6] Angle:6.017 Metr:-0.0003571 ADiff:-0.01262
GradMag = 4.536 Tol =0.0001 StepLen = 0.1
59 0.1781 [0.5858, -0.07669, 0.05057, 0.6239, -854.3, -930.5] Angle:6.005 Metr:-0.0004517 ADiff:-0.0113
GradMag = 3.345 Tol =0.0001 StepLen = 0.1
60 0.1777 [0.5861, -0.07723, 0.04978, 0.6244, -854.3, -930.4] Angle:5.99 Metr:-0.000372 ADiff:-0.0155
-....
Thanks for any clues,
Sharath Venkatesha
________________________________
From: Luis Ibanez <luis.ibanez at kitware.com>
To: Sharath Venkatesha <sharath20284 at yahoo.com>
Cc: Insight users <insight-users at itk.org>
Sent: Thursday, July 2, 2009 4:28:59 PM
Subject: Re: [Insight-users] Registration - Mutual Information + Affine - Fine tuning of parameters for optimization
Hi Sharath,
0) It is great that you have verified the Initialization
1) Translation initialization is "usually" enough...
but only if the rotation and scaling misalignments
are small....
A typical registration will not be able to correct
for more than a 30 degrees rotation, or a scaling
beyond a factor of 1.5.
So, please provide all the initialization that you can.
You may want to consider the use of the
LandmarkBasedTransformInitializer
2) If you results start diverging from the expected answer
in the very first iterations, then you should suspect:
a) Incorrect balance of the rotation vs translation
values in the parameter scaling array
or
b) Your step lenghts are too large in the optimizer
or
c) You may have set the Optimizer to Maximize, when
the Metric needs to be Minimized (or the other way
around).
3) You are correct that if the Registration is walking in the
right direction, the Joint histograms should be getting
more concentrated in the diagonal...
BUT ONLY IF
you are working with two images of the same modality.
Otherwise, if the images are of different modality, the
histogram will never be concentrated in the diagonal.
4) Unfortunately we don't have much material written on
how to analyze the histograms.
I will suggest that you save them as images using the
Histogram to Image filter, and then load them as a
series of 2D+T images into VV
http://www.creatis.insa-lyon.fr/rio/vv
so that you can play them as an animation.
5) The typical use of the vtkRegistrationMonitor will be:
#include "vtkRegistrationMonitor.h"
#include "vtkKWImage.h"
#include "vtkContourFilter.h"
#include "vtkImageData.h"
vtkRegistrationMonitor monitor;
vtkSmartPointer< vtkKWImage > fixedKWImage = vtkSmartPointer< vtkKWImage >::New();
vtkSmartPointer< vtkKWImage > movingKWImage = vtkSmartPointer< vtkKWImage >::New();
fixedKWImage->SetITKImageBase( const_cast<FixedImageType *>( fixedImage ) );
movingKWImage->SetITKImageBase( const_cast<MovingImageType *>( movingImage ) );
vtkSmartPointer< vtkContourFilter > fixedContour = vtkSmartPointer< vtkContourFilter >::New();
vtkSmartPointer< vtkContourFilter > movingContour = vtkSmartPointer< vtkContourFilter >::New();
fixedContour->SetInput( fixedKWImage->GetVTKImage() );
movingContour->SetInput( movingKWImage->GetVTKImage() );
fixedContour->SetValue( 0, 200.0 ); // level for iso-contour (DEPENDS on your data)
movingContour->SetValue( 0, 200.0 ); // level for iso-contour (DEPENDS on your data)
monitor.SetFixedSurface( fixedContour->GetOutput() );
monitor.SetMovingSurface( movingContour->GetOutput() );
monitor.SetNumberOfIterationPerUpdate( 1 ); (RECONSIDER CHANGING...)
std::string screenshotDirectory = argv[7];
itksys::SystemTools::MakeDirectory( screenshotDirectory.c_str() );
std::cout << screenshotDirectory << std::endl;
monitor.SetScreenshotOutputDirectory( screenshotDirectory.c_str() );
monitor.SetScreenshotBaseName( "registrationScreenshot" );
monitor.Observe( optimizer, rigidTransform );
You will find the vtkKWImage classes also in
InsightApplications/Auxiliary/vtk
---
Regards,
Luis
----------------------------------------------------------------------
On Thu, Jul 2, 2009 at 3:01 PM, Sharath Venkatesha <sharath20284 at yahoo.com> wrote:
>>Hi,
>
>>Thanks for the all the help and information provided in the earlier mails.
>
>>I have ensured that I am using correct initialization ( translation only) and correct values of optimizer scales.
>
>>*** Is it sufficient if I provide initialization parameters for translation (approx values) only? It is difficult for me estimate the scale and rotation parameters.
>
>>I am using MutualInformationHistogramImageToImageMetric +
>>MultiResolutionImageRegistrationMethod + AffineTransform +
>>RegularStepGradientDescentOptimizer/GradientDescentOptimizer
>
>>I am having problem with tuning of parameters (MaxStepSize and MinStepSize for RegularStepGradientDescentOptimizer, and LearningRate for GradientDescentOptimizer) . My results start diverging off the correct route in the first few iterations itself.
>
>>I have looked into
>>(1) Plotting of joint histograms (section 8.5.3 of ITK manual)
>>(2) vtkRegistrationMonitor
>
>
>>I am having trouble interpreting the results of the joint histograms.I understand that with correct parameters for registration, I should get a histogram which has highest density only along the diagonal. I have the outputs of the JointHistograms after very iteration for very level , and finding hard to follow the changes.
>
>>*** Is there is a defined procedure for understanding the results?
>
>>Can you please point me to where I can get more information on (1) and (2) ?
>
>>Thanks,
>>Sharath Venkatesha
>
>
>>------------------------------
>>Luis wrote:
>
>
>>Hi Sharath,
>
>>0) Whenever you get the exception saying that too many points
>> mapped outside of the moving image, it means that the
>> current Transform is such that when mapping the moving image
>> into the Fixed image coordinate system the overlap between
>> the two image is so small that it is unlikely that the
>> registration will recover in further iterations.
>
>> This is typically due to:
>
>
>> A) Poor initilization of the Transform
>
>> B) Poor selection of Scaling parameters
>> (the array that normalizes the dynamic
>> range of the different Transform parameters,
>> e.g. radians versus millimeters)
>
>> C) Optimizers that are set to perform jumps
>> that are too large, and bring the Transform
>> out of the range of the image.
>
>>1) You want to check this potential suspects in order.
>
>> That is.
>
>> First, verify that the initial transform
>> is reasonably placing the Moving image on top of the
>> Fixed image. You can do this by using the Resample
>> image filter, passing the moving image as input,
>> using the initial transform as Transform, and using
>> the Fixed image as reference. Then compare the
>> fixed image to the resampled moving image.
>
>> If the initial image looks ok,
>> then you want to check the values of the ParameterScaling.
>> It should be such that when you look at the Transform
>> parameters at ever iteration (using an Observer), the
>> values should change from iteration to iteration, according
>> to the expected dynamic range.
>
>> For example, transform parameters that correspond to rotations
>> should change by increments smaller than 0.01 (since they are
>> measured in Radians). While transform parameters that correspond
>> to translations should change at increments of 1 ~ 10
>
>
>> Finally, you should identify the parameter of the optimzer
>> that is responsible for selecting the size of jumps that
>> are performed in the parameteric space. (e.g. as you have
>> done for the learning rate in the GradientDescent optimizer).
>
>> You want to reduce the size of that jump, until you get the
>> Transform to have small increments at every iteration.
>
>
>>2) These parameters must be set up for every "family" of
>> registration problems. That is, the parameters that may be
>> good for registering a T1 to T2 MRI brain images, may not be
>> appropriate for registering a confocal microscopy image to
>> another.
>
>> However, once you fine tune the paramters for a pair of
>> T1-T2 images, it is likely that the same set of parameters
>> will work for another pair of the same type.
>
>> There is a need for a "smart layer" above the registration
>> framework, that could take away from the user the burden
>> of finding proper parameter settings....
>
>> any ideas are welcome :-)
>
>
>>3) Visual monitoring of the registraiton process
>> will help to make the fine-tunning process less
>> frustrating.
>
>> You may want to give it a try at the VTK helper
>> classes:
>
>> InsightApplications/Auxiliary/vtk/
>> vtkRegistrationMonitor.h
>> vtkRegistrationMonitor.cxx
>
>> they will display renderings from iso-surfaces
>> at every iteration of the registration process.
>
>> This is usually very informative...
>
>
>
>> Regards,
>
>
>> Luis
>
>
>>--------------
>>sharath v wrote:
>>> Hi,
>>>
>>> Thanks for the help. Changing the learning parameter worked...
>>>
>>>
>>For Viola MI + Affine and with learning rate of 0.01 and 100
>>iterations, I get good results on BrainProtonDensitySliceR10X13Y17
>>image. Whereas on the BrainProtonDensitySliceR10X13Y17S12 image, it
>>requires atleast 200 iterations to give correct results.
>>>
>>>
>>I want to use an optimizer which has a stopping value (i.e not fixed
>>number of iterations) like Amoeba/Evolutionary/GradientDescentStep
>>>
>>> I tried using the Amoeba optimizer with the following
>>>
>>> OptimizerType::ParametersType simplexDelta( transform->GetNumberOfParameters() );
>>> simplexDelta.Fill( 5.0 );
>>> optimizer->AutomaticInitialSimplexOff();
>>> optimizer->SetInitialSimplexDelta( simplexDelta );
>>> optimizer->SetParametersConvergenceTolerance( 0.01 ); // quarter pixel
>>> optimizer->SetFunctionConvergenceTolerance(0.001); // 0.1%
>>> optimizer->SetMaximumNumberOfIterations( 200 );
>>>
>>>
>>but I get an exception that sampled point mapped to outside of the
>>moving image, after 6-7 iterations. Similar issue happens for
>>OnePlusOne optimizer with
>>>
>>> typedef itk::Statistics::NormalVariateGenerator GeneratorType;
>>> GeneratorType::Pointer generator = GeneratorType::New();
>>> generator->Initialize(12345);
>>> optimizer->SetNormalVariateGenerator( generator );
>>> optimizer->Initialize( 10 );
>>> optimizer->SetEpsilon( 1.0 );
>>> optimizer->SetMaximumIteration( 4000 );
>>>
>>> Can you please let me know what values need to be used?
>>>
>>> And is there a way to make the registration process partially independent of these parameters?
>>>
>>> Thanks,
>>> Sharath
>
>
>
>
>
>>_____________________________________
>>Powered by www.kitware.com
>
>>Visit other Kitware open-source projects at
>http://www.kitware.com/opensource/opensource.html
>
>>Please keep messages on-topic and check the ITK FAQ at: http://www.itk.org/Wiki/ITK_FAQ
>
>>Follow this link to subscribe/unsubscribe:
>http://www.itk.org/mailman/listinfo/insight-users
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://www.itk.org/pipermail/insight-users/attachments/20090703/2db4829f/attachment-0001.htm>
More information about the Insight-users
mailing list