MNE-CPP  beta 1.0
main.cpp
1 //=============================================================================================================
36 //*************************************************************************************************************
37 //=============================================================================================================
38 // INCLUDES
39 //=============================================================================================================
40 
41 #include <fs/label.h>
42 #include <fs/surfaceset.h>
43 #include <fs/annotationset.h>
44 
45 #include <fiff/fiff_evoked.h>
46 #include <fiff/fiff.h>
47 #include <mne/mne.h>
48 
49 #include <mne/mne_sourceestimate.h>
51 
52 #include <disp3D/inverseview.h>
53 
54 #include <utils/mnemath.h>
55 
56 #include <iostream>
57 
58 
59 //*************************************************************************************************************
60 //=============================================================================================================
61 // QT INCLUDES
62 //=============================================================================================================
63 
64 #include <QGuiApplication>
65 #include <QSet>
66 
67 
68 //*************************************************************************************************************
69 //=============================================================================================================
70 // USED NAMESPACES
71 //=============================================================================================================
72 
73 using namespace MNELIB;
74 using namespace FSLIB;
75 using namespace FIFFLIB;
76 using namespace INVERSELIB;
77 using namespace DISP3DLIB;
78 using namespace UTILSLIB;
79 
80 
81 //*************************************************************************************************************
82 //=============================================================================================================
83 // MAIN
84 //=============================================================================================================
85 
86 //=============================================================================================================
95 int main(int argc, char *argv[])
96 {
97  QGuiApplication a(argc, argv);
98 
99 // QFile t_fileRaw("./MNE-sample-data/MEG/sample/sample_audvis_raw.fif");
100 // QFile t_fileFwd("./MNE-sample-data/MEG/sample/sample_audvis-eeg-oct-6-fwd.fif");
101 // AnnotationSet t_annotationSet("./MNE-sample-data/subjects/sample/label/lh.aparc.a2009s.annot", "./MNE-sample-data/subjects/sample/label/rh.aparc.a2009s.annot");
102 // SurfaceSet t_surfSet("./MNE-sample-data/subjects/sample/surf/lh.white", "./MNE-sample-data/subjects/sample/surf/rh.white");
103 
104 // QFile t_fileRaw("E:/Data/sl_data/MEG/mind006/mind006_051209_auditory01_raw.fif");
105 // QFile t_fileFwd("E:/Data/sl_data/MEG/mind006/mind006_051209_auditory01_raw-oct-6p-fwd.fif");
106 // 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");
107 // SurfaceSet t_surfSet("E:/Data/sl_data/subjects/mind006/surf/lh.white", "E:/Data/sl_data/subjects/mind006/surf/rh.white");
108 
109 // QFile t_fileRaw("E:/Data/sl_data/MEG/mind006/mind006_051209_median01_raw.fif");
110 // QFile t_fileFwd("E:/Data/sl_data/MEG/mind006/mind006_051209_median01_raw-oct-6-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("D:/Dropbox/Masterarbeit DB/Messdaten/mind006/mind006_051209_index01_raw.fif");
115 // QFile t_fileFwd("D:/Dropbox/Masterarbeit DB/Messdaten/mind006/mind006_051209_index01_raw-oct-6-fwd.fif");
116 // AnnotationSet t_annotationSet("D:/Dropbox/Masterarbeit DB/Messdaten/mind006/label/lh.aparc.a2009s.annot", "D:/Dropbox/Masterarbeit DB/Messdaten/mind006/label/rh.aparc.a2009s.annot");
117 // SurfaceSet t_surfSet("D:/Dropbox/Masterarbeit DB/Messdaten/mind006/surf/lh.white", "D:/Dropbox/Masterarbeit DB/Messdaten/mind006/surf/rh.white");
118 // QString t_sFileNameStc("D:/Dropbox/Masterarbeit DB/Messdaten/mind006/mind006_051209_index01_raw_RAP.stc");
119 
120 // QFile t_fileRaw("D:/Dropbox/Masterarbeit DB/Messdaten/EEG/2014_03_17_Lorenz_Esch_010/Processed/averaged/EEG_data_001_voluntary_left_tapping_raw_Averaged_-1_1_108.fif");
121 // QFile t_fileFwd("D:/Dropbox/Masterarbeit DB/Messdaten/Forward solutions/Lorenz-140128-Duke128-fwd.fif");
122 // AnnotationSet t_annotationSet("D:/Dropbox/Masterarbeit DB/Messdaten/Forward solutions/atlas/lh.aparc.a2009s.annot", "D:/Dropbox/Masterarbeit DB/Messdaten/Forward solutions/atlas/rh.aparc.a2009s.annot");
123 // SurfaceSet t_surfSet("D:/Dropbox/Masterarbeit DB/Messdaten/Forward solutions/surface/lh.white", "D:/Dropbox/Masterarbeit DB/Messdaten/Forward solutions/surface/rh.white");
124 // QString t_sFileNameStc("D:/Dropbox/Masterarbeit DB/Messdaten/EEG/2014_03_17_Lorenz_Esch_010/Processed/stc/EEG_data_001_voluntary_left_tapping_raw_Averaged_-1_1_108_RAP.stc");
125 
126  QFile t_fileRaw("D:/Dropbox/Masterarbeit DB/Messdaten/EEG/2014_02_04_Lorenz_Esch_005/Processed/averaged/EEG_data_002_voluntary_right_tapping_filtered_07_40_Averaged_-1_1_184.fif");
127  QFile t_fileFwd("D:/Dropbox/Masterarbeit DB/Messdaten/Forward solutions/Lorenz-140128-Duke128-fwd.fif");
128  AnnotationSet t_annotationSet("D:/Dropbox/Masterarbeit DB/Messdaten/Forward solutions/atlas/lh.aparc.a2009s.annot", "D:/Dropbox/Masterarbeit DB/Messdaten/Forward solutions/atlas/rh.aparc.a2009s.annot");
129  SurfaceSet t_surfSet("D:/Dropbox/Masterarbeit DB/Messdaten/Forward solutions/surface/lh.white", "D:/Dropbox/Masterarbeit DB/Messdaten/Forward solutions/surface/rh.white");
130  QString t_sFileNameStc("D:/Dropbox/Masterarbeit DB/Messdaten/EEG/2014_02_04_Lorenz_Esch_005/Processed/stc/EEG_data_002_voluntary_right_tapping_filtered_07_40_Averaged_-1_1_184_RAP.stc");
131 
132 // QFile t_fileRaw("D:/Dropbox/Masterarbeit DB/Messdaten/EEG/2014_03_17_Lorenz_Esch_010/Processed/filtered/EEG_data_001_voluntary_left_right_tapping_filtered_07_40_raw.fif");
133 // QFile t_fileFwd("D:/Dropbox/Masterarbeit DB/Messdaten/Forward solutions/Lorenz-140128-Duke128-fwd.fif");
134 // AnnotationSet t_annotationSet("D:/Dropbox/Masterarbeit DB/Messdaten/Forward solutions/atlas/lh.aparc.a2009s.annot", "D:/Dropbox/Masterarbeit DB/Messdaten/Forward solutions/atlas/rh.aparc.a2009s.annot");
135 // SurfaceSet t_surfSet("D:/Dropbox/Masterarbeit DB/Messdaten/Forward solutions/surface/lh.white", "D:/Dropbox/Masterarbeit DB/Messdaten/Forward solutions/surface/rh.white");
136 // QString t_sFileNameStc("D:/Dropbox/Masterarbeit DB/Messdaten/EEG/2014_03_17_Lorenz_Esch_010/Processed/stc/EEG_data_001_voluntary_left_right_tapping_filtered_07_40_raw_RAP.stc");
137 
138  qint32 numDipolePairs = 7;
139 
140  qint32 samplesStcWindow = 300;
141  float stcOverlap = 0.0f;
142 
143  qint32 startSample = 0;
144  qint32 numSample = 40000;
145 
146  bool in_samples = true;
147  bool keep_comp = true;
148 
149  // Parse command line parameters
150  for(qint32 i = 0; i < argc; ++i)
151  {
152  if(strcmp(argv[i], "-stc") == 0 || strcmp(argv[i], "--stc") == 0)
153  {
154  if(i + 1 < argc)
155  t_sFileNameStc = QString::fromUtf8(argv[i+1]);
156  }
157  }
158 
159  //
160  // Read raw data
161  //
162  MNEForwardSolution t_Fwd(t_fileFwd);
163  if(t_Fwd.isEmpty())
164  return 1;
165 
166  QList<Label> t_qListLabels;
167  QList<RowVector4i> t_qListRGBAs;
168 
169  //ToDo overload toLabels using instead of t_surfSet rr of MNESourceSpace
170  t_annotationSet.toLabels(t_surfSet, t_qListLabels, t_qListRGBAs);
171 
172  QList<Label> t_qListLabelSelection;
173 
174  //LH
175  t_qListLabelSelection << t_qListLabels[28] << t_qListLabels[29] << t_qListLabels[45] << t_qListLabels[66];
176  //RH
177  t_qListLabelSelection << t_qListLabels[103] << t_qListLabels[104] << t_qListLabels[120]<< t_qListLabels[141];
178 
179 // std::cout << "t_qListLabels.size " << t_qListLabels.size() << std::endl;
180 
181 // for(qint32 i = 0; i < t_qListLabels.size(); ++i)
182 // {
183 // qDebug() << "Num" << i << t_qListLabels[i].name;
184 // std::cout << t_qListLabels[i].hemi << std::endl;
185 // }
186 
187  MNEForwardSolution t_SelectFwd = t_Fwd.pick_regions(t_qListLabelSelection);
188 
189  //
190  // Setup for reading the raw data
191  //
192  FiffRawData raw(t_fileRaw);
193 
194  float from = raw.first_samp + startSample;
195  float to = from + numSample - 1;
196 
197  //
198  // Set up pick list: MEG + STI 014 - bad channels
199  //
200  QStringList include;
201 // include << "STI 014";
202  bool want_meg = false;
203  bool want_eeg = true;
204  bool want_stim = false;
205 
206  RowVectorXi picks = raw.info.pick_types(want_meg, want_eeg, want_stim, include);//, raw.info.bads);
207 
208  //
209  // Set up projection
210  //
211  qint32 k = 0;
212  if (raw.info.projs.size() == 0)
213  printf("No projector specified for these data\n");
214  else
215  {
216  //
217  // Activate the projection items
218  //
219  for (k = 0; k < raw.info.projs.size(); ++k)
220  raw.info.projs[k].active = true;
221 
222  printf("%d projection items activated\n",raw.info.projs.size());
223  //
224  // Create the projector
225  //
226  fiff_int_t nproj = raw.info.make_projector(raw.proj);
227 
228  if (nproj == 0)
229  printf("The projection vectors do not apply to these channels\n");
230  else
231  printf("Created an SSP operator (subspace dimension = %d)\n",nproj);
232  }
233 
234  //
235  // Set up the CTF compensator
236  //
237  qint32 current_comp = raw.info.get_current_comp();
238  qint32 dest_comp = -1;
239 
240  if (current_comp > 0)
241  printf("Current compensation grade : %d\n",current_comp);
242 
243  if (keep_comp)
244  dest_comp = current_comp;
245 
246  if (current_comp != dest_comp)
247  {
248  qDebug() << "This part needs to be debugged";
249  if(MNE::make_compensator(raw.info, current_comp, dest_comp, raw.comp))
250  {
251  raw.info.set_current_comp(dest_comp);
252  printf("Appropriate compensator added to change to grade %d.\n",dest_comp);
253  }
254  else
255  {
256  printf("Could not make the compensator\n");
257  return -1;
258  }
259  }
260  //
261  // Read a data segment
262  // times output argument is optional
263  //
264  bool readSuccessful = false;
265  MatrixXd data;
266  MatrixXd times;
267  if (in_samples)
268  readSuccessful = raw.read_raw_segment(data, times, (qint32)from, (qint32)to, picks);
269  else
270  readSuccessful = raw.read_raw_segment_times(data, times, from, to, picks);
271 
272  if (!readSuccessful)
273  {
274  printf("Could not read raw segment.\n");
275  return -1;
276  }
277 
278  printf("Read %d samples.\n",(qint32)data.cols());
279 
280  //########################################################################################
281  // RAP MUSIC Source Estimate
282 
283  //
284  // Cluster forward solution;
285  //
286  MNEForwardSolution t_clusteredFwd = t_SelectFwd.cluster_forward_solution(t_annotationSet, 20);//t_Fwd.cluster_forward_solution_ccr(t_annotationSet, 20);//40);
287 
288 
289  //
290  // Compute inverse solution
291  //
292  PwlRapMusic t_pwlRapMusic(t_clusteredFwd, false, numDipolePairs);
293 
294  MNESourceEstimate sourceEstimate;
295 
296  float tstep = 1.0f/raw.info.sfreq;
297 
298  float tmin;
299  if(in_samples)
300  tmin = from * tstep;
301  else
302  tmin = from;
303 
304  //
305  // Rap MUSIC Source estimate
306  //
307  sourceEstimate.data = MatrixXd::Zero(t_clusteredFwd.nsource, data.cols());
308 
309  //Results
310  sourceEstimate.vertices = VectorXi(t_clusteredFwd.src[0].vertno.size() + t_clusteredFwd.src[1].vertno.size());
311  sourceEstimate.vertices << t_clusteredFwd.src[0].vertno, t_clusteredFwd.src[1].vertno;
312 
313  sourceEstimate.times = RowVectorXf::Zero(data.cols());
314  sourceEstimate.times[0] = tmin;
315  for(qint32 i = 1; i < sourceEstimate.times.size(); ++i)
316  sourceEstimate.times[i] = sourceEstimate.times[i-1] + tstep;
317  sourceEstimate.tmin = tmin;
318  sourceEstimate.tstep = tstep;
319 
320 
321  bool first = true;
322  bool last = false;
323 
324  qint32 t_iNumSensors = data.rows();
325  qint32 t_iNumSteps = data.cols();
326 
327  qint32 t_iSamplesOverlap = (qint32)floor(((float)samplesStcWindow)*stcOverlap);
328  qint32 t_iSamplesDiscard = t_iSamplesOverlap/2;
329 
330  MatrixXd measData = MatrixXd::Zero(t_iNumSensors, samplesStcWindow);
331 
332  qint32 curSample = 0;
333  qint32 curResultSample = 0;
334  qint32 stcWindowSize = samplesStcWindow - 2*t_iSamplesDiscard;
335 
336  while(!last)
337  {
338  //Data
339  if(curSample + samplesStcWindow >= t_iNumSteps) //last
340  {
341  last = true;
342  measData = data.block(0, data.cols()-samplesStcWindow, t_iNumSensors, samplesStcWindow);
343  }
344  else
345  measData = data.block(0, curSample, t_iNumSensors, samplesStcWindow);
346 
347 
348  curSample += (samplesStcWindow - t_iSamplesOverlap);
349  if(first)
350  curSample -= t_iSamplesDiscard; //shift on start t_iSamplesDiscard backwards
351 
352  //Calculate
353  MNESourceEstimate stcData = t_pwlRapMusic.calculateInverse(measData, 0.0f, tstep);
354 
355  //Assign Result
356  if(last)
357  stcWindowSize = stcData.data.cols() - curResultSample;
358 
359  sourceEstimate.data.block(0,curResultSample,sourceEstimate.data.rows(),stcWindowSize) =
360  stcData.data.block(0,0,stcData.data.rows(),stcWindowSize);
361 
362  curResultSample += stcWindowSize;
363 
364  if(first)
365  first = false;
366  }
367 
368  if(sourceEstimate.isEmpty())
369  return 1;
370 
371  // Write stc to file
372  if(!t_sFileNameStc.isEmpty())
373  {
374  QFile t_fileClusteredStc(t_sFileNameStc);
375  sourceEstimate.write(t_fileClusteredStc);
376  }
377 
378  //Source Estimate end
379  //########################################################################################
380 
381  InverseView view(t_pwlRapMusic.getSourceSpace(), t_qListLabels, t_qListRGBAs, 24, true, false, false);//true);
382 
383  if (view.stereoType() != QGLView::RedCyanAnaglyph)
384  view.camera()->setEyeSeparation(0.3f);
385  QStringList args = QCoreApplication::arguments();
386  int w_pos = args.indexOf("-width");
387  int h_pos = args.indexOf("-height");
388  if (w_pos >= 0 && h_pos >= 0)
389  {
390  bool ok = true;
391  int w = args.at(w_pos + 1).toInt(&ok);
392  if (!ok)
393  {
394  qWarning() << "Could not parse width argument:" << args;
395  return 1;
396  }
397  int h = args.at(h_pos + 1).toInt(&ok);
398  if (!ok)
399  {
400  qWarning() << "Could not parse height argument:" << args;
401  return 1;
402  }
403  view.resize(w, h);
404  }
405  else
406  {
407  view.resize(800, 600);
408  }
409  view.show();
410 
411  //Push Estimate
412  view.pushSourceEstimate(sourceEstimate);
413 
414  if(!t_sFileNameStc.isEmpty())
415  {
416  QFile t_fileClusteredStc(t_sFileNameStc);
417  sourceEstimate.write(t_fileClusteredStc);
418  }
419 
420  return a.exec();
421 }
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
SurfaceSet class declaration.
3D stereoscopic labels
Definition: inverseview.h:118
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)
MNEForwardSolution pick_regions(const QList< Label > &p_qListLabels) const
Definition: fiff.h:98
MNESourceEstimate class declaration.
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