ANIMA  4.0
animaEPGSignalSimulator.cxx
Go to the documentation of this file.
1 #include <cmath>
2 
3 #include <iostream>
4 #include "animaEPGSignalSimulator.h"
5 
6 namespace anima
7 {
8 
9 EPGSignalSimulator::EPGSignalSimulator()
10 {
11  m_NumberOfEchoes = 1;
12  m_EchoSpacing = 10;
13  m_ExcitationFlipAngle = M_PI / 2.0;
14 }
15 
16 EPGSignalSimulator::RealVectorType &EPGSignalSimulator::GetValue(double t1Value, double t2Value,
17  double flipAngle, double m0Value)
18 {
19  m_SimulatedT2Values.set_size(3 * m_NumberOfEchoes + 1,m_NumberOfEchoes + 1);
20  m_OutputVector.resize(m_NumberOfEchoes);
21 
22  this->ComputeT2SignalMatrixElements(t1Value,t2Value,flipAngle);
23 
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;
28 
29  // Loop on all signals to be generated
30  for (unsigned int i = 0;i < m_NumberOfEchoes;++i)
31  {
32  // First line
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);
36 
37  // First block
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;
41 
42  if (m_NumberOfEchoes > 1)
43  {
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;
46 
47  if (m_NumberOfEchoes > 2)
48  m_SimulatedT2Values(i+1,2) += m_ThirdEPGProduct * m_SimulatedT2Values(i,8);
49  }
50 
51  //center blocks
52  unsigned int j;
53  for (j = 1;j < m_NumberOfEchoes - 1;++j)
54  {
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);
56  if (j > 1)
57  m_SimulatedT2Values(i+1,1 + j * 3) += m_ThirdEPGProduct * m_SimulatedT2Values(i,1 + (j - 2) * 3);
58  else
59  m_SimulatedT2Values(i+1,1 + j * 3) += m_ThirdEPGProduct * m_SimulatedT2Values(i,0);
60 
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;
63 
64  if ((j + 1) < m_NumberOfEchoes)
65  {
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;
68 
69  if ((j + 2) < m_NumberOfEchoes)
70  m_SimulatedT2Values(i+1,2 + j * 3) += m_ThirdEPGProduct * m_SimulatedT2Values(i,2 + (j + 2) * 3);
71  }
72  }
73 
74  // end block line
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;
79 
80  if (j > 1)
81  m_SimulatedT2Values(i+1,1 + j * 3) += m_ThirdEPGProduct * m_SimulatedT2Values(i,1 + (j - 2) * 3);
82  else
83  m_SimulatedT2Values(i+1,1 + j * 3) += m_ThirdEPGProduct * m_SimulatedT2Values(i,0);
84 
85  m_OutputVector[i] = m_SimulatedT2Values(i+1,0);
86  }
87 
88  return m_OutputVector;
89 }
90 
91 void EPGSignalSimulator::ComputeT2SignalMatrixElements(double t1Value, double t2Value,
92  double flipAngle)
93 {
94  double espT2Value = 0.0;
95  if (t2Value != 0.0)
96  espT2Value = std::exp(- m_EchoSpacing / (2 * t2Value));
97 
98  double espT1Value = 0.0;
99  if (t1Value != 0.0)
100  espT1Value = std::exp(- m_EchoSpacing / (2 * t1Value));
101 
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);
106 
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;
112 
113  m_FirstDerivativeProduct = cosB1alpha2 * sinB1alpha2 * espT2Value * espT2Value;
114  m_SecondDerivativeProduct = cosB1alpha * espT1Value * espT2Value;
115  m_ThirdDerivativeProduct = - sinB1alpha * espT1Value * espT1Value;
116 }
117 
118 EPGSignalSimulator::RealVectorType &EPGSignalSimulator::GetFADerivative()
119 {
120  m_OutputB1Derivative.resize(m_NumberOfEchoes);
121 
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;
125 
126  // Loop on all signals to be generated
127  for (unsigned int i = 0;i < m_NumberOfEchoes;++i)
128  {
129  // Start by doing dE * simulatedValues
130  // First line
131  m_SimulatedDerivativeT2Values(i+1,0) = m_FirstDerivativeProduct * m_SimulatedT2Values(i,0) - m_SecondDerivativeProduct * m_SimulatedT2Values(i,3);
132 
133  // First block
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;
137 
138  if (m_NumberOfEchoes > 1)
139  {
140  m_SimulatedDerivativeT2Values(i+1,0) -= m_FirstDerivativeProduct * m_SimulatedT2Values(i,5);
141  m_SimulatedDerivativeT2Values(i+1,2) -= m_SecondDerivativeProduct * m_SimulatedT2Values(i,6);
142 
143  if (m_NumberOfEchoes > 2)
144  m_SimulatedDerivativeT2Values(i+1,2) -= m_FirstDerivativeProduct * m_SimulatedT2Values(i,8);
145 
146  m_SimulatedDerivativeT2Values(i+1,3) += m_SecondDerivativeProduct * m_SimulatedT2Values(i,5) / 2.0;
147  }
148 
149  //center blocks
150  unsigned int j;
151  for (j = 1;j < m_NumberOfEchoes - 1;++j)
152  {
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);
154  if (j > 1)
155  m_SimulatedDerivativeT2Values(i+1,1 + j * 3) -= m_FirstDerivativeProduct * m_SimulatedT2Values(i,1 + (j - 2) * 3);
156  else
157  m_SimulatedDerivativeT2Values(i+1,1 + j * 3) -= m_FirstDerivativeProduct * m_SimulatedT2Values(i,0);
158 
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;
161 
162  if ((j + 1) < m_NumberOfEchoes)
163  {
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;
166 
167  if ((j + 2) < m_NumberOfEchoes)
168  m_SimulatedDerivativeT2Values(i+1,2 + j * 3) -= m_FirstDerivativeProduct * m_SimulatedT2Values(i,2 + (j + 2) * 3);
169  }
170  }
171 
172  // end block line
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);
175  if (j > 1)
176  m_SimulatedDerivativeT2Values(i+1,1 + j * 3) -= m_FirstDerivativeProduct * m_SimulatedT2Values(i,1 + (j - 2) * 3);
177  else
178  m_SimulatedDerivativeT2Values(i+1,1 + j * 3) -= m_FirstDerivativeProduct * m_SimulatedT2Values(i,0);
179 
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;
182 
183  // Now adding E * simulatedDerivative[i]
184  // First line
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);
188 
189  // First block
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)
194  {
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;
197 
198  if (m_NumberOfEchoes > 2)
199  m_SimulatedDerivativeT2Values(i+1,2) += m_ThirdEPGProduct * m_SimulatedDerivativeT2Values(i,8);
200  }
201 
202  //center blocks
203  for (j = 1;j < m_NumberOfEchoes - 1;++j)
204  {
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);
206  if (j > 1)
207  m_SimulatedDerivativeT2Values(i+1,1 + j * 3) += m_ThirdEPGProduct * m_SimulatedDerivativeT2Values(i,1 + (j - 2) * 3);
208  else
209  m_SimulatedDerivativeT2Values(i+1,1 + j * 3) += m_ThirdEPGProduct * m_SimulatedDerivativeT2Values(i,0);
210 
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;
213 
214  if ((j + 1) < m_NumberOfEchoes)
215  {
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;
218 
219  if ((j + 2) < m_NumberOfEchoes)
220  m_SimulatedDerivativeT2Values(i+1,2 + j * 3) += m_ThirdEPGProduct * m_SimulatedDerivativeT2Values(i,2 + (j + 2) * 3);
221  }
222  }
223 
224  // end block line
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;
228 
229  if (j > 1)
230  m_SimulatedDerivativeT2Values(i+1,1 + j * 3) += m_ThirdEPGProduct * m_SimulatedDerivativeT2Values(i,1 + (j - 2) * 3);
231  else
232  m_SimulatedDerivativeT2Values(i+1,1 + j * 3) += m_ThirdEPGProduct * m_SimulatedDerivativeT2Values(i,0);
233 
234  m_OutputB1Derivative[i] = m_SimulatedDerivativeT2Values(i+1,0);
235  }
236 
237  return m_OutputB1Derivative;
238 }
239 
240 } // end of namespace anima
anima::EPGSignalSimulator::RealVectorType
std::vector< double > RealVectorType
Definition: animaEPGSignalSimulator.h:17
anima::EPGSignalSimulator::ComputeT2SignalMatrixElements
void ComputeT2SignalMatrixElements(double t1Value, double t2Value, double flipAngle)
Definition: animaEPGSignalSimulator.cxx:91
anima::EPGSignalSimulator::GetValue
RealVectorType & GetValue(double t1Value, double t2Value, double flipAngle, double m0Value)
Get EPG values at given point.
Definition: animaEPGSignalSimulator.cxx:16
anima::EPGSignalSimulator::GetFADerivative
RealVectorType & GetFADerivative()
Get EPG derivative values at same point that was used for getting EPG values. Requires a run of GetVa...
Definition: animaEPGSignalSimulator.cxx:118
Generated by   doxygen 1.8.13