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