13 m_ExcitationFlipAngle = M_PI / 2.0;
17 double flipAngle,
double m0Value)
19 m_SimulatedT2Values.set_size(3 * m_NumberOfEchoes + 1,m_NumberOfEchoes + 1);
20 m_OutputVector.resize(m_NumberOfEchoes);
24 double baseValue = m0Value * std::sin(m_ExcitationFlipAngle);
25 m_SimulatedT2Values(0,0) = baseValue;
26 for (
unsigned int i = 1;i <= m_NumberOfEchoes;++i)
27 m_SimulatedT2Values(0,i) = 0.0;
30 for (
unsigned int i = 0;i < m_NumberOfEchoes;++i)
33 m_SimulatedT2Values(i+1,0) = m_FirstEPGProduct * m_SimulatedT2Values(i,0) - m_SecondEPGProduct * m_SimulatedT2Values(i,3);
34 if (m_NumberOfEchoes > 1)
35 m_SimulatedT2Values(i+1,0) += m_ThirdEPGProduct * m_SimulatedT2Values(i,5);
38 m_SimulatedT2Values(i+1,1) = m_FourthEPGProduct * m_SimulatedT2Values(i,2);
39 m_SimulatedT2Values(i+1,2) = m_FirstEPGProduct * m_SimulatedT2Values(i,1);
40 m_SimulatedT2Values(i+1,3) = m_FifthEPGProduct * m_SimulatedT2Values(i,3) - m_SecondEPGProduct * m_SimulatedT2Values(i,0) / 2.0;
42 if (m_NumberOfEchoes > 1)
44 m_SimulatedT2Values(i+1,2) -= m_SecondEPGProduct * m_SimulatedT2Values(i,6);
45 m_SimulatedT2Values(i+1,3) += m_SecondEPGProduct * m_SimulatedT2Values(i,5) / 2.0;
47 if (m_NumberOfEchoes > 2)
48 m_SimulatedT2Values(i+1,2) += m_ThirdEPGProduct * m_SimulatedT2Values(i,8);
53 for (j = 1;j < m_NumberOfEchoes - 1;++j)
55 m_SimulatedT2Values(i+1,1 + j * 3) = m_SecondEPGProduct * m_SimulatedT2Values(i,3 + (j - 1) * 3) + m_FirstEPGProduct * m_SimulatedT2Values(i,2 + j * 3);
57 m_SimulatedT2Values(i+1,1 + j * 3) += m_ThirdEPGProduct * m_SimulatedT2Values(i,1 + (j - 2) * 3);
59 m_SimulatedT2Values(i+1,1 + j * 3) += m_ThirdEPGProduct * m_SimulatedT2Values(i,0);
61 m_SimulatedT2Values(i+1,2 + j * 3) = m_FirstEPGProduct * m_SimulatedT2Values(i,1 + j * 3);
62 m_SimulatedT2Values(i+1,3 + j * 3) = m_FifthEPGProduct * m_SimulatedT2Values(i,3 + j * 3) - m_SecondEPGProduct * m_SimulatedT2Values(i,1 + (j - 1) * 3) / 2.0;
64 if ((j + 1) < m_NumberOfEchoes)
66 m_SimulatedT2Values(i+1,2 + j * 3) -= m_SecondEPGProduct * m_SimulatedT2Values(i,3 + (j + 1) * 3);
67 m_SimulatedT2Values(i+1,3 + j * 3) += m_SecondEPGProduct * m_SimulatedT2Values(i,2 + (j + 1) * 3) / 2.0;
69 if ((j + 2) < m_NumberOfEchoes)
70 m_SimulatedT2Values(i+1,2 + j * 3) += m_ThirdEPGProduct * m_SimulatedT2Values(i,2 + (j + 2) * 3);
75 j = m_NumberOfEchoes - 1;
76 m_SimulatedT2Values(i+1,1 + j * 3) = m_SecondEPGProduct * m_SimulatedT2Values(i,3 + (j - 1) * 3) + m_FirstEPGProduct * m_SimulatedT2Values(i,2 + j * 3);
77 m_SimulatedT2Values(i+1,2 + j * 3) = 0.0;
78 m_SimulatedT2Values(i+1,3 + j * 3) = m_FifthEPGProduct * m_SimulatedT2Values(i,3 + j * 3) - m_SecondEPGProduct * m_SimulatedT2Values(i,1 + (j - 1) * 3) / 2.0;
81 m_SimulatedT2Values(i+1,1 + j * 3) += m_ThirdEPGProduct * m_SimulatedT2Values(i,1 + (j - 2) * 3);
83 m_SimulatedT2Values(i+1,1 + j * 3) += m_ThirdEPGProduct * m_SimulatedT2Values(i,0);
85 m_OutputVector[i] = m_SimulatedT2Values(i+1,0);
88 return m_OutputVector;
94 double espT2Value = 0.0;
96 espT2Value = std::exp(- m_EchoSpacing / (2 * t2Value));
98 double espT1Value = 0.0;
100 espT1Value = std::exp(- m_EchoSpacing / (2 * t1Value));
102 double cosB1alpha = std::cos(flipAngle);
103 double cosB1alpha2 = std::cos(flipAngle / 2.0);
104 double sinB1alpha = std::sin(flipAngle);
105 double sinB1alpha2 = std::sin(flipAngle / 2.0);
107 m_FirstEPGProduct = sinB1alpha2 * sinB1alpha2 * espT2Value * espT2Value;
108 m_SecondEPGProduct = sinB1alpha * espT1Value * espT2Value;
109 m_ThirdEPGProduct = cosB1alpha2 * cosB1alpha2 * espT2Value * espT2Value;
110 m_FourthEPGProduct = espT2Value * espT2Value;
111 m_FifthEPGProduct = cosB1alpha * espT1Value * espT1Value;
113 m_FirstDerivativeProduct = cosB1alpha2 * sinB1alpha2 * espT2Value * espT2Value;
114 m_SecondDerivativeProduct = cosB1alpha * espT1Value * espT2Value;
115 m_ThirdDerivativeProduct = - sinB1alpha * espT1Value * espT1Value;
120 m_OutputB1Derivative.resize(m_NumberOfEchoes);
122 m_SimulatedDerivativeT2Values.set_size(3 * m_NumberOfEchoes + 1,m_NumberOfEchoes + 1);
123 for (
unsigned int i = 0;i <= m_NumberOfEchoes;++i)
124 m_SimulatedDerivativeT2Values(0,i) = 0.0;
127 for (
unsigned int i = 0;i < m_NumberOfEchoes;++i)
131 m_SimulatedDerivativeT2Values(i+1,0) = m_FirstDerivativeProduct * m_SimulatedT2Values(i,0) - m_SecondDerivativeProduct * m_SimulatedT2Values(i,3);
134 m_SimulatedDerivativeT2Values(i+1,1) = 0.0;
135 m_SimulatedDerivativeT2Values(i+1,2) = m_FirstDerivativeProduct * m_SimulatedT2Values(i,1);
136 m_SimulatedDerivativeT2Values(i+1,3) = m_ThirdDerivativeProduct * m_SimulatedT2Values(i,3) - m_SecondDerivativeProduct * m_SimulatedT2Values(i,0) / 2.0;
138 if (m_NumberOfEchoes > 1)
140 m_SimulatedDerivativeT2Values(i+1,0) -= m_FirstDerivativeProduct * m_SimulatedT2Values(i,5);
141 m_SimulatedDerivativeT2Values(i+1,2) -= m_SecondDerivativeProduct * m_SimulatedT2Values(i,6);
143 if (m_NumberOfEchoes > 2)
144 m_SimulatedDerivativeT2Values(i+1,2) -= m_FirstDerivativeProduct * m_SimulatedT2Values(i,8);
146 m_SimulatedDerivativeT2Values(i+1,3) += m_SecondDerivativeProduct * m_SimulatedT2Values(i,5) / 2.0;
151 for (j = 1;j < m_NumberOfEchoes - 1;++j)
153 m_SimulatedDerivativeT2Values(i+1,1 + j * 3) = m_SecondDerivativeProduct * m_SimulatedT2Values(i,3 + (j - 1) * 3) + m_FirstDerivativeProduct * m_SimulatedT2Values(i,2 + j * 3);
155 m_SimulatedDerivativeT2Values(i+1,1 + j * 3) -= m_FirstDerivativeProduct * m_SimulatedT2Values(i,1 + (j - 2) * 3);
157 m_SimulatedDerivativeT2Values(i+1,1 + j * 3) -= m_FirstDerivativeProduct * m_SimulatedT2Values(i,0);
159 m_SimulatedDerivativeT2Values(i+1,2 + j * 3) = m_FirstDerivativeProduct * m_SimulatedT2Values(i,1 + j * 3);
160 m_SimulatedDerivativeT2Values(i+1,3 + j * 3) = m_ThirdDerivativeProduct * m_SimulatedT2Values(i,3 + j * 3) - m_SecondDerivativeProduct * m_SimulatedT2Values(i,1 + (j - 1) * 3) / 2.0;
162 if ((j + 1) < m_NumberOfEchoes)
164 m_SimulatedDerivativeT2Values(i+1,2 + j * 3) -= m_SecondDerivativeProduct * m_SimulatedT2Values(i,3 + (j + 1) * 3);
165 m_SimulatedDerivativeT2Values(i+1,3 + j * 3) += m_SecondDerivativeProduct * m_SimulatedT2Values(i,2 + (j + 1) * 3) / 2.0;
167 if ((j + 2) < m_NumberOfEchoes)
168 m_SimulatedDerivativeT2Values(i+1,2 + j * 3) -= m_FirstDerivativeProduct * m_SimulatedT2Values(i,2 + (j + 2) * 3);
173 j = m_NumberOfEchoes - 1;
174 m_SimulatedDerivativeT2Values(i+1,1 + j * 3) = m_SecondDerivativeProduct * m_SimulatedT2Values(i,3 + (j - 1) * 3) + m_FirstDerivativeProduct * m_SimulatedT2Values(i,2 + j * 3);
176 m_SimulatedDerivativeT2Values(i+1,1 + j * 3) -= m_FirstDerivativeProduct * m_SimulatedT2Values(i,1 + (j - 2) * 3);
178 m_SimulatedDerivativeT2Values(i+1,1 + j * 3) -= m_FirstDerivativeProduct * m_SimulatedT2Values(i,0);
180 m_SimulatedDerivativeT2Values(i+1,2 + j * 3) = 0.0;
181 m_SimulatedDerivativeT2Values(i+1,3 + j * 3) = m_ThirdDerivativeProduct * m_SimulatedT2Values(i,3 + j * 3) - m_SecondDerivativeProduct * m_SimulatedT2Values(i,1 + (j - 1) * 3) / 2.0;
185 m_SimulatedDerivativeT2Values(i+1,0) += m_FirstEPGProduct * m_SimulatedDerivativeT2Values(i,0) - m_SecondEPGProduct * m_SimulatedDerivativeT2Values(i,3);
186 if (m_NumberOfEchoes > 1)
187 m_SimulatedDerivativeT2Values(i+1,0) += m_ThirdEPGProduct * m_SimulatedDerivativeT2Values(i,5);
190 m_SimulatedDerivativeT2Values(i+1,1) += m_FourthEPGProduct * m_SimulatedDerivativeT2Values(i,2);
191 m_SimulatedDerivativeT2Values(i+1,2) += m_FirstEPGProduct * m_SimulatedDerivativeT2Values(i,1);
192 m_SimulatedDerivativeT2Values(i+1,3) += m_FifthEPGProduct * m_SimulatedDerivativeT2Values(i,3) - m_SecondEPGProduct * m_SimulatedDerivativeT2Values(i,0) / 2.0;
193 if (m_NumberOfEchoes > 1)
195 m_SimulatedDerivativeT2Values(i+1,2) -= m_SecondEPGProduct * m_SimulatedDerivativeT2Values(i,6);
196 m_SimulatedDerivativeT2Values(i+1,3) += m_SecondEPGProduct * m_SimulatedDerivativeT2Values(i,5) / 2.0;
198 if (m_NumberOfEchoes > 2)
199 m_SimulatedDerivativeT2Values(i+1,2) += m_ThirdEPGProduct * m_SimulatedDerivativeT2Values(i,8);
203 for (j = 1;j < m_NumberOfEchoes - 1;++j)
205 m_SimulatedDerivativeT2Values(i+1,1 + j * 3) += m_SecondEPGProduct * m_SimulatedDerivativeT2Values(i,3 + (j - 1) * 3) + m_FirstEPGProduct * m_SimulatedDerivativeT2Values(i,2 + j * 3);
207 m_SimulatedDerivativeT2Values(i+1,1 + j * 3) += m_ThirdEPGProduct * m_SimulatedDerivativeT2Values(i,1 + (j - 2) * 3);
209 m_SimulatedDerivativeT2Values(i+1,1 + j * 3) += m_ThirdEPGProduct * m_SimulatedDerivativeT2Values(i,0);
211 m_SimulatedDerivativeT2Values(i+1,2 + j * 3) += m_FirstEPGProduct * m_SimulatedDerivativeT2Values(i,1 + j * 3);
212 m_SimulatedDerivativeT2Values(i+1,3 + j * 3) += m_FifthEPGProduct * m_SimulatedDerivativeT2Values(i,3 + j * 3) - m_SecondEPGProduct * m_SimulatedDerivativeT2Values(i,1 + (j - 1) * 3) / 2.0;
214 if ((j + 1) < m_NumberOfEchoes)
216 m_SimulatedDerivativeT2Values(i+1,2 + j * 3) -= m_SecondEPGProduct * m_SimulatedDerivativeT2Values(i,3 + (j + 1) * 3);
217 m_SimulatedDerivativeT2Values(i+1,3 + j * 3) += m_SecondEPGProduct * m_SimulatedDerivativeT2Values(i,2 + (j + 1) * 3) / 2.0;
219 if ((j + 2) < m_NumberOfEchoes)
220 m_SimulatedDerivativeT2Values(i+1,2 + j * 3) += m_ThirdEPGProduct * m_SimulatedDerivativeT2Values(i,2 + (j + 2) * 3);
225 j = m_NumberOfEchoes - 1;
226 m_SimulatedDerivativeT2Values(i+1,1 + j * 3) += m_SecondEPGProduct * m_SimulatedDerivativeT2Values(i,3 + (j - 1) * 3) + m_FirstEPGProduct * m_SimulatedDerivativeT2Values(i,2 + j * 3);
227 m_SimulatedDerivativeT2Values(i+1,3 + j * 3) += m_FifthEPGProduct * m_SimulatedDerivativeT2Values(i,3 + j * 3) - m_SecondEPGProduct * m_SimulatedDerivativeT2Values(i,1 + (j - 1) * 3) / 2.0;
230 m_SimulatedDerivativeT2Values(i+1,1 + j * 3) += m_ThirdEPGProduct * m_SimulatedDerivativeT2Values(i,1 + (j - 2) * 3);
232 m_SimulatedDerivativeT2Values(i+1,1 + j * 3) += m_ThirdEPGProduct * m_SimulatedDerivativeT2Values(i,0);
234 m_OutputB1Derivative[i] = m_SimulatedDerivativeT2Values(i+1,0);
237 return m_OutputB1Derivative;
std::vector< double > RealVectorType
void ComputeT2SignalMatrixElements(double t1Value, double t2Value, double flipAngle)
RealVectorType & GetValue(double t1Value, double t2Value, double flipAngle, double m0Value)
Get EPG values at given point.
RealVectorType & GetFADerivative()
Get EPG derivative values at same point that was used for getting EPG values. Requires a run of GetVa...