1 /*=========================================================================
2  *
3  *  Copyright Insight Software Consortium
4  *
5  *  Licensed under the Apache License, Version 2.0 (the "License");
6  *  you may not use this file except in compliance with the License.
7  *  You may obtain a copy of the License at
8  *
9  *         http://www.apache.org/licenses/LICENSE-2.0.txt
10  *
11  *  Unless required by applicable law or agreed to in writing, software
12  *  distributed under the License is distributed on an "AS IS" BASIS,
13  *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14  *  See the License for the specific language governing permissions and
15  *  limitations under the License.
16  *
17  *=========================================================================*/
18 #ifndef _itkQuaternionRigidTransformGradientDescentOptimizer_hxx
19 #define _itkQuaternionRigidTransformGradientDescentOptimizer_hxx
20 
21 #include "itkQuaternionRigidTransformGradientDescentOptimizer.h"
22 #include "vnl/vnl_quaternion.h"
23 
24 namespace itk
25 {
26 /**
27  * Advance one Step following the gradient direction
28  */
29 void
30 QuaternionRigidTransformGradientDescentOptimizer
AdvanceOneStep()31 ::AdvanceOneStep()
32 {
33   const double direction = ( m_Maximize ) ? 1.0 : -1.0;
34   const ScalesType & scales = this->GetScales();
35 
36   const unsigned int spaceDimension =  m_CostFunction->GetNumberOfParameters();
37 
38   // Make sure the scales have been set
39   if ( scales.size() != spaceDimension )
40     {
41     itkExceptionMacro(<< "The size of Scales is "
42                       << scales.size()
43                       << ", but the NumberOfParameters is "
44                       << spaceDimension
45                       << ".");
46     }
47 
48   DerivativeType transformedGradient(spaceDimension);
49   for ( unsigned int i = 0; i < spaceDimension; i++ )
50     {
51     transformedGradient[i] = m_Gradient[i] / scales[i];
52     }
53 
54   ParametersType currentPosition = this->GetCurrentPosition();
55 
56   // compute new quaternion value
57   vnl_quaternion< double > newQuaternion;
58   for ( unsigned int j = 0; j < 4; j++ )
59     {
60     newQuaternion[j] = currentPosition[j] + direction * m_LearningRate
61                        * transformedGradient[j];
62     }
63 
64   newQuaternion.normalize();
65 
66   ParametersType newPosition(spaceDimension);
67   // update quaternion component of currentPosition
68   for ( unsigned int j = 0; j < 4; j++ )
69     {
70     newPosition[j] = newQuaternion[j];
71     }
72 
73   // update the translation component
74   for ( unsigned int j = 4; j < spaceDimension; j++ )
75     {
76     newPosition[j] = currentPosition[j]
77                      + direction * m_LearningRate * transformedGradient[j];
78     }
79 
80   // First invoke the event, so the current position
81   // still corresponds to the metric values.
82   this->InvokeEvent( IterationEvent() );
83 
84   this->SetCurrentPosition(newPosition);
85 }
86 } // end namespace itk
87 
88 #endif
89