7 #ifndef TOMATO_OXCALCULATORT1SHMOLLI_HXX 8 #define TOMATO_OXCALCULATORT1SHMOLLI_HXX 12 template<
typename MeasureType >
18 if (this->getNSamples() != 7){
19 throw std::runtime_error(
"T1_SHMOLLI reconstruction available only for 7 samples.");
28 if (this->getFitter()->getMaxFunctionEvals() == 0){
33 for (
int i = 0; i < this->getNSamples()-1; i++){
34 if (this->getInvTimes()[i] > this->getInvTimes()[i+1]){
35 throw std::runtime_error(
"InvTimes have to be sorted!");
40 this->getSignCalculator()->setNSamples(this->getNSamples());
41 this->getSignCalculator()->setInvTimes(this->getInvTimes());
42 this->getSignCalculator()->setSigMag(this->getSigMag());
43 this->getSignCalculator()->setSigPha(this->getSigPha());
44 this->getSignCalculator()->setSignal(this->_Signal);
45 this->getSignCalculator()->setSigns(this->_Signs);
47 this->getSignCalculator()->calculateSign();
50 this->getStartPointCalculator()->setNDims(this->_nDims);
51 if (!this->getStartPointCalculator()->getInputStartPoint()){
52 if (this->_nDims == 3){
53 MeasureType
const temp[] = {100, 200, 1000};
54 this->getStartPointCalculator()->setInputStartPoint(temp);
56 throw std::runtime_error(
"Calculator: Set InputStartPoint in StartPointCalculator");
59 this->getStartPointCalculator()->setNSamples(this->getNSamples());
60 this->getStartPointCalculator()->setInvTimes(this->getInvTimes());
61 this->getStartPointCalculator()->setSigMag(this->getSigMag());
62 this->getStartPointCalculator()->setSigns(this->getSigns());
63 this->getStartPointCalculator()->setCalculatedStartPoint(this->_StartPoint);
65 this->getStartPointCalculator()->calculateStartPoint();
68 calculateTRRaverageHB();
73 template<
typename MeasureType >
77 const MeasureType * sigMag5,
const MeasureType * sigMag6,
const MeasureType * sigMag7,
78 const MeasureType * signs5,
const MeasureType * signs6,
const MeasureType * signs7){
80 int nStartFitPoints = 0;
83 startPointCalculator->setCalculatedStartPoint(this->_StartPoint);
85 if (this->getSigPha()){
86 startPointCalculator->setNSamples(5);
87 startPointCalculator->setInvTimes(invTimes5);
88 startPointCalculator->setSigMag(sigMag5);
89 startPointCalculator->setSigns(signs5);
92 if( nStartFitPoints < 3 ) {
93 startPointCalculator->setNSamples(6);
94 startPointCalculator->setInvTimes(invTimes6);
95 startPointCalculator->setSigMag(sigMag6);
96 startPointCalculator->setSigns(signs6);
99 if( nStartFitPoints < 3 ) {
100 startPointCalculator->setNSamples(7);
101 startPointCalculator->setInvTimes(invTimes7);
102 startPointCalculator->setSigMag(sigMag7);
103 startPointCalculator->setSigns(signs7);
107 if( nStartFitPoints < 3 ) {
108 startPointCalculator->setNSamples(7);
109 startPointCalculator->setInvTimes(invTimes7);
110 startPointCalculator->setSigMag(sigMag7);
111 startPointCalculator->setSigns(signs7);
119 template<
typename MeasureType >
125 const MeasureType *invTimes = this->getInvTimes();
126 MeasureType invTimesDifferences[4];
127 if (this->_nSamples == 7) {
128 invTimesDifferences[0] = invTimes[3] - invTimes[0];
129 invTimesDifferences[1] = invTimes[4] - invTimes[3];
130 invTimesDifferences[2] = invTimes[5] - invTimes[4];
131 invTimesDifferences[3] = invTimes[6] - invTimes[5];
134 _TRRaverageHB = KWUtil::calcMeanArray(4, invTimesDifferences);
140 template<
typename MeasureType >
145 this->_Results = std::map <std::string, MeasureType>();
148 if (KWUtil::calcMeanArray(this->getNSamples(), this->getSigMag()) < this->getMeanCutOff()) {
153 if(this->prepareToCalculate() == 1){
158 const double MaxTIForSignInvert = (this->MAX_T1_TRESHOLD * 0.67);
159 MeasureType lastValue = 1e32;
160 MeasureType lastValueTemp = lastValue;
168 double chiTemp = 1e32;
169 unsigned int nShmolliSamplesUsed = 0;
170 MeasureType TRRaverageHB = this->_TRRaverageHB;
172 std::map <std::string, MeasureType> results0, results5, results6, results7;
173 MeasureType signal5[5], signal6[6], signal7[7];
174 MeasureType invTimes5[5], invTimes6[6], invTimes7[7];
175 MeasureType signs5[5], signs6[6], signs7[7];
177 getShmolliSamples(this->_InvTimes, invTimes5, invTimes6, invTimes7);
178 getShmolliSamples(this->_Signal, signal5, signal6, signal7);
179 getShmolliSamples(this->_Signs, signs5, signs6, signs7);
181 this->getStartPointSKPShmolli(
182 invTimes5, invTimes6, invTimes7,
183 signal5, signal6, signal7,
184 signs5, signs6, signs7);
187 signal5[0] = -fabs(signal5[0]);
190 results5 = this->calculateWithSignCheck(5, invTimes5, signal5, signs5);
192 if (results5[
"A"] > 1){
193 nShmolliSamplesUsed = 5;
194 T1temp = results5[
"T1"];
195 chiTemp = results5[
"ChiSqrt"];
197 chiTemp = results5[
"LastValue"];
201 results6 = this->calculateWithSignCheck(6, invTimes6, signal6, signs6);
203 if ((T1temp <= TRRaverageHB)
204 || (results6[
"A"] <= 1)
205 || (results6[
"T1"] <= 0.4 * TRRaverageHB)
209 if ((results6[
"A"] > 1)
210 && (results6[
"ChiSqrt"] * results6[
"T1"] < chiTemp * TRRaverageHB))
212 nShmolliSamplesUsed = 6;
213 T1temp = results6[
"T1"];
217 results7 = this->calculateWithSignCheck(7, invTimes7, signal7, signs7);
219 if ((results7[
"A"] > 1)
220 && (results7[
"ChiSqrt"] * results7[
"T1"] < results6[
"ChiSqrt"] * TRRaverageHB * 0.4)
221 && (T1temp < TRRaverageHB)
224 nShmolliSamplesUsed = 7;
229 if (nShmolliSamplesUsed == 5) {
230 this->_Results = results5;
232 else if (nShmolliSamplesUsed == 6) {
233 this->_Results = results6;
235 else if (nShmolliSamplesUsed == 7) {
236 this->_Results = results7;
238 this->_Results = results0;
241 this->_Results[
"NShmolliSamplesUsed"] = nShmolliSamplesUsed;
242 this->_Results[
"hasBeenCalculated"] = 1;
243 this->_ParametersAfterFitting[0] = this->_Results[
"A"];
244 this->_ParametersAfterFitting[1] = this->_Results[
"B"];
245 this->_ParametersAfterFitting[2] = this->_Results[
"T1star"];
250 template<
typename MeasureType >
253 ::getShmolliSamples(
const MeasureType* inArray, MeasureType* outArray5, MeasureType* outArray6, MeasureType* outArray7){
254 int indexes7[] = {0,1,2,3,4,5,6};
255 int indexes6[] = {0,1,3,4,5,6};
256 int indexes5[] = {0,3,4,5,6};
257 for (
int i = 0; i < 7; i++) outArray7[i] = inArray[indexes7[i]];
258 for (
int i = 0; i < 6; i++) outArray6[i] = inArray[indexes6[i]];
259 for (
int i = 0; i < 5; i++) outArray5[i] = inArray[indexes5[i]];
265 #endif //TOMATO_OXCALCULATORT1SHMOLLI_HXX int calculateTRRaverageHB()
Definition: OxCalculatorT1Shmolli.hxx:122
virtual int calculate()
Definition: OxCalculatorT1Shmolli.hxx:143
virtual int prepareToCalculate()
Definition: OxCalculatorT1Shmolli.hxx:15
Definition: OxCalculatorT1Shmolli.h:22
Definition: OxStartPointCalculator.h:21
Definition: OxCalculator.h:19
virtual int calculateStartPoint()=0