MNE-CPP  beta 1.0
main.cpp
1 //=============================================================================================================
37 //*************************************************************************************************************
38 //=============================================================================================================
39 // INCLUDES
40 //=============================================================================================================
41 
42 #include <disp3DNew/brainview.h>
43 
44 #include <fs/label.h>
45 #include <fs/surface.h>
46 #include <fs/annotationset.h>
47 
48 #include <fiff/fiff_evoked.h>
49 #include <fiff/fiff.h>
50 #include <mne/mne.h>
51 
53 
54 #include <mne/mne_sourceestimate.h>
56 
57 #include <utils/mnemath.h>
58 
59 #include <iostream>
60 
61 //*************************************************************************************************************
62 //=============================================================================================================
63 // QT INCLUDES
64 //=============================================================================================================
65 
66 #include <QGuiApplication>
67 
68 
69 //*************************************************************************************************************
70 //=============================================================================================================
71 // USED NAMESPACES
72 //=============================================================================================================
73 
74 using namespace DISP3DNEWLIB;
75 using namespace MNELIB;
76 using namespace FSLIB;
77 using namespace FIFFLIB;
78 using namespace INVERSELIB;
79 using namespace UTILSLIB;
80 
81 
82 //*************************************************************************************************************
83 //=============================================================================================================
84 // MAIN
85 //=============================================================================================================
86 
87 //=============================================================================================================
96 int main(int argc, char *argv[])
97 {
98  QGuiApplication a(argc, argv);
99 
100  //*************************** Create source estimate *****************************************//
101  QFile t_fileFwd("./MNE-sample-data/MEG/sample/sample_audvis-meg-eeg-oct-6-fwd.fif");
102  QFile t_fileCov("./MNE-sample-data/MEG/sample/sample_audvis-cov.fif");
103  QFile t_fileRaw("./MNE-sample-data/MEG/sample/sample_audvis_raw.fif");
104  QString t_sEventName = "./MNE-sample-data/MEG/sample/sample_audvis_raw-eve.fif";
105 
106  qint32 event = 1;
107 
108  float tmin = -0.2f;
109  float tmax = 0.4f;
110 
111  bool keep_comp = false;
112  fiff_int_t dest_comp = 0;
113  bool pick_all = true;
114 
115  qint32 k, p;
116 
117  //
118  // Setup for reading the raw data
119  //
120  FiffRawData raw(t_fileRaw);
121 
122  RowVectorXi picks;
123  if (pick_all)
124  {
125  //
126  // Pick all
127  //
128  picks.resize(raw.info.nchan);
129 
130  for(k = 0; k < raw.info.nchan; ++k)
131  picks(k) = k;
132  //
133  }
134  else
135  {
136  QStringList include;
137  include << "STI 014";
138  bool want_meg = true;
139  bool want_eeg = false;
140  bool want_stim = false;
141 
142 // picks = Fiff::pick_types(raw.info, want_meg, want_eeg, want_stim, include, raw.info.bads);
143  picks = raw.info.pick_types(want_meg, want_eeg, want_stim, include, raw.info.bads);//prefer member function
144  }
145 
146  QStringList ch_names;
147  for(k = 0; k < picks.cols(); ++k)
148  ch_names << raw.info.ch_names[picks(0,k)];
149 
150  //
151  // Set up projection
152  //
153  if (raw.info.projs.size() == 0)
154  printf("No projector specified for these data\n");
155  else
156  {
157  //
158  // Activate the projection items
159  //
160  for (k = 0; k < raw.info.projs.size(); ++k)
161  raw.info.projs[k].active = true;
162 
163  printf("%d projection items activated\n",raw.info.projs.size());
164  //
165  // Create the projector
166  //
167 // fiff_int_t nproj = MNE::make_projector_info(raw.info, raw.proj); Using the member function instead
168  fiff_int_t nproj = raw.info.make_projector(raw.proj);
169 
170  if (nproj == 0)
171  {
172  printf("The projection vectors do not apply to these channels\n");
173  }
174  else
175  {
176  printf("Created an SSP operator (subspace dimension = %d)\n",nproj);
177  }
178  }
179 
180  //
181  // Set up the CTF compensator
182  //
183 // qint32 current_comp = MNE::get_current_comp(raw.info);
184  qint32 current_comp = raw.info.get_current_comp();
185  if (current_comp > 0)
186  printf("Current compensation grade : %d\n",current_comp);
187 
188  if (keep_comp)
189  dest_comp = current_comp;
190 
191  if (current_comp != dest_comp)
192  {
193  qDebug() << "This part needs to be debugged";
194  if(MNE::make_compensator(raw.info, current_comp, dest_comp, raw.comp))
195  {
196 // raw.info.chs = MNE::set_current_comp(raw.info.chs,dest_comp);
197  raw.info.set_current_comp(dest_comp);
198  printf("Appropriate compensator added to change to grade %d.\n",dest_comp);
199  }
200  else
201  {
202  printf("Could not make the compensator\n");
203  return 0;
204  }
205  }
206 
207  //
208  // Read the events
209  //
210  QFile t_EventFile;
211  MatrixXi events;
212  if (t_sEventName.size() == 0)
213  {
214  p = t_fileRaw.fileName().indexOf(".fif");
215  if (p > 0)
216  {
217  t_sEventName = t_fileRaw.fileName().replace(p, 4, "-eve.fif");
218  }
219  else
220  {
221  printf("Raw file name does not end properly\n");
222  return 0;
223  }
224 // events = mne_read_events(t_sEventName);
225 
226  t_EventFile.setFileName(t_sEventName);
227  MNE::read_events(t_EventFile, events);
228  printf("Events read from %s\n",t_sEventName.toUtf8().constData());
229  }
230  else
231  {
232  //
233  // Binary file
234  //
235  p = t_fileRaw.fileName().indexOf(".fif");
236  if (p > 0)
237  {
238  t_EventFile.setFileName(t_sEventName);
239  if(!MNE::read_events(t_EventFile, events))
240  {
241  printf("Error while read events.\n");
242  return 0;
243  }
244  printf("Binary event file %s read\n",t_sEventName.toUtf8().constData());
245  }
246  else
247  {
248  //
249  // Text file
250  //
251  printf("Text file %s is not supported jet.\n",t_sEventName.toUtf8().constData());
252  }
253  }
254 
255  //
256  // Select the desired events
257  //
258  qint32 count = 0;
259  MatrixXi selected = MatrixXi::Zero(1, events.rows());
260  for (p = 0; p < events.rows(); ++p)
261  {
262  if (events(p,1) == 0 && events(p,2) == event)
263  {
264  selected(0,count) = p;
265  ++count;
266  }
267  }
268  selected.conservativeResize(1, count);
269  if (count > 0)
270  printf("%d matching events found\n",count);
271  else
272  {
273  printf("No desired events found.\n");
274  return 0;
275  }
276 
277 
278  fiff_int_t event_samp, from, to;
279  MatrixXd timesDummy;
280 
281  MNEEpochDataList data;
282 
283  MNEEpochData* epoch = NULL;
284 
285  MatrixXd times;
286 
287  for (p = 0; p < count; ++p)
288  {
289  //
290  // Read a data segment
291  //
292  event_samp = events(selected(p),0);
293  from = event_samp + tmin*raw.info.sfreq;
294  to = event_samp + floor(tmax*raw.info.sfreq + 0.5);
295 
296  epoch = new MNEEpochData();
297 
298  if(raw.read_raw_segment(epoch->epoch, timesDummy, from, to, picks))
299  {
300  if (p == 0)
301  {
302  times.resize(1, to-from+1);
303  for (qint32 i = 0; i < times.cols(); ++i)
304  times(0, i) = ((float)(from-event_samp+i)) / raw.info.sfreq;
305  }
306 
307  epoch->event = event;
308  epoch->tmin = ((float)(from)-(float)(raw.first_samp))/raw.info.sfreq;
309  epoch->tmax = ((float)(to)-(float)(raw.first_samp))/raw.info.sfreq;
310 
311  data.append(MNEEpochData::SPtr(epoch));//List takes ownwership of the pointer - no delete need
312  }
313  else
314  {
315  printf("Can't read the event data segments");
316  return 0;
317  }
318  }
319 
320  if(data.size() > 0)
321  {
322  printf("Read %d epochs, %d samples each.\n",data.size(),(qint32)data[0]->epoch.cols());
323 
324  //DEBUG
325  std::cout << data[0]->epoch.block(0,0,10,10) << std::endl;
326  qDebug() << data[0]->epoch.rows() << " x " << data[0]->epoch.cols();
327 
328  std::cout << times.block(0,0,1,10) << std::endl;
329  qDebug() << times.rows() << " x " << times.cols();
330  }
331 
332  //
333  // calculate the average
334 
335  //Option 1
336  qint32 numAverages = 99;
337  VectorXi vecSel(numAverages);
338  srand (time(NULL)); // initialize random seed
339 
340  for(qint32 i = 0; i < vecSel.size(); ++i)
341  {
342  qint32 val = rand() % data.size();
343  vecSel(i) = val;
344  }
345 
346  std::cout << "Select following epochs to average:\n" << vecSel << std::endl;
347 
348  FiffEvoked evoked = data.average(raw.info, tmin*raw.info.sfreq, floor(tmax*raw.info.sfreq + 0.5), vecSel);
349 
350  //Source estimation
351  double snr = 1.0f;//0.1f;//1.0f;//3.0f;//0.1f;//3.0f;
352  QString method("dSPM"); //"MNE" | "dSPM" | "sLORETA"
353 
354  QString t_sFileNameClusteredInv("");
355  QString t_sFileNameStc("test_mind006_051209_auditory01.stc");
356 
357  // Parse command line parameters
358  for(qint32 i = 0; i < argc; ++i)
359  {
360  if(strcmp(argv[i], "-snr") == 0 || strcmp(argv[i], "--snr") == 0)
361  {
362  if(i + 1 < argc)
363  snr = atof(argv[i+1]);
364  }
365  else if(strcmp(argv[i], "-method") == 0 || strcmp(argv[i], "--method") == 0)
366  {
367  if(i + 1 < argc)
368  method = QString::fromUtf8(argv[i+1]);
369  }
370  else if(strcmp(argv[i], "-inv") == 0 || strcmp(argv[i], "--inv") == 0)
371  {
372  if(i + 1 < argc)
373  t_sFileNameClusteredInv = QString::fromUtf8(argv[i+1]);
374  }
375  else if(strcmp(argv[i], "-stc") == 0 || strcmp(argv[i], "--stc") == 0)
376  {
377  if(i + 1 < argc)
378  t_sFileNameStc = QString::fromUtf8(argv[i+1]);
379  }
380  }
381 
382  double lambda2 = 1.0 / pow(snr, 2);
383  qDebug() << "Start calculation with: SNR" << snr << "; Lambda" << lambda2 << "; Method" << method << "; stc:" << t_sFileNameStc;
384 
385  MNEForwardSolution t_Fwd(t_fileFwd);
386  if(t_Fwd.isEmpty())
387  return 1;
388 
389  FiffCov noise_cov(t_fileCov);
390 
391  // regularize noise covariance
392  noise_cov = noise_cov.regularize(evoked.info, 0.05, 0.05, 0.1, true);
393 
394  //
395  // make an inverse operators
396  //
397  FiffInfo info = evoked.info;
398 
399  MNEInverseOperator inverse_operator(info, t_Fwd, noise_cov, 0.2f, 0.8f);
400 
401  //
402  // Compute inverse solution
403  //
404  MinimumNorm minimumNorm(inverse_operator, lambda2, method);
405 
406  MNESourceEstimate sourceEstimate = minimumNorm.calculateInverse(evoked);
407 
408  if(sourceEstimate.isEmpty())
409  return 1;
410  // Create the test view
411  std::cout<<"Creating BrainView"<<std::endl;
412 
413  BrainView testWindow("sample", 2, "inflated", "./MNE-sample-data/subjects");
414 
415  testWindow.initStcDataModel("sample", 2, "inflated", "./MNE-sample-data/subjects", "aparc.a2009s", t_Fwd);
416 
417  testWindow.addSourceEstimate(sourceEstimate);
418 
419  testWindow.show();
420 
421  return a.exec();
422 }
static bool make_compensator(const FiffInfo &info, fiff_int_t from, fiff_int_t to, FiffCtfComp &ctf_comp, bool exclude_comp_chs=false)
Definition: mne.h:225
FIFF measurement file information.
Definition: fiff_info.h:96
MNEEpochDataList class declaration.
FreeSurfer surface visualisation.
Definition: brainview.h:129
evoked data
Definition: fiff_evoked.h:91
static bool read_events(QIODevice &p_IODevice, MatrixXi &eventlist)
Definition: mne.cpp:68
FIFF raw measurement data.
Definition: fiff_raw_data.h:94
Minimum norm class declaration.
covariance data
Definition: fiff_cov.h:94
Definition: fiff.h:98
Surface class declaration.
MNESourceEstimate class declaration.
FiffEvoked average(FiffInfo &p_info, fiff_int_t first, fiff_int_t last, VectorXi sel=defaultVectorXi, bool proj=false)
QSharedPointer< MNEEpochData > SPtr
Minimum norm estimation.
Definition: minimumnorm.h:82
AnnotationSet class declaration.
Label class declaration.
MNEMath class declaration.
FIFF class declaration, which provides static wrapper functions to stay consistent with mne matlab to...