Tomato
OxModelT1TwoParam.hxx
Go to the documentation of this file.
1 
7 #ifndef Tomato_OXModelT1TwoParam_HXX
8 #define Tomato_OXModelT1TwoParam_HXX
9 
10 namespace Ox {
11 
12 
13  template< typename MeasureType >
14  MeasureType
16  ::calcModelValue(const MeasureType* parameters, MeasureType time){
17 
18  MeasureType A = parameters[0];
19  MeasureType T1 = parameters[1];
20 
21  if (fabs(T1) < std::numeric_limits<MeasureType>::min())
22  return A;
23 
24  return A * (1 - exp( -time / T1 ));
25  }
26 
27  template< typename MeasureType >
28  void
30  ::calcLSResiduals(const MeasureType* parameters, MeasureType* residuals){
31 
32  unsigned int nSamples = this->_nSamples;
33 
34  for (unsigned int i = 0; i < nSamples; i++) {
35  MeasureType invTime = this->_InvTimes[i];
36  MeasureType measured = this->_Signal[i];
37  MeasureType calculated = 0;
38 
39  calculated = calcModelValue(parameters, invTime);
40 
41  residuals[i] = calculated - measured;
42  }
43  }
44 
45  template< typename MeasureType >
46  void
48  ::calcLSJacobian(const MeasureType* parameters, MeasureType* jacobian){
49 
50  int nSamples = this->_nSamples;
51 
52  MeasureType A = parameters[0];
53  MeasureType T1 = parameters[1];
54  MeasureType invTime, myexp;
55 
56  for (int i = 0; i < nSamples; i++) {
57  invTime = this->_InvTimes[i];
58  myexp = exp(-invTime/T1);
59 
60  // calculated in matlab (syms A T1 t), f=A*(1-exp(-t./T1)); diff(f,A), diff(f,T1)
61  jacobian[i*2+0] = 1.0 - myexp;
62  jacobian[i*2+1] = -( A * invTime * myexp )/ ( T1 * T1 );
63  }
64  }
65 
66  template< typename MeasureType >
67  MeasureType
69  ::calcCostValue(const MeasureType* parameters){
70 
71  int nSamples = this->_nSamples;
72 
73  calcLSResiduals(parameters, this->_Residuals);
74  MeasureType result = 0;
75 
76  for (int i = 0; i < nSamples; ++i) {
77  result = result + this->_Residuals[i] * this->_Residuals[i];
78  }
79 
80  return result;
81  }
82 
83  template< typename MeasureType >
84  void
86  ::calcCostDerivative(const MeasureType* parameters, MeasureType* derivative){
87 
88  int nSamples = this->_nSamples;
89 
90  derivative[0] = 0;
91  derivative[1] = 0;
92 
93  MeasureType measured, invTime, myexp;
94 
95  MeasureType A = parameters[0];
96  MeasureType T1 = parameters[1];
97 
98  // calculated in matlab (syms A B T1 t y), f=(A-B*exp(-t./T1)-y).^2; diff(f,A), diff(f,B), diff(f,T1)
99  for (int i = 0; i < nSamples; ++i){
100  measured = this->getSignal()[i];
101  invTime = this->getInvTimes()[i];
102  myexp = exp(-invTime / T1);
103 
104  derivative[0] = derivative[0] + 2 * (myexp - 1)*(measured + A*(myexp - 1));
105  derivative[1] = derivative[1] + ( 2 * A * invTime * myexp * ( measured + A * ( myexp - 1 ))) / ( T1 * T1 );
106  }
107  }
108 
109 } //namespace Ox
110 
111 
112 #endif //Tomato_OXModelT1TwoParam_HXX
virtual MeasureType calcCostValue(const MeasureType *parameters)
Definition: OxModelT1TwoParam.hxx:69
virtual MeasureType calcModelValue(const MeasureType *parameters, MeasureType time)
Definition: OxModelT1TwoParam.hxx:16
virtual void calcCostDerivative(const MeasureType *parameters, MeasureType *derivative)
Definition: OxModelT1TwoParam.hxx:86
virtual void calcLSResiduals(const MeasureType *parameters, MeasureType *residuals)
Definition: OxModelT1TwoParam.hxx:30
Definition: OxCalculator.h:19
virtual void calcLSJacobian(const MeasureType *parameters, MeasureType *jacobian)
Definition: OxModelT1TwoParam.hxx:48