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