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