7 #ifndef Tomato_OXModelT1Shmolli_HXX 8 #define Tomato_OXModelT1Shmolli_HXX 13 template<
typename MeasureType >
17 MeasureType A = parameters[0];
18 MeasureType B = parameters[1];
19 MeasureType T1star = parameters[2];
21 MeasureType calculated = 0;
23 if (_preventUnderOverFlow){
24 if (fabs(T1star) <= (time*0.005 ))
26 else if (fabs(T1star) >= (time*500 ))
27 calculated = (A - B*exp((MeasureType)-500));
36 if (fabs(T1star) < std::numeric_limits<MeasureType>::min()){
41 calculated = A - B * exp((MeasureType) -fabs( time / T1star ) );
43 calculated = A - B * exp((MeasureType) -time / T1star );
49 template<
typename MeasureType >
55 unsigned int nSamples = this->_nSamples;
57 MeasureType A = parameters[0];
58 MeasureType B = parameters[1];
59 MeasureType T1star = parameters[2];
61 for (
unsigned int i = 0; i < nSamples; i++) {
62 MeasureType invTime = this->_InvTimes[i];
63 MeasureType measured = this->_Signal[i];
64 MeasureType calculated = 0;
66 calculated = calcModelValue(parameters, invTime);
68 residuals[i] = calculated - measured;
72 template<
typename MeasureType >
76 int nSamples = this->_nSamples;
79 MeasureType B = parameters[1];
80 MeasureType T1star = parameters[2];
81 MeasureType invTime, myexp;
83 for (
int i = 0; i < nSamples; i++) {
84 invTime = this->_InvTimes[i];
85 myexp = exp(-invTime/T1star);
88 jacobian[i*3+0] = 1.0;
89 jacobian[i*3+1] = -myexp;
90 jacobian[i*3+2] = -B*invTime*myexp/(T1star*T1star);
94 template<
typename MeasureType >
100 int nSamples = this->_nSamples;
102 calcLSResiduals(parameters, this->_Residuals);
103 MeasureType result = 0;
105 if (_rootMedianSquareCost){
107 for (
int i = 0; i < nSamples; ++i) {
108 this->_Residuals[i] = this->_Residuals[i] * this->_Residuals[i];
111 result = KWUtil::calcMedianArray(nSamples, this->_Residuals);
116 for (
int i = 0; i < nSamples; ++i) {
117 result = result + this->_Residuals[i] * this->_Residuals[i];
124 if(parameters[0] < 0. )
125 result -= parameters[0];
126 if(parameters[0] > parameters[1])
127 result += (parameters[0] - parameters[1]);
128 if(parameters[1] > (4. * parameters[0]))
129 result += (parameters[1]-4. * parameters[0]);
130 if(parameters[2] < 0.)
131 result -= parameters[2];
139 template<
typename MeasureType >
145 int nSamples = this->_nSamples;
151 MeasureType measured, invTime, myexp;
153 MeasureType A = parameters[0];
154 MeasureType B = parameters[1];
155 MeasureType T1star = parameters[2];
158 for (
int i = 0; i < nSamples; ++i){
159 measured = this->getSignal()[i];
160 invTime = this->getInvTimes()[i];
161 myexp = exp(-invTime/T1star);
163 derivative[0] = derivative[0] + A*2 - measured*2 - myexp*B*2;;
164 derivative[1] = derivative[1] + myexp*2*(measured - A + myexp*B);
165 derivative[2] = derivative[2] + (invTime*myexp*B*2*(measured - A + myexp*B))/(T1star*T1star);
172 #endif //Tomato_OXModelT1Shmolli_HXX virtual MeasureType calcModelValue(const MeasureType *parameters, MeasureType time)
Definition: OxModelT1Shmolli.hxx:16
virtual void calcLSResiduals(const MeasureType *parameters, MeasureType *residuals)
Definition: OxModelT1Shmolli.hxx:52
virtual MeasureType calcCostValue(const MeasureType *parameters)
Definition: OxModelT1Shmolli.hxx:97
virtual void calcLSJacobian(const MeasureType *parameters, MeasureType *jacobian)
Definition: OxModelT1Shmolli.hxx:75
Definition: OxCalculator.h:19
virtual void calcCostDerivative(const MeasureType *parameters, MeasureType *derivative)
Definition: OxModelT1Shmolli.hxx:142