MNE-CPP  beta 1.0
minimumnorm.cpp
Go to the documentation of this file.
1 //=============================================================================================================
36 //*************************************************************************************************************
37 //=============================================================================================================
38 // INCLUDES
39 //=============================================================================================================
40 
41 #include "minimumnorm.h"
42 
43 #include <mne/mne_sourceestimate.h>
44 #include <fiff/fiff_evoked.h>
45 
46 
47 //*************************************************************************************************************
48 //=============================================================================================================
49 // Eigen INCLUDES
50 //=============================================================================================================
51 
52 #include <Eigen/Core>
53 
54 
55 //*************************************************************************************************************
56 //=============================================================================================================
57 // USED NAMESPACES
58 //=============================================================================================================
59 
60 using namespace Eigen;
61 using namespace MNELIB;
62 using namespace INVERSELIB;
63 
64 
65 //*************************************************************************************************************
66 //=============================================================================================================
67 // DEFINE MEMBER METHODS
68 //=============================================================================================================
69 
70 MinimumNorm::MinimumNorm(const MNEInverseOperator &p_inverseOperator, float lambda, const QString method)
71 : m_inverseOperator(p_inverseOperator)
72 , inverseSetup(false)
73 {
74  this->setRegularization(lambda);
75  this->setMethod(method);
76 }
77 
78 //*************************************************************************************************************
79 
80 MinimumNorm::MinimumNorm(const MNEInverseOperator &p_inverseOperator, float lambda, bool dSPM, bool sLORETA)
81 : m_inverseOperator(p_inverseOperator)
82 , inverseSetup(false)
83 {
84  this->setRegularization(lambda);
85  this->setMethod(dSPM, sLORETA);
86 }
87 
88 
89 //*************************************************************************************************************
90 
91 MNESourceEstimate MinimumNorm::calculateInverse(const FiffEvoked &p_fiffEvoked, bool pick_normal)
92 {
93  //
94  // Set up the inverse according to the parameters
95  //
96  qint32 nave = p_fiffEvoked.nave;
97 
98  if(!m_inverseOperator.check_ch_names(p_fiffEvoked.info))
99  {
100  qWarning("Channel name check failed.");
101  return MNESourceEstimate();
102  }
103 
104  doInverseSetup(nave,pick_normal);
105 
106  //
107  // Pick the correct channels from the data
108  //
109  FiffEvoked t_fiffEvoked = p_fiffEvoked.pick_channels(inv.noise_cov->names);
110 
111  printf("Picked %d channels from the data\n",t_fiffEvoked.info.nchan);
112 
113  //Results
114  float tmin = ((float)t_fiffEvoked.first) / t_fiffEvoked.info.sfreq;
115  float tstep = 1/t_fiffEvoked.info.sfreq;
116 
117  return calculateInverse(t_fiffEvoked.data, tmin, tstep);
118 
119 // //
120 // // Set up the inverse according to the parameters
121 // //
122 // qint32 nave = p_fiffEvoked.nave;
123 
124 // if(!m_inverseOperator.check_ch_names(p_fiffEvoked.info))
125 // {
126 // qWarning("Channel name check failed.");
127 // return SourceEstimate();
128 // }
129 
130 // //ToDo his could be heavily accelerated for real time calculation -> ToDo calculate inverse RT
131 // MNEInverseOperator inv = m_inverseOperator.prepare_inverse_operator(nave, m_fLambda, m_bdSPM, m_bsLORETA);
132 // //
133 // // Pick the correct channels from the data
134 // //
135 // FiffEvoked t_fiffEvoked = p_fiffEvoked.pick_channels(inv.noise_cov->names);
136 
137 // printf("Picked %d channels from the data\n",t_fiffEvoked.info.nchan);
138 // printf("Computing inverse...");
139 
140 // MatrixXd K;
141 // SparseMatrix<double> noise_norm;
142 // QList<VectorXi> vertno;
143 // Label label;
144 // inv.assemble_kernel(label, m_sMethod, pick_normal, K, noise_norm, vertno);
145 
146 // MatrixXd sol = K * t_fiffEvoked.data; //apply imaging kernel
147 
148 // if (inv.source_ori == FIFFV_MNE_FREE_ORI)
149 // {
150 // printf("combining the current components...");
151 // MatrixXd sol1(sol.rows()/3,sol.cols());
152 // for(qint32 i = 0; i < sol.cols(); ++i)
153 // {
154 // VectorXd* tmp = MNEMath::combine_xyz(sol.block(0,i,sol.rows(),1));
155 // sol1.block(0,i,sol.rows()/3,1) = tmp->cwiseSqrt();
156 // delete tmp;
157 // }
158 // sol.resize(sol1.rows(),sol1.cols());
159 // sol = sol1;
160 // }
161 
162 // if (m_bdSPM)
163 // {
164 // printf("(dSPM)...");
165 // sol = inv.noisenorm*sol;
166 // }
167 // else if (m_bsLORETA)
168 // {
169 // printf("(sLORETA)...");
170 // sol = inv.noisenorm*sol;
171 // }
172 // printf("[done]\n");
173 
174 // //Results
175 // float tmin = ((float)t_fiffEvoked.first) / t_fiffEvoked.info.sfreq;
176 // float tstep = 1/t_fiffEvoked.info.sfreq;
177 
178 // QList<VectorXi> t_qListVertices;
179 // for(qint32 h = 0; h < inv.src.size(); ++h)
180 // t_qListVertices.push_back(inv.src[h].vertno);
181 
182 // return SourceEstimate(sol, t_qListVertices, tmin, tstep);
183 }
184 
185 
186 //*************************************************************************************************************
187 
188 MNESourceEstimate MinimumNorm::calculateInverse(const MatrixXd &data, float tmin, float tstep) const
189 {
190  if(!inverseSetup)
191  {
192  qWarning("Inverse not setup -> call doInverseSetup first!");
193  return MNESourceEstimate();
194  }
195 
196  MatrixXd sol = K * data; //apply imaging kernel
197 
198  if (inv.source_ori == FIFFV_MNE_FREE_ORI)
199  {
200  printf("combining the current components...");
201  MatrixXd sol1(sol.rows()/3,sol.cols());
202  for(qint32 i = 0; i < sol.cols(); ++i)
203  {
204  VectorXd* tmp = MNEMath::combine_xyz(sol.block(0,i,sol.rows(),1));
205  sol1.block(0,i,sol.rows()/3,1) = tmp->cwiseSqrt();
206  delete tmp;
207  }
208  sol.resize(sol1.rows(),sol1.cols());
209  sol = sol1;
210  }
211 
212  if (m_bdSPM)
213  {
214  printf("(dSPM)...");
215  sol = inv.noisenorm*sol;
216  }
217  else if (m_bsLORETA)
218  {
219  printf("(sLORETA)...");
220  sol = inv.noisenorm*sol;
221  }
222  printf("[done]\n");
223 
224  //Results
225  VectorXi p_vecVertices(inv.src[0].vertno.size() + inv.src[1].vertno.size());
226  p_vecVertices << inv.src[0].vertno, inv.src[1].vertno;
227 
228 // VectorXi p_vecVertices();
229 // for(qint32 h = 0; h < inv.src.size(); ++h)
230 // t_qListVertices.push_back(inv.src[h].vertno);
231 
232  return MNESourceEstimate(sol, p_vecVertices, tmin, tstep);
233 
234 }
235 
236 
237 //*************************************************************************************************************
238 
239 void MinimumNorm::doInverseSetup(qint32 nave, bool pick_normal)
240 {
241  //
242  // Set up the inverse according to the parameters
243  //
244  inv = m_inverseOperator.prepare_inverse_operator(nave, m_fLambda, m_bdSPM, m_bsLORETA);
245 
246  printf("Computing inverse...");
247  inv.assemble_kernel(label, m_sMethod, pick_normal, K, noise_norm, vertno);
248 
249  std::cout << "K " << K.rows() << " x " << K.cols() << std::endl;
250 
251  inverseSetup = true;
252 }
253 
254 
255 //*************************************************************************************************************
256 
257 const char* MinimumNorm::getName() const
258 {
259  return "Minimum Norm Estimate";
260 }
261 
262 
263 //*************************************************************************************************************
264 
266 {
267  return m_inverseOperator.src;
268 }
269 
270 //*************************************************************************************************************
271 
272 void MinimumNorm::setMethod(QString method)
273 {
274  if(method.compare("MNE") == 0)
275  setMethod(false, false);
276  else if(method.compare("dSPM") == 0)
277  setMethod(true, false);
278  else if(method.compare("sLORETA") == 0)
279  setMethod(false, true);
280  else
281  {
282  qWarning("Method not recognized!");
283  method = "dSPM";
284  setMethod(true, false);
285  }
286 
287  printf("\tSet minimum norm method to %s.\n", method.toLatin1().constData());
288 }
289 
290 
291 //*************************************************************************************************************
292 
293 void MinimumNorm::setMethod(bool dSPM, bool sLORETA)
294 {
295  if(dSPM && sLORETA)
296  {
297  qWarning("Cant activate dSPM and sLORETA at the same time! - Activating dSPM");
298  m_bdSPM = true;
299  m_bsLORETA = false;
300  }
301  else
302  {
303  m_bdSPM = dSPM;
304  m_bsLORETA = sLORETA;
305  if(dSPM)
306  m_sMethod = QString("dSPM");
307  else if(sLORETA)
308  m_sMethod = QString("sLORETA");
309  else
310  m_sMethod = QString("MNE");
311 
312  }
313 }
314 
315 
316 //*************************************************************************************************************
317 
319 {
320  m_fLambda = lambda;
321 }
virtual const char * getName() const
Source Space descritpion.
bool assemble_kernel(const Label &label, QString method, bool pick_normal, MatrixXd &K, SparseMatrix< double > &noise_norm, QList< VectorXi > &vertno)
void setMethod(QString method)
evoked data
Definition: fiff_evoked.h:91
static VectorXd * combine_xyz(const VectorXd &vec)
Definition: mnemath.cpp:79
SparseMatrix< double > noisenorm
Minimum norm class declaration.
MNEInverseOperator prepare_inverse_operator(qint32 nave, float lambda2, bool dSPM, bool sLORETA=false) const
bool check_ch_names(const FiffInfo &info) const
FiffEvoked pick_channels(const QStringList &include=defaultQStringList, const QStringList &exclude=defaultQStringList) const
MNESourceEstimate class declaration.
void setRegularization(float lambda)
virtual const MNESourceSpace & getSourceSpace() const
MinimumNorm(const MNEInverseOperator &p_inverseOperator, float lambda, const QString method)
Definition: minimumnorm.cpp:70
virtual MNESourceEstimate calculateInverse(const FiffEvoked &p_fiffEvoked, bool pick_normal=false)
Definition: minimumnorm.cpp:91