MantisBT - ITK
View Issue Details
IDProjectCategoryView StatusDate SubmittedLast Update
0004827ITKpublic2007-04-10 09:302008-06-03 04:44
ReporterTom Vercauteren 
Assigned ToStephen Aylward 
PriorityhighSeveritymajorReproducibilityalways
StatusclosedResolutionfixed 
PlatformOSOS Version
Product Version 
Target VersionFixed in Version 
Resolution Date
Sprint
Sprint Status
Summary0004827: Line search routine fails in itkPowellOptimizer.cxx
DescriptionI have successfully been using the FRPR optimizer for my registration
problems for a while now.

The FRPROptimizer class depends on PowellOptimizer.

Since the replacement in itkPowellOptimizer.cxx of the code from
numerical recipes by an open-source implementation, I am having
problems.

Specifically, the line search procedure does not always optimize the
metric anymore.

Below is an example of the evolution of the metric (minus correlation
coefficient squared) I observe on an 2D rigid body registration
problem. I show the iteration number, the current metric and the
current position. Note that with the old code (from ITK 2.2 to ITK
3.0.1), an IterationEvent was invoked for each line iteration, so that
one iteration appears several time.

Old Code (from ITK 2.2 to ITK 3.0.1):
==========================
0 -0.769804 [-0.0168971, 25.7384, 15.1015]
0 -0.769804 [-0.0168971, 25.7384, 15.1015]
0 -0.770011 [-0.0140876, 25.6878, 15.0462]
0 -0.770057 [-0.0143894, 25.6933, 15.0521]
0 -0.770059 [-0.015085, 25.7058, 15.0658]
0 -0.770059 [-0.015085, 25.7058, 15.0658]
0 -0.770059 [-0.015085, 25.7058, 15.0658]
0 -0.770059 [-0.015085, 25.7058, 15.0658]
0 -0.770059 [-0.015085, 25.7058, 15.0658]
0 -0.77006 [-0.015186, 25.7076, 15.0678]
0 -0.77006 [-0.015186, 25.7076, 15.0678]
0 -0.77006 [-0.015186, 25.7076, 15.0678]
0 -0.77006 [-0.015186, 25.7076, 15.0678]
0 -0.77006 [-0.0151956, 25.7078, 15.068]
0 -0.77006 [-0.0151972, 25.7078, 15.068]
0 -0.77006 [-0.0151972, 25.7078, 15.068]
1 -0.77091 [-0.0195981, 26.3341, 15.7091]
1 -0.77091 [-0.0195981, 26.3341, 15.7091]
1 -0.770982 [-0.0188923, 26.2337, 15.6063]
1 -0.770982 [-0.0188923, 26.2337, 15.6063]
1 -0.770989 [-0.0190474, 26.2558, 15.6289]
1 -0.770989 [-0.0190474, 26.2558, 15.6289]
1 -0.770989 [-0.0190474, 26.2558, 15.6289]
1 -0.770989 [-0.0190474, 26.2558, 15.6289]
1 -0.770989 [-0.0190474, 26.2558, 15.6289]
1 -0.770989 [-0.0190474, 26.2558, 15.6289]
1 -0.770989 [-0.0190474, 26.2558, 15.6289]
1 -0.770989 [-0.0190474, 26.2558, 15.6289]
1 -0.770989 [-0.0190474, 26.2558, 15.6289]
1 -0.770989 [-0.0190474, 26.2558, 15.6289]
2 -0.772777 [-0.0142605, 27.2713, 16.7834]
2 -0.773462 [-0.0160889, 26.8834, 16.3424]
2 -0.773559 [-0.0159783, 26.9069, 16.3691]
2 -0.773559 [-0.0159783, 26.9069, 16.3691]
2 -0.773581 [-0.0157697, 26.9512, 16.4194]
2 -0.773588 [-0.0158524, 26.9336, 16.3995]
2 -0.773634 [-0.0158404, 26.9362, 16.4024]
2 -0.773641 [-0.0158119, 26.9422, 16.4092]
2 -0.773645 [-0.0157958, 26.9456, 16.4131]
2 -0.773645 [-0.0157958, 26.9456, 16.4131]
2 -0.773645 [-0.0157958, 26.9456, 16.4131]
2 -0.773646 [-0.015792, 26.9464, 16.414]
2 -0.773647 [-0.0157896, 26.9469, 16.4146]
2 -0.773647 [-0.0157896, 26.9469, 16.4146]
3 -0.773647 [-0.0157896, 26.9469, 16.4146]
3 -0.773898 [-0.0129291, 26.8306, 16.361]
3 -0.773976 [-0.0137756, 26.865, 16.3769]
3 -0.774031 [-0.013866, 26.8687, 16.3786]
3 -0.774031 [-0.013866, 26.8687, 16.3786]
3 -0.774031 [-0.013866, 26.8687, 16.3786]
3 -0.774031 [-0.013866, 26.8687, 16.3786]
3 -0.774031 [-0.013866, 26.8687, 16.3786]
3 -0.774032 [-0.0138554, 26.8683, 16.3784]
3 -0.774032 [-0.0138554, 26.8683, 16.3784]
3 -0.774032 [-0.0138554, 26.8683, 16.3784]
3 -0.774032 [-0.0138554, 26.8683, 16.3784]
3 -0.774032 [-0.0138554, 26.8683, 16.3784]
3 -0.774032 [-0.0138554, 26.8683, 16.3784]
4 -0.774032 [-0.0138554, 26.8683, 16.3784]
4 -0.774032 [-0.0138554, 26.8683, 16.3784]
4 -0.774032 [-0.0138554, 26.8683, 16.3784]
4 -0.774032 [-0.0138554, 26.8683, 16.3784]
4 -0.774032 [-0.0138554, 26.8683, 16.3784]
4 -0.774032 [-0.0138554, 26.8683, 16.3784]
4 -0.774032 [-0.0138554, 26.8683, 16.3784]
4 -0.774032 [-0.0138544, 26.8683, 16.3784]
4 -0.774032 [-0.0138544, 26.8683, 16.3784]
4 -0.774032 [-0.0138544, 26.8683, 16.3784]
4 -0.774032 [-0.0138544, 26.8683, 16.3784]
5 -0.774032 [-0.0138544, 26.8683, 16.3784]
5 -0.774032 [-0.0138544, 26.8683, 16.3784]
5 -0.774032 [-0.0138544, 26.8683, 16.3784]
5 -0.774032 [-0.0138544, 26.8683, 16.3784]
5 -0.774032 [-0.0138544, 26.8683, 16.3784]
5 -0.774032 [-0.0138544, 26.8683, 16.3784]
5 -0.774032 [-0.0138544, 26.8683, 16.3784]
5 -0.774032 [-0.0138544, 26.8683, 16.3784]
5 -0.774032 [-0.0138544, 26.8683, 16.3784]
5 -0.774032 [-0.0138544, 26.8683, 16.3784]
5 -0.774032 [-0.0138544, 26.8683, 16.3784]
6 -0.774032 [-0.0138544, 26.8683, 16.3784]
6 -0.774032 [-0.0138544, 26.8683, 16.3784]
6 -0.774032 [-0.0138544, 26.8683, 16.3784]
6 -0.774032 [-0.0138544, 26.8683, 16.3784]
6 -0.774293 [-0.0137035, 26.332, 16.3866]
6 -0.774407 [-0.0137361, 26.4477, 16.3848]
6 -0.774441 [-0.0137594, 26.5307, 16.3835]
6 -0.774441 [-0.0137594, 26.5307, 16.3835]
6 -0.774441 [-0.0137594, 26.5307, 16.3835]
7 -0.774441 [-0.0137594, 26.5307, 16.3835]
7 -0.774441 [-0.0137594, 26.5307, 16.3835]
7 -0.774441 [-0.0137594, 26.5307, 16.3835]
7 -0.774441 [-0.0137594, 26.5307, 16.3835]
7 -0.774441 [-0.0137594, 26.5307, 16.3835]
7 -0.774441 [-0.0137594, 26.5307, 16.3835]
7 -0.774444 [-0.0137552, 26.5418, 16.3835]
7 -0.774448 [-0.0137527, 26.5487, 16.3835]
7 -0.77445 [-0.0137511, 26.553, 16.3835]
7 -0.774451 [-0.0137497, 26.5565, 16.3835]
7 -0.774451 [-0.0137497, 26.5565, 16.3835]
8 -0.774451 [-0.0137497, 26.5565, 16.3835]
8 -0.774451 [-0.0137497, 26.5565, 16.3835]
8 -0.774451 [-0.0137497, 26.5565, 16.3835]
8 -0.774451 [-0.0137497, 26.5565, 16.3835]
8 -0.774451 [-0.0137497, 26.5565, 16.3835]
8 -0.774451 [-0.0137497, 26.5565, 16.3835]
8 -0.774451 [-0.0137497, 26.5565, 16.3835]
8 -0.774451 [-0.0137497, 26.5565, 16.3835]
8 -0.774451 [-0.0137497, 26.5565, 16.3835]
8 -0.774451 [-0.0137434, 26.5573, 16.3836]
8 -0.774451 [-0.0137396, 26.5578, 16.3837]
8 -0.774452 [-0.0137372, 26.5581, 16.3838]
8 -0.774452 [-0.0137372, 26.5581, 16.3838]
9 -0.774452 [-0.0137372, 26.5581, 16.3838]
9 -0.774452 [-0.0137372, 26.5581, 16.3838]
9 -0.774452 [-0.0137372, 26.5581, 16.3838]
9 -0.774452 [-0.0137372, 26.5581, 16.3838]
9 -0.774452 [-0.0137372, 26.5581, 16.3838]
9 -0.774452 [-0.0137372, 26.5581, 16.3838]
9 -0.774452 [-0.0137372, 26.5581, 16.3838]
9 -0.774452 [-0.0137372, 26.5581, 16.3838]
9 -0.774452 [-0.0137372, 26.5581, 16.3838]
9 -0.774452 [-0.0137372, 26.5581, 16.3838]
9 -0.774452 [-0.0137372, 26.5581, 16.3838]
10 -0.774452 [-0.0137372, 26.5581, 16.3838]
10 -0.774452 [-0.0137372, 26.5581, 16.3838]
10 -0.774452 [-0.0137372, 26.5581, 16.3838]
10 -0.774452 [-0.0137372, 26.5581, 16.3838]
10 -0.774452 [-0.0137372, 26.5581, 16.3838]
10 -0.774452 [-0.0137372, 26.5581, 16.3838]
10 -0.774452 [-0.0137372, 26.5581, 16.3838]
10 -0.774452 [-0.0137372, 26.5581, 16.3838]
10 -0.774452 [-0.0137372, 26.5581, 16.3838]
11 -0.774452 [-0.0137372, 26.5581, 16.3838]
11 -0.774452 [-0.0137372, 26.5581, 16.3838]
11 -0.774452 [-0.0137372, 26.5581, 16.3838]
11 -0.774452 [-0.0137372, 26.5581, 16.3838]
11 -0.774452 [-0.0137372, 26.5581, 16.3838]
11 -0.774452 [-0.0137372, 26.5581, 16.3838]
11 -0.774452 [-0.0137372, 26.5581, 16.3838]
11 -0.774452 [-0.0137372, 26.5581, 16.3838]


New Code (from ITK 3.2):
==================
0 -0.769804 [-0.0168971, 25.7384, 15.1015]
1 -0.769837 [-0.0167342, 25.7253, 15.0865]
2 -0.771157 [-0.0119474, 25.8603, 15.2465]
3 -0.772384 [-0.00949709, 26.199, 15.716]
4 -0.773165 [-0.00986638, 26.5148, 16.4328]
5 -0.773165 [-0.00986638, 26.5148, 16.4328]
6 -0.773165 [-0.00986638, 26.5148, 16.4328]
7 -0.773169 [-0.0098776, 26.5201, 16.4522]
8 -0.773317 [-0.0101162, 26.4585, 16.2409]
9 -0.77335 [-0.0105894, 26.5696, 16.8105]
10 -0.773414 [-0.0106002, 26.5713, 16.8194]
11 -0.774268 [-0.013778, 26.5061, 16.6986]
12 -0.774325 [-0.0137774, 26.5069, 16.672]
13 -0.774364 [-0.0133232, 26.5137, 16.5445]
14 -0.774395 [-0.013719, 26.5101, 16.6186]
15 -0.774426 [-0.0138309, 26.5132, 16.5643]
16 -0.774426 [-0.0138309, 26.5132, 16.5643]
17 -0.774426 [-0.0138309, 26.5132, 16.5643]
18 -0.774426 [-0.0138309, 26.5132, 16.5643]
19 -0.0404348 [-0.14626, 148.764, 129.749]
20 -0.0501792 [-0.127289, 134.878, 116.987]
21 -0.0521457 [-0.128033, 139.064, 120.873]
22 -0.0608914 [-0.0356036, 143.264, 126.153]
23 -0.0610972 [-0.023946, 143.146, 126.355]
24 -0.0611195 [-0.0241748, 143.15, 126.352]
25 -0.0611195 [-0.0241748, 143.15, 126.352]
26 -0.0618224 [0.00164232, 142.13, 126.99]
27 -0.0618529 [-0.000840117, 142.214, 126.937]
28 -0.061793 [-0.0090229, 142.361, 126.84]
29 -0.0622588 [-0.0189155, 142.549, 126.691]
30 -0.0622588 [-0.0189189, 142.55, 126.691]
31 -0.0624076 [-0.000970234, 141.906, 127.156]
32 -0.0626338 [-0.0138566, 142.314, 126.866]
33 -0.0626375 [-0.0141352, 142.32, 126.862]
34 -0.0626376 [-0.0141538, 142.32, 126.861]
35 -0.0626439 [-0.0138298, 142.309, 126.871]
36 -0.0626439 [-0.0138298, 142.309, 126.871]
37 -0.0626439 [-0.0138298, 142.309, 126.871]
38 -0.0713286 [0.00488931, 133.753, 127.396]
39 -0.0713313 [0.00484949, 133.775, 127.395]
40 -0.0713314 [0.00485464, 133.776, 127.396]
41 -0.0713315 [0.00490002, 133.776, 127.398]
42 -0.0713315 [0.00490114, 133.776, 127.398]
43 -0.0713315 [0.00490114, 133.776, 127.398]
44 -0.0160748 [-0.193237, 257.032, 117.469]
45 -0.0808865 [-0.0480491, 213.052, 119.576]
46 -0.0875697 [-0.0139101, 205.018, 120.018]
47 -0.0889239 [-0.0171101, 206.342, 120]
48 -0.157121 [0.223701, 216.295, 124.404]
49 -0.164172 [0.232336, 212.851, 114.938]
50 -0.164213 [0.232305, 212.836, 114.894]
51 -0.165212 [0.2276, 212.64, 114.229]
52 -0.165285 [0.227216, 212.617, 114.251]
53 -0.165519 [0.223366, 212.497, 114.091]
54 -0.16557 [0.223824, 212.503, 114.118]
55 -0.16557 [0.223824, 212.503, 114.118]
56 -0.165639 [0.224715, 212.584, 114.246]
57 -0.165697 [0.224897, 212.601, 114.273]
58 -0.165698 [0.224891, 212.601, 114.272]
59 -0.165698 [0.224889, 212.601, 114.272]
60 -0.165698 [0.224889, 212.601, 114.272]
61 -0.165702 [0.225115, 212.744, 114.285]
62 -0.165702 [0.225113, 212.743, 114.284]
63 -0.165842 [0.222427, 212.661, 114.147]
64 -0.165846 [0.222191, 212.659, 114.137]
65 -0.165846 [0.222191, 212.659, 114.137]
66 -0.165846 [0.222191, 212.659, 114.137]
67 -0.165846 [0.222191, 212.659, 114.137]
Steps To Reproduce
Additional Information
TagsNo tags attached.
Relationships
Attached Fileszip TestFRPR.zip (285,497) 1969-12-31 19:00
https://public.kitware.com/Bug/file/1010/TestFRPR.zip
? itkPowellOptimizer.h (7,945) 2007-10-18 12:48
https://public.kitware.com/Bug/file/1193/itkPowellOptimizer.h
cxx itkPowellOptimizer.cxx (11,902) 2007-10-18 12:48
https://public.kitware.com/Bug/file/1194/itkPowellOptimizer.cxx
Issue History
Date ModifiedUsernameFieldChange
2007-10-18 12:48Tom VercauterenFile Added: itkPowellOptimizer.h
2007-10-18 12:48Tom VercauterenFile Added: itkPowellOptimizer.cxx
2007-10-18 12:49Tom VercauterenNote Added: 0009510
2008-04-29 11:29Luis IbanezAssigned ToLuis Ibanez => Stephen Aylward
2008-04-30 11:07Stephen AylwardStatusassigned => resolved
2008-04-30 11:07Stephen AylwardResolutionopen => fixed
2008-04-30 11:07Stephen AylwardNote Added: 0011602
2008-06-03 04:44Tom VercauterenStatusresolved => closed

Notes
(0007695)
Luis Ibanez   
2007-05-19 07:16   
From Tom's email on 04/16/07:

Hi Luis,

Please find attached a small example that shows the problem. I didn't
include the CMakeLists file as I don't usually build my stuff within
the ITK tree. I didn't submit the code to the insight journal because
it is not a decent quality code ;)

You can simply run it as:
TestFRPR -f colon_fix.mhd -m colon_mov.mhd

You can try to build it with the open-source line search and with the
numerical recipes one and see that the open-source one diverges.

Best,
Tom
(0009510)
Tom Vercauteren   
2007-10-18 12:48   
The problem seems to arise from the LineBracket function in
itk::PowellOptimizer.cxx

I saw several bugs in this class. First of all there are two unused
functions Swap and Shift.

Then in the first part of LineBracket, the goal seems to be to get an
initial bracket. This is done by a recursive
 *x3 = *x1 + goldenRatio * ( *x2 - *x1 );
From what I have seen in the litterature, it should be
 *x3 = *x2 + goldenRatio * ( *x2 - *x1 );

In the second part of this function, it is said that a Golden Section
search is used. I am not sure why this is useful as the bracket is
always followed by a line search. It is also not used in NRC for
example. Finally, this seems buggy as it assumes fa > fb < fc while
from the if above we can have fc > fa < fb. Instead of an if, a swap
should for example be used.

Attached are two patched files that fixes the problems. My solution
was to simply use the first initial bracketing and skipping the Golden
Section search.
(0011602)
Stephen Aylward   
2008-04-30 11:07   
Committed fix on 2008.04.29