]> Creatis software - clitk.git/blob - registration/clitkGenericOptimizer.h
Add rotate and translate option to clitkAffineTransform
[clitk.git] / registration / clitkGenericOptimizer.h
1 /*=========================================================================
2   Program:   vv                     http://www.creatis.insa-lyon.fr/rio/vv
3
4   Authors belong to:
5   - University of LYON              http://www.universite-lyon.fr/
6   - Léon Bérard cancer center       http://www.centreleonberard.fr
7   - CREATIS CNRS laboratory         http://www.creatis.insa-lyon.fr
8
9   This software is distributed WITHOUT ANY WARRANTY; without even
10   the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
11   PURPOSE.  See the copyright notices for more information.
12
13   It is distributed under dual licence
14
15   - BSD        See included LICENSE.txt file
16   - CeCILL-B   http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html
17 ===========================================================================**/
18
19 #ifndef __clitkGenericOptimizer_h
20 #define __clitkGenericOptimizer_h
21
22 //clitk include
23 #include "clitkLBFGSBOptimizer.h"
24
25 //itk include
26 #include "itkAmoebaOptimizer.h"
27 #include "itkPowellOptimizer.h"
28 #include "itkFRPROptimizer.h"
29 #include "itkRegularStepGradientDescentOptimizer.h"
30 #include "itkVersorRigid3DTransformOptimizer.h"
31 #include "itkConjugateGradientOptimizer.h"
32 #include "itkLBFGSOptimizer.h"
33 #include "itkSPSAOptimizer.h"
34
35 // general
36 #include <iomanip>
37
38 /*Requires at least this section in the associated *.ggo file: (modify defaults to fit application)
39
40 option "optimizer"      -       "0=Simplex, 1=Powell, 2=FRPR, 3=Regular Step GD, 4=VersorRigid3D, 5=Conjugated Gradient, 6=L-BFGS, 7=L-BFGS-B" int no default="0"
41 option "delta"          -       "0: Initial delta, otherwise automatic"                                                 double  no
42 option "step"           -       "1,2,3,4: Initial stepsize (to be multiplied with the gradient)"                        double  no      default="2.0"
43 option "relax"          -       "3,4: Relaxation of the stepsize (multiplied each time the gradient changes sign)"      double  no      default="0.7"
44 option "valueTol"       -       "0,1,2: Tolerance on the function"                                                      double  no      default="0.01"
45 option "stepTol"        -       "0,1,3,4: Tolerance on the step size"                                                   double  no      default="0.1"
46 option "gradTol"        -       "3,4,6,7: Tolerance on the (projected) gradient magnitude (7: 1=low->1e-10=high precision)"     double  no      default="1e-5"
47 option "lineAcc"        -       "6: Line accuracy (eg: high=0.1, low=0.9)"      double  no      default="0.9"
48 option "convFactor"     -       "7: Convergence factor: terminate if factor*machine_precision>reduction in cost (1e+12 low -> 1e+1 high precision) "    double          no      default="1e+7"
49 option "maxIt"          -       "0-7: Maximum number of iterations"                                     int     no      default="500"
50 option "maxLineIt"      -       "1,2: Maximum number of line iterations"                                        int     no      default="50"
51 option "maxEval"        -       "6,7: Maximum number of evaluations"                                            int     no      default="500"
52 option "maxCorr"        -       "7: Maximum number of corrections"                                              int     no      default="5"
53 option "selectBound"    -       "7: Select the type of bound: 0=none, 1=u, 2=u&l, 3=l"                  int     no      default="0"
54 option "lowerBound"     -       "7: The lower bound"                                                    double  no
55 option "upperBound"     -       "7: The upper bound"                                                    double  no
56
57 */
58
59 namespace clitk
60 {
61 template<class args_info_type>
62 class GenericOptimizer : public itk::LightObject
63 {
64 public:
65   //==============================================
66   typedef GenericOptimizer     Self;
67   typedef itk::LightObject     Superclass;
68   typedef itk::SmartPointer<Self>            Pointer;
69   typedef itk::SmartPointer<const Self>      ConstPointer;
70
71   /** Method for creation through the object factory. */
72   itkNewMacro(Self);
73
74
75   typedef itk::SingleValuedNonLinearOptimizer OptimizerType;
76   typedef OptimizerType::Pointer OptimizerPointer;
77
78   //==============================================
79   //Set members
80   void SetArgsInfo(args_info_type args_info) {
81     m_ArgsInfo= args_info;
82     m_Verbose=m_ArgsInfo.verbose_flag;
83   }
84   void SetMaximize(bool f) {
85     m_Maximize=f;
86     m_MaximizeGiven=true;
87   }
88
89   void SetNumberOfParameters(unsigned int n) {
90     m_NumberOfParameters=n;
91   }
92
93   //==============================================
94   //Get members
95   OptimizerPointer GetOptimizerPointer(void) {
96
97     if (!m_MaximizeGiven)std::cerr<<"Warning: Maximize/Minimize was not given to the optimizer!!"<<std::endl;
98     else if(m_Verbose) {
99       if (m_Maximize)  std::cout<<"Maximizing similarity measure..."<<std::endl;
100       else std::cout<<"Minimizing similarity measure..."<<std::endl;
101     }
102
103
104     //============================================================================
105     //switch on  the  metric type chosen and set metric specific members
106     switch (m_ArgsInfo.optimizer_arg) {
107
108     case 0: {
109       itk::AmoebaOptimizer::Pointer m =
110         itk::AmoebaOptimizer::New();
111       m_OptimizerIsVNL=true;
112
113       //Set Parameters for this optimizer
114       m->SetMaximize(m_Maximize);
115       if (m_ArgsInfo.delta_given) {
116         itk::Array<double> delta(m_NumberOfParameters);
117         for (unsigned int i =0 ; i < m_NumberOfParameters; i ++)
118           delta[i]=m_ArgsInfo.delta_arg;
119         m->SetInitialSimplexDelta(delta);
120         m->SetAutomaticInitialSimplex(false);
121       } else m->SetAutomaticInitialSimplex(true);
122       m->SetMaximumNumberOfIterations(m_ArgsInfo.maxIt_arg);
123       m->SetFunctionConvergenceTolerance(m_ArgsInfo.valueTol_arg);
124       m->SetParametersConvergenceTolerance(m_ArgsInfo.stepTol_arg);
125       if (m_Verbose) std::cout<<"Using Simplex Optimizer..."<<std::endl;
126       m_Optimizer=m;
127       break;
128     }
129
130     case 1: {
131       itk::PowellOptimizer::Pointer m =
132         itk::PowellOptimizer::New();
133       m_OptimizerIsVNL=false;
134
135       //Set Parameters for this optimizer
136       m->SetMaximize(m_Maximize);
137       m->SetStepLength(m_ArgsInfo.step_arg);
138       m->SetStepTolerance(m_ArgsInfo.stepTol_arg);
139       m->SetValueTolerance (m_ArgsInfo.valueTol_arg);
140       m->SetMaximumIteration(m_ArgsInfo.maxIt_arg);
141       m->SetMaximumLineIteration(m_ArgsInfo.maxLineIt_arg);
142       if (m_Verbose) std::cout<<"Using Powell Optimizer..."<<std::endl;
143       m_Optimizer=m;
144       break;
145     }
146
147     case 2: {
148       itk::FRPROptimizer::Pointer m =
149         itk::FRPROptimizer::New();
150       m_OptimizerIsVNL=false;
151
152       //Set Parameters for this optimizer
153       m->SetMaximize(m_Maximize);
154       m->SetStepLength(m_ArgsInfo.step_arg);
155       m->SetStepTolerance(m_ArgsInfo.stepTol_arg);
156       m->SetValueTolerance (m_ArgsInfo.valueTol_arg);
157       m->SetMaximumIteration(m_ArgsInfo.maxIt_arg);
158       m->SetMaximumLineIteration(m_ArgsInfo.maxLineIt_arg);
159       if (m_Verbose) std::cout<<"Using Powell Optimizer..."<<std::endl;
160       m_Optimizer=m;
161       break;
162     }
163
164     case 3: {
165       itk::RegularStepGradientDescentOptimizer::Pointer m =
166         itk::RegularStepGradientDescentOptimizer::New();
167       m_OptimizerIsVNL=false;
168
169       //Set Parameters for this optimizer
170       m->SetMaximize(m_Maximize);
171       m->SetMinimumStepLength( m_ArgsInfo.stepTol_arg );
172       m->SetRelaxationFactor(  m_ArgsInfo.relax_arg );
173       m->SetMaximumStepLength( m_ArgsInfo.step_arg );
174       m->SetNumberOfIterations(m_ArgsInfo.maxIt_arg  );
175       m->SetGradientMagnitudeTolerance(m_ArgsInfo.gradTol_arg  );
176       if (m_Verbose) std::cout<<"Using Regular Step Gradient Descent Optimizer..."<<std::endl;
177       m_Optimizer=m;
178       break;
179     }
180
181     case 4: {
182       itk::VersorRigid3DTransformOptimizer::Pointer m
183       = itk::VersorRigid3DTransformOptimizer::New();
184       m_OptimizerIsVNL=false;
185
186       //Set Parameters for this optimizer
187       m->SetMaximize(m_Maximize);
188       m->SetMinimumStepLength( m_ArgsInfo.stepTol_arg );
189       m->SetRelaxationFactor(  m_ArgsInfo.relax_arg );
190       m->SetMaximumStepLength( m_ArgsInfo.step_arg );
191       m->SetNumberOfIterations(m_ArgsInfo.maxIt_arg  );
192       m->SetGradientMagnitudeTolerance(m_ArgsInfo.gradTol_arg  );
193       if (m_Verbose) std::cout<<"Using VersorRigid3DTransform Optimizer..."<<std::endl;
194       m_Optimizer=m;
195       break;
196     }
197
198     case 5: {
199
200       itk::ConjugateGradientOptimizer::Pointer m
201       = itk::ConjugateGradientOptimizer::New();
202       m_OptimizerIsVNL=true;
203
204       //Set Parameters for this optimizer
205       m->SetMaximize(m_Maximize);
206       if (m_Verbose) std::cout<<"Using Conjugated Gradient Optimizer..."<<std::endl;
207       m_Optimizer=m;
208       break;
209     }
210
211     case 6: {
212       itk::LBFGSOptimizer::Pointer m
213       = itk::LBFGSOptimizer::New();
214       m_OptimizerIsVNL=true;
215
216       //Set Parameters for this optimizer
217       m->SetMaximize(m_Maximize);
218       m->SetGradientConvergenceTolerance(m_ArgsInfo.gradTol_arg);
219       m->SetLineSearchAccuracy(m_ArgsInfo.lineAcc_arg);
220       m->SetMaximumNumberOfFunctionEvaluations(m_ArgsInfo.maxEval_arg);
221       m->SetTrace(true);
222       if (m_Verbose) std::cout<<"Using L-BFGS Optimizer..."<<std::endl;
223       m_Optimizer=m;
224       break;
225     }
226
227     case 7: {
228       LBFGSBOptimizer::Pointer m
229       = LBFGSBOptimizer::New();
230       m_OptimizerIsVNL=true;
231
232       //Set Parameters for this optimizer
233       m->SetMaximize(m_Maximize);
234       m->SetMaximumNumberOfEvaluations(m_ArgsInfo.maxEval_arg);
235       m->SetMaximumNumberOfIterations(m_ArgsInfo.maxIt_arg);
236       m->SetMaximumNumberOfCorrections(m_ArgsInfo.maxCorr_arg);
237       m->SetCostFunctionConvergenceFactor( m_ArgsInfo.convFactor_arg);
238       m->SetProjectedGradientTolerance( m_ArgsInfo.gradTol_arg);
239
240       // Bounds
241       clitk::LBFGSBOptimizer::BoundSelectionType boundSelect( m_NumberOfParameters );
242       clitk::LBFGSBOptimizer::BoundValueType upperBound( m_NumberOfParameters );
243       clitk::LBFGSBOptimizer::BoundValueType lowerBound( m_NumberOfParameters );
244       boundSelect.Fill( m_ArgsInfo.selectBound_arg );
245       upperBound.Fill( m_ArgsInfo.upperBound_arg );
246       lowerBound.Fill( m_ArgsInfo.lowerBound_arg );
247       m->SetBoundSelection( boundSelect );
248       m->SetUpperBound( upperBound );
249       m->SetLowerBound( lowerBound );
250
251       if (m_Verbose) std::cout<<"Using L-BFGS-B Optimizer..."<<std::endl;
252       m_Optimizer=m;
253       break;
254     }
255
256     //     case 8:
257     //       {
258     //  itk::SPSAOptimizer::Pointer m =
259     //    itk::SPSAOptimizerr::Pointer::New();
260     //      m_OptimizerIsVNL=false;
261
262     //  //Set Parameters for this optimizer
263     //  m->SetMaximize(m_Maximize);
264     //  m->SetGamma(m_ArgsInfo.gamma_arg);
265     //  m->SetAlpha(m_ArgsInfo.gamma_arg);
266     //  m->SetA(m_ArgsInfo.bigA_arg);
267     //  m->Seta(m_ArgsInfo.littleA_arg);
268     //  m->Setc(m_ArgsInfo.littleC_arg);
269     //  m->SetMaximumNumberOfIterations(m_ArgsInfo.maxIt_arg);
270     //  m->SetMinimumNumberOfIterations(m_ArgsInfo.minIt_arg);
271     //  m->SetNumberOfPerturbations(m_ArgsInfo.per_arg);
272     //  m->SetStateOfConvergenceDecayRate(args_inf.conv_arg);
273     //  if (m_Verbose) std::cout<<"Using SPSA Optimizer..."<<std::endl;
274     //  m_Optimizer=m;
275     //  break;
276     //       }
277
278     }
279
280
281     //============================================================================
282     //return the pointer
283     return m_Optimizer;
284   }
285
286
287   void OutputIterationInfo(void) const {
288     switch (m_ArgsInfo.optimizer_arg) {
289     case 0: {
290       itk::AmoebaOptimizer* o=dynamic_cast<itk::AmoebaOptimizer*>(m_Optimizer.GetPointer());
291       std::cout<<std::setprecision(6);
292       //if (m_OutputIteration)  std::cout<<"Iteration: "<< o->GetIteration()<<std::endl;
293       if (m_OutputPosition)  std::cout<<"Position: "<< o->GetCachedCurrentPosition()<<std::endl;
294       if (m_OutputValue)  std::cout<<"Value: "<< o->GetCachedValue()<<std::endl;
295       // if (m_OutputGradient)  std::cout<<"Gradient: "<< o->GetCachedDerivative()<<std::endl;
296       std::cout<<std::endl;
297       break;
298     }
299     case 1: {
300       itk::PowellOptimizer* o=dynamic_cast<itk::PowellOptimizer*>(m_Optimizer.GetPointer());
301       std::cout<<std::setprecision(6);
302       if (m_OutputIteration)  std::cout<<"Iteration: "<< o->GetCurrentIteration()<<std::endl;
303       if (m_OutputIteration)  std::cout<<"Line Iteration: "<< o->GetCurrentLineIteration()<<std::endl;
304       if (m_OutputPosition)  std::cout<<"Position: "<< o->GetCurrentPosition()<<std::endl;
305       if (m_OutputValue)  std::cout<<"Value: "<< o->GetCurrentCost()<<std::endl;
306       //if (m_OutputGradient)  std::cout<<"Gradient: "<< o->GetDerivative()<<std::endl;
307       std::cout<<std::endl;
308       break;
309     }
310     case 2: {
311       itk::FRPROptimizer* o=dynamic_cast<itk::FRPROptimizer*>(m_Optimizer.GetPointer());
312       std::cout<<std::setprecision(6);
313       if (m_OutputIteration)  std::cout<<"Iteration: "<< o->GetCurrentIteration()<<std::endl;
314       if (m_OutputIteration)  std::cout<<"Line Iteration: "<< o->GetCurrentLineIteration()<<std::endl;
315       if (m_OutputPosition)  std::cout<<"Position: "<< o->GetCurrentPosition()<<std::endl;
316       if (m_OutputValue)  std::cout<<"Value: "<< o->GetCurrentCost()<<std::endl;
317       //if (m_OutputGradient)  std::cout<<"Gradient: "<< o->GetDerivative()<<std::endl;
318       std::cout<<std::endl;
319       break;
320     }
321     case 3: {
322       itk::RegularStepGradientDescentOptimizer* o=dynamic_cast<itk::RegularStepGradientDescentOptimizer*>(m_Optimizer.GetPointer());
323       std::cout<<std::setprecision(6);
324       if (m_OutputIteration)  std::cout<<"Iteration: "<< o->GetCurrentIteration()<<std::endl;
325       if (m_OutputPosition)  std::cout<<"Position: "<< o->GetCurrentPosition()<<std::endl;
326       if (m_OutputValue)  std::cout<<"Step length: "<< o->GetCurrentStepLength()<<std::endl;
327       if (m_OutputValue)  std::cout<<"Value: "<< o->GetValue()<<std::endl;
328       if (m_OutputGradient)  std::cout<<"Gradient: "<< o->GetGradient()<<std::endl;
329       std::cout<<std::endl;
330       break;
331     }
332     case 4: {
333       itk::VersorRigid3DTransformOptimizer* o=dynamic_cast<itk::VersorRigid3DTransformOptimizer*>(m_Optimizer.GetPointer());
334       std::cout<<std::setprecision(6);
335       if (m_OutputIteration)  std::cout<<"Iteration: "<< o->GetCurrentIteration()<<std::endl;
336       if (m_OutputPosition)  std::cout<<"Position: "<< o->GetCurrentPosition()<<std::endl;
337       if (m_OutputPosition)  std::cout<<"Step length: "<< o->GetCurrentStepLength()<<std::endl;
338       if (m_OutputValue)  std::cout<<"Value: "<< o->GetValue()<<std::endl;
339       if (m_OutputGradient)  std::cout<<"Gradient: "<< o->GetGradient()<<std::endl;
340       std::cout<<std::endl;
341       break;
342     }
343     case 5: {
344       itk::ConjugateGradientOptimizer* o=dynamic_cast<itk::ConjugateGradientOptimizer*>(m_Optimizer.GetPointer());
345       std::cout<<std::setprecision(6);
346       if (m_OutputIteration)  std::cout<<"Iteration: "<< o->GetCurrentIteration()<<std::endl;
347       if (m_OutputPosition)  std::cout<<"Position: "<< o->GetCachedCurrentPosition()<<std::endl;
348       if (m_OutputValue)  std::cout<<"Value: "<< o->GetCachedValue()<<std::endl;
349       if (m_OutputGradient)  std::cout<<"Gradient: "<< o->GetCachedDerivative()<<std::endl;
350       std::cout<<std::endl;
351       break;
352     }
353     case 6: {
354       itk::LBFGSOptimizer* o=dynamic_cast<itk::LBFGSOptimizer*>(m_Optimizer.GetPointer());
355       std::cout<<std::setprecision(6);
356       //if (m_OutputIteration)  std::cout<<"Iteration: "<< o->GetCurrentIteration()<<std::endl;
357       if (m_OutputPosition)  std::cout<<"Position: "<< o->GetCachedCurrentPosition()<<std::endl;
358       if (m_OutputValue)  std::cout<<"Value: "<< o->GetCachedValue()<<std::endl;
359       //          if (m_OutputValue)  std::cout<<"Norm Projected Gradient: "<< o->GetInfinityNormOfProjectedGradient()<<std::endl;
360       if (m_OutputGradient)  std::cout<<"Gradient: "<< o->GetCachedDerivative()<<std::endl;
361       std::cout<<std::endl;
362       break;
363     }
364     case 7: {
365       clitk::LBFGSBOptimizer* o=dynamic_cast<clitk::LBFGSBOptimizer*>(m_Optimizer.GetPointer());
366       std::cout<<std::setprecision(6);
367       if (m_OutputIteration)  std::cout<<"Iteration: "<< o->GetCurrentIteration()<<std::endl;
368       if (m_OutputPosition)  std::cout<<"Position: "<< o->GetCachedCurrentPosition()<<std::endl;
369       if (m_OutputValue)  std::cout<<"Value: "<< o->GetValue()<<std::endl;
370       if (m_OutputValue)  std::cout<<"Norm Projected Gradient:"<< o->GetInfinityNormOfProjectedGradient()<<std::endl;
371       if (m_OutputGradient)  std::cout<<"Gradient: "<< o->GetCachedDerivative()<<std::endl;
372       std::cout<<std::endl;
373       break;
374     }
375     //       case 8:
376     //  {
377     //    itk::VersorRigid3DTransformOptimizer* o=dynamic_cast<itk::VersorRigid3DTransformOptimizer*>(m_Optimizer.GetPointer());
378     //    std::cout<<std::setprecision(6);
379     //    if (m_OutputIteration)  std::cout<<"Iteration: "<< o->GetIteration()<<std::endl;
380     //    if (m_OutputPosition)  std::cout<<"Position: "<< o->GetCurrentPosition()<<std::endl;
381     //    if (m_OutputValue)  std::cout<<"Value: "<< o->GetCurrentCost()<<std::endl;
382     //    if (m_OutputGradient)  std::cout<<"Gradient: "<< o->GetDerivative()<<std::endl;
383     //    std::cout<<std::endl;
384     //  }
385     }
386   }
387
388   void SetOutputIteration(bool o) {
389     m_OutputIteration=o;
390   }
391   void SetOutputPosition(bool o) {
392     m_OutputPosition=o;
393   }
394   void SetOutputValue(bool o) {
395     m_OutputValue=o;
396   }
397   void SetOutputGradient(bool o) {
398     m_OutputGradient=o;
399   }
400
401
402   //==============================================
403 protected:
404   GenericOptimizer() {
405     m_Optimizer=NULL;
406     m_Maximize=false;
407     m_Verbose=false;
408     m_MaximizeGiven=false;
409     m_NumberOfParameters=1;
410     m_OptimizerIsVNL=false;
411     m_OutputIteration=true;
412     m_OutputPosition=false;
413     m_OutputValue=true;
414     m_OutputGradient=false;
415
416   }
417   ~GenericOptimizer() {};
418
419 private:
420   args_info_type m_ArgsInfo;
421   OptimizerPointer m_Optimizer;
422   bool m_Maximize;
423   bool m_Verbose;
424   bool m_MaximizeGiven;
425   unsigned int m_NumberOfParameters;
426   bool m_OptimizerIsVNL;
427   bool m_OutputIteration;
428   bool m_OutputPosition;
429   bool m_OutputValue;
430   bool m_OutputGradient;
431 };
432
433 } // end namespace clitk
434
435 #endif // #define __clitkGenericOptimizer_h