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_fileFwd("./MNE-sample-data/MEG/sample/sample_audvis-meg-eeg-oct-6-fwd.fif");
106 // QFile t_fileCov("./MNE-sample-data/MEG/sample/sample_audvis-cov.fif");
107 // QFile t_fileRaw("./MNE-sample-data/MEG/sample/sample_audvis_raw.fif");
108 // QString t_sEventName = "./MNE-sample-data/MEG/sample/sample_audvis_raw-eve.fif";
109 // AnnotationSet t_annotationSet("./MNE-sample-data/subjects/sample/label/lh.aparc.a2009s.annot", "./MNE-sample-data/subjects/sample/label/rh.aparc.a2009s.annot");
110 // SurfaceSet t_surfSet("./MNE-sample-data/subjects/sample/surf/lh.white", "./MNE-sample-data/subjects/sample/surf/rh.white");
111 
112  QFile t_fileFwd("D:/Data/MEG/mind006/mind006_051209_auditory01_raw-oct-6p-fwd.fif");
113  QFile t_fileCov("D:/Data/MEG/mind006/mind006_051209_auditory01_raw-cov.fif");
114  QFile t_fileRaw("D:/Data/MEG/mind006/mind006_051209_auditory01_raw.fif");
115  QString t_sEventName = "D:/Data/MEG/mind006/mind006_051209_auditory01_raw-eve.fif";
116  AnnotationSet t_annotationSet("mind006", 2, "aparc.a2009s", "D:/Data/subjects");
117  SurfaceSet t_surfSet("mind006", 2, "white", "D:/Data/subjects");
118 
119 
120  qint32 event = 1;
121 
122  float tmin = -0.2f;
123  float tmax = 0.4f;
124 
125  bool keep_comp = false;
126  fiff_int_t dest_comp = 0;
127  bool pick_all = true;
128 
129  qint32 k, p;
130 
131  //
132  // Setup for reading the raw data
133  //
134  FiffRawData raw(t_fileRaw);
135 
136  RowVectorXi picks;
137  if (pick_all)
138  {
139  //
140  // Pick all
141  //
142  picks.resize(raw.info.nchan);
143 
144  for(k = 0; k < raw.info.nchan; ++k)
145  picks(k) = k;
146  //
147  }
148  else
149  {
150  QStringList include;
151  include << "STI 014";
152  bool want_meg = true;
153  bool want_eeg = false;
154  bool want_stim = false;
155 
156 // picks = Fiff::pick_types(raw.info, want_meg, want_eeg, want_stim, include, raw.info.bads);
157  picks = raw.info.pick_types(want_meg, want_eeg, want_stim, include, raw.info.bads);//prefer member function
158  }
159 
160  QStringList ch_names;
161  for(k = 0; k < picks.cols(); ++k)
162  ch_names << raw.info.ch_names[picks(0,k)];
163 
164  //
165  // Set up projection
166  //
167  if (raw.info.projs.size() == 0)
168  printf("No projector specified for these data\n");
169  else
170  {
171  //
172  // Activate the projection items
173  //
174  for (k = 0; k < raw.info.projs.size(); ++k)
175  raw.info.projs[k].active = true;
176 
177  printf("%d projection items activated\n",raw.info.projs.size());
178  //
179  // Create the projector
180  //
181 // fiff_int_t nproj = MNE::make_projector_info(raw.info, raw.proj); Using the member function instead
182  fiff_int_t nproj = raw.info.make_projector(raw.proj);
183 
184  if (nproj == 0)
185  {
186  printf("The projection vectors do not apply to these channels\n");
187  }
188  else
189  {
190  printf("Created an SSP operator (subspace dimension = %d)\n",nproj);
191  }
192  }
193 
194  //
195  // Set up the CTF compensator
196  //
197 // qint32 current_comp = MNE::get_current_comp(raw.info);
198  qint32 current_comp = raw.info.get_current_comp();
199  if (current_comp > 0)
200  printf("Current compensation grade : %d\n",current_comp);
201 
202  if (keep_comp)
203  dest_comp = current_comp;
204 
205  if (current_comp != dest_comp)
206  {
207  qDebug() << "This part needs to be debugged";
208  if(MNE::make_compensator(raw.info, current_comp, dest_comp, raw.comp))
209  {
210 // raw.info.chs = MNE::set_current_comp(raw.info.chs,dest_comp);
211  raw.info.set_current_comp(dest_comp);
212  printf("Appropriate compensator added to change to grade %d.\n",dest_comp);
213  }
214  else
215  {
216  printf("Could not make the compensator\n");
217  return 0;
218  }
219  }
220  //
221  // Read the events
222  //
223  QFile t_EventFile;
224  MatrixXi events;
225  if (t_sEventName.size() == 0)
226  {
227  p = t_fileRaw.fileName().indexOf(".fif");
228  if (p > 0)
229  {
230  t_sEventName = t_fileRaw.fileName().replace(p, 4, "-eve.fif");
231  }
232  else
233  {
234  printf("Raw file name does not end properly\n");
235  return 0;
236  }
237 // events = mne_read_events(t_sEventName);
238 
239  t_EventFile.setFileName(t_sEventName);
240  MNE::read_events(t_EventFile, events);
241  printf("Events read from %s\n",t_sEventName.toUtf8().constData());
242  }
243  else
244  {
245  //
246  // Binary file
247  //
248  p = t_fileRaw.fileName().indexOf(".fif");
249  if (p > 0)
250  {
251  t_EventFile.setFileName(t_sEventName);
252  if(!MNE::read_events(t_EventFile, events))
253  {
254  printf("Error while read events.\n");
255  return 0;
256  }
257  printf("Binary event file %s read\n",t_sEventName.toUtf8().constData());
258  }
259  else
260  {
261  //
262  // Text file
263  //
264  printf("Text file %s is not supported jet.\n",t_sEventName.toUtf8().constData());
265 // try
266 // events = load(eventname);
267 // catch
268 // error(me,mne_omit_first_line(lasterr));
269 // end
270 // if size(events,1) < 1
271 // error(me,'No data in the event file');
272 // end
273 // //
274 // // Convert time to samples if sample number is negative
275 // //
276 // for p = 1:size(events,1)
277 // if events(p,1) < 0
278 // events(p,1) = events(p,2)*raw.info.sfreq;
279 // end
280 // end
281 // //
282 // // Select the columns of interest (convert to integers)
283 // //
284 // events = int32(events(:,[1 3 4]));
285 // //
286 // // New format?
287 // //
288 // if events(1,2) == 0 && events(1,3) == 0
289 // fprintf(1,'The text event file %s is in the new format\n',eventname);
290 // if events(1,1) ~= raw.first_samp
291 // error(me,'This new format event file is not compatible with the raw data');
292 // end
293 // else
294 // fprintf(1,'The text event file %s is in the old format\n',eventname);
295 // //
296 // // Offset with first sample
297 // //
298 // events(:,1) = events(:,1) + raw.first_samp;
299 // end
300  }
301  }
302 
303  //
304  // Select the desired events
305  //
306  qint32 count = 0;
307  MatrixXi selected = MatrixXi::Zero(1, events.rows());
308  for (p = 0; p < events.rows(); ++p)
309  {
310  if (events(p,1) == 0 && events(p,2) == event)
311  {
312  selected(0,count) = p;
313  ++count;
314  }
315  }
316  selected.conservativeResize(1, count);
317  if (count > 0)
318  printf("%d matching events found\n",count);
319  else
320  {
321  printf("No desired events found.\n");
322  return 0;
323  }
324 
325 
326  fiff_int_t event_samp, from, to;
327  MatrixXd timesDummy;
328 
329  MNEEpochDataList data;
330 
331  MNEEpochData* epoch = NULL;
332 
333  MatrixXd times;
334 
335  for (p = 0; p < count; ++p)
336  {
337  //
338  // Read a data segment
339  //
340  event_samp = events(selected(p),0);
341  from = event_samp + tmin*raw.info.sfreq;
342  to = event_samp + floor(tmax*raw.info.sfreq + 0.5);
343 
344  epoch = new MNEEpochData();
345 
346  if(raw.read_raw_segment(epoch->epoch, timesDummy, from, to, picks))
347  {
348  if (p == 0)
349  {
350  times.resize(1, to-from+1);
351  for (qint32 i = 0; i < times.cols(); ++i)
352  times(0, i) = ((float)(from-event_samp+i)) / raw.info.sfreq;
353  }
354 
355  epoch->event = event;
356  epoch->tmin = ((float)(from)-(float)(raw.first_samp))/raw.info.sfreq;
357  epoch->tmax = ((float)(to)-(float)(raw.first_samp))/raw.info.sfreq;
358 
359  data.append(MNEEpochData::SPtr(epoch));//List takes ownwership of the pointer - no delete need
360  }
361  else
362  {
363  printf("Can't read the event data segments");
364  return 0;
365  }
366  }
367 
368  if(data.size() > 0)
369  {
370  printf("Read %d epochs, %d samples each.\n",data.size(),(qint32)data[0]->epoch.cols());
371 
372  //DEBUG
373  std::cout << data[0]->epoch.block(0,0,10,10) << std::endl;
374  qDebug() << data[0]->epoch.rows() << " x " << data[0]->epoch.cols();
375 
376  std::cout << times.block(0,0,1,10) << std::endl;
377  qDebug() << times.rows() << " x " << times.cols();
378  }
379 
380  //
381  // calculate the average
382  //
383 // //Option 1
384 // qint32 numAverages = 99;
385 // VectorXi vecSel(numAverages);
386 // srand (time(NULL)); // initialize random seed
387 
388 // for(qint32 i = 0; i < vecSel.size(); ++i)
389 // {
390 // qint32 val = rand() % data.size();
391 // vecSel(i) = val;
392 // }
393 
394 // //Option 2
395 // VectorXi vecSel(20);
396 
398 
399 // vecSel << 65, 22, 47, 55, 16, 29, 14, 36, 57, 97, 89, 46, 9, 93, 83, 52, 71, 52, 3, 96;
400 
401  //Option 3 Newest
402  VectorXi vecSel(10);
403 
404  vecSel << 0, 96, 80, 55, 66, 25, 26, 2, 55, 58, 6, 88;
405 
406 
407  std::cout << "Select following epochs to average:\n" << vecSel << std::endl;
408 
409  FiffEvoked evoked = data.average(raw.info, tmin*raw.info.sfreq, floor(tmax*raw.info.sfreq + 0.5), vecSel);
410 
411 
412 
413  //########################################################################################
414  // Source Estimate
415 
416  double snr = 1.0f;//0.1f;//1.0f;//3.0f;//0.1f;//3.0f;
417  QString method("dSPM"); //"MNE" | "dSPM" | "sLORETA"
418 
419  QString t_sFileNameClusteredInv("");
420  QString t_sFileNameStc("test_mind006_051209_auditory01.stc");
421 
422  // Parse command line parameters
423  for(qint32 i = 0; i < argc; ++i)
424  {
425  if(strcmp(argv[i], "-snr") == 0 || strcmp(argv[i], "--snr") == 0)
426  {
427  if(i + 1 < argc)
428  snr = atof(argv[i+1]);
429  }
430  else if(strcmp(argv[i], "-method") == 0 || strcmp(argv[i], "--method") == 0)
431  {
432  if(i + 1 < argc)
433  method = QString::fromUtf8(argv[i+1]);
434  }
435  else if(strcmp(argv[i], "-inv") == 0 || strcmp(argv[i], "--inv") == 0)
436  {
437  if(i + 1 < argc)
438  t_sFileNameClusteredInv = QString::fromUtf8(argv[i+1]);
439  }
440  else if(strcmp(argv[i], "-stc") == 0 || strcmp(argv[i], "--stc") == 0)
441  {
442  if(i + 1 < argc)
443  t_sFileNameStc = QString::fromUtf8(argv[i+1]);
444  }
445  }
446 
447  double lambda2 = 1.0 / pow(snr, 2);
448  qDebug() << "Start calculation with: SNR" << snr << "; Lambda" << lambda2 << "; Method" << method << "; stc:" << t_sFileNameStc;
449 
450 // // Load data
451 // fiff_int_t setno = 1;
452 // QPair<QVariant, QVariant> baseline(QVariant(), 0);
453 // FiffEvoked evoked(t_fileEvoked, setno, baseline);
454 // if(evoked.isEmpty())
455 // return 1;
456 
457  MNEForwardSolution t_Fwd(t_fileFwd);
458  if(t_Fwd.isEmpty())
459  return 1;
460 
461 
462  FiffCov noise_cov(t_fileCov);
463 
464  // regularize noise covariance
465  noise_cov = noise_cov.regularize(evoked.info, 0.05, 0.05, 0.1, true);
466 
467  //
468  // Cluster forward solution;
469  //
470  MatrixXd D;
471  MNEForwardSolution t_clusteredFwd = t_Fwd.cluster_forward_solution(t_annotationSet, 20, D, noise_cov, evoked.info);
472 
473  //
474  // make an inverse operators
475  //
476  FiffInfo info = evoked.info;
477 
478  MNEInverseOperator inverse_operator(info, t_clusteredFwd, noise_cov, 0.2f, 0.8f);
479 
480  //
481  // save clustered inverse
482  //
483  if(!t_sFileNameClusteredInv.isEmpty())
484  {
485  QFile t_fileClusteredInverse(t_sFileNameClusteredInv);
486  inverse_operator.write(t_fileClusteredInverse);
487  }
488 
489  //
490  // Compute inverse solution
491  //
492  MinimumNorm minimumNorm(inverse_operator, lambda2, method);
493 
494 #ifdef BENCHMARK
495  //
496  // Set up the inverse according to the parameters
497  //
498  minimumNorm.doInverseSetup(vecSel.size(),false);
499 
500  MNESourceEstimate sourceEstimate;
501  QList<qint64> qVecElapsedTime;
502  for(qint32 i = 0; i < 100; ++i)
503  {
504  //Benchmark time
505  QElapsedTimer timer;
506  timer.start();
507  sourceEstimate = minimumNorm.calculateInverse(evoked.data, evoked.times(0), evoked.times(1)-evoked.times(0));
508  qVecElapsedTime.append(timer.elapsed());
509  }
510 
511  double meanTime = 0.0;
512  qint32 offset = 19;
513  qint32 c = 0;
514  for(qint32 i = offset; i < qVecElapsedTime.size(); ++i)
515  {
516  meanTime += qVecElapsedTime[i];
517  ++c;
518  }
519 
520  meanTime /= (double)c;
521 
522  double varTime = 0;
523  for(qint32 i = offset; i < qVecElapsedTime.size(); ++i)
524  varTime += pow(qVecElapsedTime[i] - meanTime,2);
525 
526  varTime /= (double)c - 1.0f;
527  varTime = sqrt(varTime);
528 
529  qDebug() << "MNE calculation took" << meanTime << "+-" << varTime << "ms in average";
530 
531 #else
532  MNESourceEstimate sourceEstimate = minimumNorm.calculateInverse(evoked);
533 #endif
534 
535  if(sourceEstimate.isEmpty())
536  return 1;
537 
538  // View activation time-series
539  std::cout << "\nsourceEstimate:\n" << sourceEstimate.data.block(0,0,10,10) << std::endl;
540  std::cout << "time\n" << sourceEstimate.times.block(0,0,1,10) << std::endl;
541  std::cout << "timeMin\n" << sourceEstimate.times[0] << std::endl;
542  std::cout << "timeMax\n" << sourceEstimate.times[sourceEstimate.times.size()-1] << std::endl;
543  std::cout << "time step\n" << sourceEstimate.tstep << std::endl;
544 
545  //Condition Numbers
546 // MatrixXd mags(102, t_Fwd.sol->data.cols());
547 // qint32 count = 0;
548 // for(qint32 i = 2; i < 306; i += 3)
549 // {
550 // mags.row(count) = t_Fwd.sol->data.row(i);
551 // ++count;
552 // }
553 // MatrixXd magsClustered(102, t_clusteredFwd.sol->data.cols());
554 // count = 0;
555 // for(qint32 i = 2; i < 306; i += 3)
556 // {
557 // magsClustered.row(count) = t_clusteredFwd.sol->data.row(i);
558 // ++count;
559 // }
560 
561 // MatrixXd grads(204, t_Fwd.sol->data.cols());
562 // count = 0;
563 // for(qint32 i = 0; i < 306; i += 3)
564 // {
565 // grads.row(count) = t_Fwd.sol->data.row(i);
566 // ++count;
567 // grads.row(count) = t_Fwd.sol->data.row(i+1);
568 // ++count;
569 // }
570 // MatrixXd gradsClustered(204, t_clusteredFwd.sol->data.cols());
571 // count = 0;
572 // for(qint32 i = 0; i < 306; i += 3)
573 // {
574 // gradsClustered.row(count) = t_clusteredFwd.sol->data.row(i);
575 // ++count;
576 // gradsClustered.row(count) = t_clusteredFwd.sol->data.row(i+1);
577 // ++count;
578 // }
579 
580  VectorXd s;
581 
582  double t_dConditionNumber = MNEMath::getConditionNumber(t_Fwd.sol->data, s);
583  double t_dConditionNumberClustered = MNEMath::getConditionNumber(t_clusteredFwd.sol->data, s);
584 
585 
586  std::cout << "Condition Number:\n" << t_dConditionNumber << std::endl;
587  std::cout << "Clustered Condition Number:\n" << t_dConditionNumberClustered << std::endl;
588 
589  std::cout << "ForwardSolution" << t_Fwd.sol->data.block(0,0,10,10) << std::endl;
590 
591  std::cout << "Clustered ForwardSolution" << t_clusteredFwd.sol->data.block(0,0,10,10) << std::endl;
592 
593 
594 // double t_dConditionNumberMags = MNEMath::getConditionNumber(mags, s);
595 // double t_dConditionNumberMagsClustered = MNEMath::getConditionNumber(magsClustered, s);
596 
597 // std::cout << "Condition Number Magnetometers:\n" << t_dConditionNumberMags << std::endl;
598 // std::cout << "Clustered Condition Number Magnetometers:\n" << t_dConditionNumberMagsClustered << std::endl;
599 
600 // double t_dConditionNumberGrads = MNEMath::getConditionNumber(grads, s);
601 // double t_dConditionNumberGradsClustered = MNEMath::getConditionNumber(gradsClustered, s);
602 
603 // std::cout << "Condition Number Gradiometers:\n" << t_dConditionNumberGrads << std::endl;
604 // std::cout << "Clustered Condition Number Gradiometers:\n" << t_dConditionNumberGradsClustered << std::endl;
605 
606 
607  //Source Estimate end
608  //########################################################################################
609 
610 // //only one time point - P100
611 // qint32 sample = 0;
612 // for(qint32 i = 0; i < sourceEstimate.times.size(); ++i)
613 // {
614 // if(sourceEstimate.times(i) >= 0)
615 // {
616 // sample = i;
617 // break;
618 // }
619 // }
620 // sample += (qint32)ceil(0.106/sourceEstimate.tstep); //100ms
621 // sourceEstimate = sourceEstimate.reduce(sample, 1);
622 
623  QList<Label> t_qListLabels;
624  QList<RowVector4i> t_qListRGBAs;
625 
626  //ToDo overload toLabels using instead of t_surfSet rr of MNESourceSpace
627  t_annotationSet.toLabels(t_surfSet, t_qListLabels, t_qListRGBAs);
628 
629  InverseView view(minimumNorm.getSourceSpace(), t_qListLabels, t_qListRGBAs, 24, true, false, true);
630 
631  if (view.stereoType() != QGLView::RedCyanAnaglyph)
632  view.camera()->setEyeSeparation(0.3f);
633  QStringList args = QCoreApplication::arguments();
634  int w_pos = args.indexOf("-width");
635  int h_pos = args.indexOf("-height");
636  if (w_pos >= 0 && h_pos >= 0)
637  {
638  bool ok = true;
639  int w = args.at(w_pos + 1).toInt(&ok);
640  if (!ok)
641  {
642  qWarning() << "Could not parse width argument:" << args;
643  return 1;
644  }
645  int h = args.at(h_pos + 1).toInt(&ok);
646  if (!ok)
647  {
648  qWarning() << "Could not parse height argument:" << args;
649  return 1;
650  }
651  view.resize(w, h);
652  }
653  else
654  {
655  view.resize(800, 600);
656  }
657  view.show();
658 
659  //Push Estimate
660  view.pushSourceEstimate(sourceEstimate);
661 
662  if(!t_sFileNameStc.isEmpty())
663  {
664  QFile t_fileClusteredStc(t_sFileNameStc);
665  sourceEstimate.write(t_fileClusteredStc);
666  }
667 
668 //*/
669  return a.exec();//1;//a.exec();
670 }
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
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
RowVectorXf times
Definition: fiff_evoked.h:216
static bool read_events(QIODevice &p_IODevice, MatrixXi &eventlist)
Definition: mne.cpp:68
FIFF raw measurement data.
Definition: fiff_raw_data.h:94
static double getConditionNumber(const MatrixXd &A, VectorXd &s)
Definition: mnemath.cpp:103
Minimum norm class declaration.
Annotation set.
Definition: annotationset.h:96
bool write(QIODevice &p_IODevice)
covariance data
Definition: fiff_cov.h:94
FiffNamedMatrix::SDPtr sol
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.
InverseView class declaration.
Label class declaration.
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