MNE-CPP  beta 1.0
main.cpp
1 //=============================================================================================================
37 //*************************************************************************************************************
38 //=============================================================================================================
39 // INCLUDES
40 //=============================================================================================================
41 
42 #include <fs/label.h>
43 #include <fs/surface.h>
44 #include <fs/annotationset.h>
45 
46 #include <fiff/fiff_evoked.h>
47 #include <fiff/fiff.h>
48 #include <mne/mne.h>
49 
51 
52 #include <mne/mne_sourceestimate.h>
54 
55 #include <disp3D/inverseview.h>
56 
57 #include <utils/mnemath.h>
58 
59 #include <iostream>
60 
61 
62 //*************************************************************************************************************
63 //=============================================================================================================
64 // QT INCLUDES
65 //=============================================================================================================
66 
67 #include <QGuiApplication>
68 #include <QSet>
69 #include <QElapsedTimer>
70 
71 //#define BENCHMARK
72 
73 
74 //*************************************************************************************************************
75 //=============================================================================================================
76 // USED NAMESPACES
77 //=============================================================================================================
78 
79 using namespace MNELIB;
80 using namespace FSLIB;
81 using namespace FIFFLIB;
82 using namespace INVERSELIB;
83 using namespace DISP3DLIB;
84 using namespace UTILSLIB;
85 
86 
87 //*************************************************************************************************************
88 //=============================================================================================================
89 // MAIN
90 //=============================================================================================================
91 
92 //=============================================================================================================
101 int main(int argc, char *argv[])
102 {
103  QGuiApplication a(argc, argv);
104 
105 // QFile t_fileRaw("./MNE-sample-data/MEG/sample/sample_audvis_raw.fif");
106 // QString t_sEventName = "./MNE-sample-data/MEG/sample/sample_audvis_raw-eve.fif";
107 // QFile t_fileFwd("./MNE-sample-data/MEG/sample/sample_audvis-meg-eeg-oct-6-fwd.fif");
108 // AnnotationSet t_annotationSet("sample", 2, "aparc.a2009s", "./MNE-sample-data/subjects");
109 // SurfaceSet t_surfSet("sample", 2, "white", "./MNE-sample-data/subjects");
110 
111 
112 // QFile t_fileRaw("D:/Users/Christoph/SkyDrive/Thesis_Data/MIND/mind004_050924_median01_raw.fif");
113 // QString t_sEventName = "D:/Users/Christoph/SkyDrive/Thesis_Data/MIND/mind004_050924_median01_raw-eve.fif";
114 // QFile t_fileFwd("D:/Users/Christoph/SkyDrive/Thesis_Data/MIND/mind004_050924_median01_raw-oct-6-fwd.fif");
115 // AnnotationSet t_annotationSet("mind004", 2, "aparc.a2009s", "D:/Users/Christoph/SkyDrive/Thesis_Data/subjects");
116 // SurfaceSet t_surfSet("mind004", 2, "white", "D:/Users/Christoph/SkyDrive/Thesis_Data/subjects");
117 
118 
119  QFile t_fileRaw("D:/Users/Christoph/SkyDrive/Thesis_Data/MIND/mind006_051209_auditory01_raw.fif");
120  QString t_sEventName = "D:/Users/Christoph/SkyDrive/Thesis_Data/MIND/mind006_051209_auditory01_raw-eve.fif";
121  QFile t_fileFwd("D:/Users/Christoph/SkyDrive/Thesis_Data/MIND/mind006_051209_auditory01_raw-oct-6-fwd.fif");
122  AnnotationSet t_annotationSet("mind006", 2, "aparc.a2009s", "D:/Users/Christoph/SkyDrive/Thesis_Data/subjects");
123  SurfaceSet t_surfSet("mind006", 2, "white", "D:/Users/Christoph/SkyDrive/Thesis_Data/subjects");
124 
125 // QFile t_fileRaw("D:/Users/Christoph/SkyDrive/Thesis_Data/MIND/mind006_051210_median02_raw.fif");
126 // QString t_sEventName = "D:/Users/Christoph/SkyDrive/Thesis_Data/MIND/mind006_051210_median02_raw-eve.fif";
127 // QFile t_fileFwd("D:/Users/Christoph/SkyDrive/Thesis_Data/MIND/mind006_051210_median02_raw-oct-6-fwd.fif");
128 // AnnotationSet t_annotationSet("mind006", 2, "aparc.a2009s", "D:/Users/Christoph/SkyDrive/Thesis_Data/subjects");
129 // SurfaceSet t_surfSet("mind006", 2, "white", "D:/Users/Christoph/SkyDrive/Thesis_Data/subjects");
130 
131  QString t_sFileNameStc("");//("mind006_051209_auditory01.stc");
132 
133 
134  bool doMovie = false;
135 
136  qint32 numDipolePairs = 1;//7;
137 
138 // //Right Medianus MIND004
139 // qint32 event = 65;
140 // float tmin = 0.0f;
141 // float tmax = 0.1f;
142 
143  //Right500 Auditory MIND006
144  qint32 event = 2;
145 
146  float tmin = 0.1f;
147  float tmax = 0.2f;
148 
149  bool keep_comp = false;
150  fiff_int_t dest_comp = 0;
151  bool pick_all = true;
152 
153  qint32 k, p;
154 
155 
156  // Parse command line parameters
157  for(qint32 i = 0; i < argc; ++i)
158  {
159  if(strcmp(argv[i], "-stc") == 0 || strcmp(argv[i], "--stc") == 0)
160  {
161  if(i + 1 < argc)
162  t_sFileNameStc = QString::fromUtf8(argv[i+1]);
163  }
164  }
165 
166  //
167  // Load data
168  //
169  MNEForwardSolution t_Fwd(t_fileFwd);
170  if(t_Fwd.isEmpty())
171  return 1;
172 
173  //
174  // Setup for reading the raw data
175  //
176  FiffRawData raw(t_fileRaw);
177 
178  RowVectorXi picks;
179  if (pick_all)
180  {
181  //
182  // Pick all
183  //
184  picks.resize(raw.info.nchan);
185 
186  for(k = 0; k < raw.info.nchan; ++k)
187  picks(k) = k;
188  //
189  }
190  else
191  {
192  QStringList include;
193  include << "STI 014";
194  bool want_meg = true;
195  bool want_eeg = false;
196  bool want_stim = false;
197 
198  picks = raw.info.pick_types(want_meg, want_eeg, want_stim, include, raw.info.bads);//prefer member function
199  }
200 
201  QStringList ch_names;
202  for(k = 0; k < picks.cols(); ++k)
203  ch_names << raw.info.ch_names[picks(0,k)];
204 
205  //
206  // Set up projection
207  //
208  if (raw.info.projs.size() == 0)
209  printf("No projector specified for these data\n");
210  else
211  {
212  //
213  // Activate the projection items
214  //
215  for (k = 0; k < raw.info.projs.size(); ++k)
216  raw.info.projs[k].active = true;
217 
218  printf("%d projection items activated\n",raw.info.projs.size());
219  //
220  // Create the projector
221  //
222 // fiff_int_t nproj = MNE::make_projector_info(raw.info, raw.proj); Using the member function instead
223  fiff_int_t nproj = raw.info.make_projector(raw.proj);
224 
225  if (nproj == 0)
226  {
227  printf("The projection vectors do not apply to these channels\n");
228  }
229  else
230  {
231  printf("Created an SSP operator (subspace dimension = %d)\n",nproj);
232  }
233  }
234 
235  //
236  // Set up the CTF compensator
237  //
238  qint32 current_comp = raw.info.get_current_comp();
239  if (current_comp > 0)
240  printf("Current compensation grade : %d\n",current_comp);
241 
242  if (keep_comp)
243  dest_comp = current_comp;
244 
245  if (current_comp != dest_comp)
246  {
247  qDebug() << "This part needs to be debugged";
248  if(MNE::make_compensator(raw.info, current_comp, dest_comp, raw.comp))
249  {
250  raw.info.set_current_comp(dest_comp);
251  printf("Appropriate compensator added to change to grade %d.\n",dest_comp);
252  }
253  else
254  {
255  printf("Could not make the compensator\n");
256  return 0;
257  }
258  }
259  //
260  // Read the events
261  //
262  QFile t_EventFile;
263  MatrixXi events;
264  if (t_sEventName.size() == 0)
265  {
266  p = t_fileRaw.fileName().indexOf(".fif");
267  if (p > 0)
268  {
269  t_sEventName = t_fileRaw.fileName().replace(p, 4, "-eve.fif");
270  }
271  else
272  {
273  printf("Raw file name does not end properly\n");
274  return 0;
275  }
276 // events = mne_read_events(t_sEventName);
277 
278  t_EventFile.setFileName(t_sEventName);
279  MNE::read_events(t_EventFile, events);
280  printf("Events read from %s\n",t_sEventName.toUtf8().constData());
281  }
282  else
283  {
284  //
285  // Binary file
286  //
287  p = t_fileRaw.fileName().indexOf(".fif");
288  if (p > 0)
289  {
290  t_EventFile.setFileName(t_sEventName);
291  if(!MNE::read_events(t_EventFile, events))
292  {
293  printf("Error while read events.\n");
294  return 0;
295  }
296  printf("Binary event file %s read\n",t_sEventName.toUtf8().constData());
297  }
298  else
299  {
300  //
301  // Text file
302  //
303  printf("Text file %s is not supported jet.\n",t_sEventName.toUtf8().constData());
304 // try
305 // events = load(eventname);
306 // catch
307 // error(me,mne_omit_first_line(lasterr));
308 // end
309 // if size(events,1) < 1
310 // error(me,'No data in the event file');
311 // end
312 // //
313 // // Convert time to samples if sample number is negative
314 // //
315 // for p = 1:size(events,1)
316 // if events(p,1) < 0
317 // events(p,1) = events(p,2)*raw.info.sfreq;
318 // end
319 // end
320 // //
321 // // Select the columns of interest (convert to integers)
322 // //
323 // events = int32(events(:,[1 3 4]));
324 // //
325 // // New format?
326 // //
327 // if events(1,2) == 0 && events(1,3) == 0
328 // fprintf(1,'The text event file %s is in the new format\n',eventname);
329 // if events(1,1) ~= raw.first_samp
330 // error(me,'This new format event file is not compatible with the raw data');
331 // end
332 // else
333 // fprintf(1,'The text event file %s is in the old format\n',eventname);
334 // //
335 // // Offset with first sample
336 // //
337 // events(:,1) = events(:,1) + raw.first_samp;
338 // end
339  }
340  }
341 
342  //
343  // Select the desired events
344  //
345  qint32 count = 0;
346  MatrixXi selected = MatrixXi::Zero(1, events.rows());
347  for (p = 0; p < events.rows(); ++p)
348  {
349  if (events(p,1) == 0 && events(p,2) == event)
350  {
351  selected(0,count) = p;
352  ++count;
353  }
354  }
355  selected.conservativeResize(1, count);
356  if (count > 0)
357  printf("%d matching events found\n",count);
358  else
359  {
360  printf("No desired events found.\n");
361  return 0;
362  }
363 
364 
365  fiff_int_t event_samp, from, to;
366  MatrixXd timesDummy;
367 
368  MNEEpochDataList data;
369 
370  MNEEpochData* epoch = NULL;
371 
372  MatrixXd times;
373 
374  for (p = 0; p < count; ++p)
375  {
376  //
377  // Read a data segment
378  //
379  event_samp = events(selected(p),0);
380  from = event_samp + tmin*raw.info.sfreq;
381  to = event_samp + floor(tmax*raw.info.sfreq + 0.5);
382 
383  epoch = new MNEEpochData();
384 
385  if(raw.read_raw_segment(epoch->epoch, timesDummy, from, to, picks))
386  {
387  if (p == 0)
388  {
389  times.resize(1, to-from+1);
390  for (qint32 i = 0; i < times.cols(); ++i)
391  times(0, i) = ((float)(from-event_samp+i)) / raw.info.sfreq;
392  }
393 
394  epoch->event = event;
395  epoch->tmin = ((float)(from)-(float)(raw.first_samp))/raw.info.sfreq;
396  epoch->tmax = ((float)(to)-(float)(raw.first_samp))/raw.info.sfreq;
397 
398  data.append(MNEEpochData::SPtr(epoch));//List takes ownwership of the pointer - no delete need
399  }
400  else
401  {
402  printf("Can't read the event data segments");
403  return 0;
404  }
405  }
406 
407  if(data.size() > 0)
408  {
409  printf("Read %d epochs, %d samples each.\n",data.size(),(qint32)data[0]->epoch.cols());
410 
411 // //DEBUG
412 // std::cout << data[0]->epoch.block(0,0,10,10) << std::endl;
413 // qDebug() << data[0]->epoch.rows() << " x " << data[0]->epoch.cols();
414 
415 // std::cout << times.block(0,0,1,10) << std::endl;
416 // qDebug() << times.rows() << " x " << times.cols();
417  }
418 
419  //
420  // calculate the average
421  //
422 // //Option 1
423 // qint32 numAverages = 99;
424 // VectorXi vecSel(numAverages);
425 // srand (time(NULL)); // initialize random seed
426 
427 // for(qint32 i = 0; i < vecSel.size(); ++i)
428 // {
429 // qint32 val = rand() % data.size();
430 // vecSel(i) = val;
431 // }
432 
433  //Option 2
434 // VectorXi vecSel(20);
435 
437 
438 // vecSel << 65, 22, 47, 55, 16, 29, 14, 36, 57, 97, 89, 46, 9, 93, 83, 52, 71, 52, 3, 96;
439 
440 // //Option 3 Newest
441 // VectorXi vecSel(10);
442 
443 // vecSel << 0, 96, 80, 55, 66, 25, 26, 2, 55, 58, 6, 88;
444 
445 
446 // VectorXi vecSel(1);
447 
448 // vecSel << 0;//2,3;
449 
450 
451  VectorXi vecSel(2);//153, 147 or 146, 113
452 // //MIND 004 medianus 01
453 // vecSel << 147, 153;
454 
455  //MIND 006 auditory 01
456 // vecSel << 42, 48; //1
457 // vecSel << 77, 113; //Perfect!
458 // vecSel << 48, 85; //3
459 // vecSel << 56, 83; //3
460 
461 
462  srand (time(NULL)); // initialize random seed
463 
464  for(qint32 i = 0; i < vecSel.size(); ++i)
465  {
466  qint32 val = rand() % count;
467  vecSel(i) = val;
468  }
469 
470 
471  std::cout << "Select following epochs to average:\n" << vecSel << std::endl;
472 
473  FiffEvoked evoked = data.average(raw.info, tmin*raw.info.sfreq, floor(tmax*raw.info.sfreq + 0.5), vecSel);
474 
475  QStringList ch_sel_names = t_Fwd.info.ch_names;
476  FiffEvoked pickedEvoked = evoked.pick_channels(ch_sel_names);
477 
478 
479 
480  //########################################################################################
481  // RAP MUSIC Source Estimate
482 
483  //
484  // Cluster forward solution;
485  //
486  MNEForwardSolution t_clusteredFwd = t_Fwd.cluster_forward_solution(t_annotationSet, 20);//40);
487 
488  //
489  // Compute inverse solution
490  //
491  PwlRapMusic t_pwlRapMusic(t_clusteredFwd, false, numDipolePairs);
492 
493  if(doMovie)
494  t_pwlRapMusic.setStcAttr(200,0.5);
495 
496 
497 #ifdef BENCHMARK
498  MNESourceEstimate sourceEstimate;
499  QList<qint64> qVecElapsedTime;
500  for(qint32 i = 0; i < 100; ++i)
501  {
502  //Benchmark time
503  QElapsedTimer timer;
504  timer.start();
505  sourceEstimate = t_pwlRapMusic.calculateInverse(pickedEvoked);
506  qVecElapsedTime.append(timer.elapsed());
507  }
508 
509  double meanTime = 0.0;
510  qint32 offset = 19;
511  qint32 c = 0;
512  for(qint32 i = offset; i < qVecElapsedTime.size(); ++i)
513  {
514  meanTime += qVecElapsedTime[i];
515  ++c;
516  }
517 
518  meanTime /= (double)c;
519 
520  double varTime = 0;
521  for(qint32 i = offset; i < qVecElapsedTime.size(); ++i)
522  varTime += pow(qVecElapsedTime[i] - meanTime,2);
523 
524  varTime /= (double)c - 1.0f;
525  varTime = sqrt(varTime);
526 
527  qDebug() << "RAP-MUSIC calculation took" << meanTime << "+-" << varTime << "ms in average";
528 
529 #else
530  MNESourceEstimate sourceEstimate = t_pwlRapMusic.calculateInverse(pickedEvoked);
531 #endif
532 
533  if(sourceEstimate.isEmpty())
534  return 1;
535 
536 // // View activation time-series
537 // std::cout << "\nsourceEstimate:\n" << sourceEstimate.data.block(0,0,10,10) << std::endl;
538 // std::cout << "time\n" << sourceEstimate.times.block(0,0,1,10) << std::endl;
539 // std::cout << "timeMin\n" << sourceEstimate.times[0] << std::endl;
540 // std::cout << "timeMax\n" << sourceEstimate.times[sourceEstimate.times.size()-1] << std::endl;
541 // std::cout << "time step\n" << sourceEstimate.tstep << std::endl;
542 
543  //Source Estimate end
544  //########################################################################################
545 
546 // //only one time point - P100
547 // qint32 sample = 0;
548 // for(qint32 i = 0; i < sourceEstimate.times.size(); ++i)
549 // {
550 // if(sourceEstimate.times(i) >= 0)
551 // {
552 // sample = i;
553 // break;
554 // }
555 // }
556 // sample += (qint32)ceil(0.106/sourceEstimate.tstep); //100ms
557 // sourceEstimate = sourceEstimate.reduce(sample, 1);
558 
559  QList<Label> t_qListLabels;
560  QList<RowVector4i> t_qListRGBAs;
561 
562  //ToDo overload toLabels using instead of t_surfSet rr of MNESourceSpace
563  t_annotationSet.toLabels(t_surfSet, t_qListLabels, t_qListRGBAs);
564 
565  InverseView view(t_pwlRapMusic.getSourceSpace(), t_qListLabels, t_qListRGBAs, 24, true, false, false);//true);
566 
567  if (view.stereoType() != QGLView::RedCyanAnaglyph)
568  view.camera()->setEyeSeparation(0.3f);
569  QStringList args = QCoreApplication::arguments();
570  int w_pos = args.indexOf("-width");
571  int h_pos = args.indexOf("-height");
572  if (w_pos >= 0 && h_pos >= 0)
573  {
574  bool ok = true;
575  int w = args.at(w_pos + 1).toInt(&ok);
576  if (!ok)
577  {
578  qWarning() << "Could not parse width argument:" << args;
579  return 1;
580  }
581  int h = args.at(h_pos + 1).toInt(&ok);
582  if (!ok)
583  {
584  qWarning() << "Could not parse height argument:" << args;
585  return 1;
586  }
587  view.resize(w, h);
588  }
589  else
590  {
591  view.resize(800, 600);
592  }
593  view.show();
594 
595  //Push Estimate
596  view.pushSourceEstimate(sourceEstimate);
597 
598  if(!t_sFileNameStc.isEmpty())
599  {
600  QFile t_fileClusteredStc(t_sFileNameStc);
601  sourceEstimate.write(t_fileClusteredStc);
602  }
603 
604  return a.exec();//1;//a.exec();
605 }
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
MNEForwardSolution cluster_forward_solution(const AnnotationSet &p_AnnotationSet, qint32 p_iClusterSize, MatrixXd &p_D=defaultD, const FiffCov &p_pNoise_cov=defaultCov, const FiffInfo &p_pInfo=defaultInfo, QString p_sMethod="cityblock") const
MNEEpochDataList class declaration.
3D stereoscopic labels
Definition: inverseview.h:118
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
PwlRapMusic algorithm class declaration.
Annotation set.
Definition: annotationset.h:96
bool write(QIODevice &p_IODevice)
Definition: fiff.h:98
FiffEvoked pick_channels(const QStringList &include=defaultQStringList, const QStringList &exclude=defaultQStringList) const
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
AnnotationSet class declaration.
InverseView class declaration.
Label class declaration.
The PwlRapMusic class provides the POWELL RAP MUSIC Algorithm CPU implementation. ToDo: Paper referen...
Definition: pwlrapmusic.h:96
MNEMath class declaration.
FIFF class declaration, which provides static wrapper functions to stay consistent with mne matlab to...
A hemisphere set of surfaces.
Definition: surfaceset.h:83