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