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 <mne/mne_sourceestimate.h>
49 
50 #include <disp3D/inverseview.h>
51 
52 #include <utils/mnemath.h>
53 
54 #include <iostream>
55 
56 
57 //*************************************************************************************************************
58 //=============================================================================================================
59 // QT INCLUDES
60 //=============================================================================================================
61 
62 #include <QGuiApplication>
63 #include <QSet>
64 
65 
66 //*************************************************************************************************************
67 //=============================================================================================================
68 // USED NAMESPACES
69 //=============================================================================================================
70 
71 using namespace MNELIB;
72 using namespace FSLIB;
73 using namespace FIFFLIB;
74 using namespace INVERSELIB;
75 using namespace DISP3DLIB;
76 using namespace UTILSLIB;
77 
78 
79 //*************************************************************************************************************
80 //=============================================================================================================
81 // MAIN
82 //=============================================================================================================
83 
84 //=============================================================================================================
93 int main(int argc, char *argv[])
94 {
95  QGuiApplication a(argc, argv);
96 
97  //########################################################################################
98  // Source Estimate
99 
100 // QFile t_fileFwd("./MNE-sample-data/MEG/sample/sample_audvis-meg-eeg-oct-6-fwd.fif");
101 // QFile t_fileCov("./MNE-sample-data/MEG/sample/sample_audvis-cov.fif");
102 // QFile t_fileEvoked("./MNE-sample-data/MEG/sample/sample_audvis-ave.fif");
103 
104 // QFile t_fileFwd("/home/chdinh/sl_data/MEG/mind006/mind006_051209_auditory01_raw-oct-6p-fwd.fif");
105 // QFile t_fileCov("/home/chdinh/sl_data/MEG/mind006/mind006_051209_auditory01_raw-cov.fif");
106 // QFile t_fileEvoked("/home/chdinh/sl_data/MEG/mind006/mind006_051209_auditory01_raw-ave.fif");
107 
108  QFile t_fileFwd("E:/Data/sl_data/MEG/mind006/mind006_051209_auditory01_raw-oct-6p-fwd.fif");
109  QFile t_fileCov("E:/Data/sl_data/MEG/mind006/mind006_051209_auditory01_raw-cov.fif");
110  QFile t_fileEvoked("E:/Data/sl_data/MEG/mind006/mind006_051209_auditory01_raw-ave.fif");
111 
112 
113  double snr = 1.0f;//3.0f;//0.1f;//3.0f;
114  QString method("dSPM"); //"MNE" | "dSPM" | "sLORETA"
115 
116  QString t_sFileNameClusteredInv("");
117  QString t_sFileNameStc("");
118 
119  // Parse command line parameters
120  for(qint32 i = 0; i < argc; ++i)
121  {
122  if(strcmp(argv[i], "-snr") == 0 || strcmp(argv[i], "--snr") == 0)
123  {
124  if(i + 1 < argc)
125  snr = atof(argv[i+1]);
126  }
127  else if(strcmp(argv[i], "-method") == 0 || strcmp(argv[i], "--method") == 0)
128  {
129  if(i + 1 < argc)
130  method = QString::fromUtf8(argv[i+1]);
131  }
132  else if(strcmp(argv[i], "-inv") == 0 || strcmp(argv[i], "--inv") == 0)
133  {
134  if(i + 1 < argc)
135  t_sFileNameClusteredInv = QString::fromUtf8(argv[i+1]);
136  }
137  else if(strcmp(argv[i], "-stc") == 0 || strcmp(argv[i], "--stc") == 0)
138  {
139  if(i + 1 < argc)
140  t_sFileNameStc = QString::fromUtf8(argv[i+1]);
141  }
142  }
143 
144  double lambda2 = 1.0 / pow(snr, 2);
145  qDebug() << "Start calculation with: SNR" << snr << "; Lambda" << lambda2 << "; Method" << method << "; stc:" << t_sFileNameStc;
146 
147  // Load data
148  fiff_int_t setno = 0;
149  QPair<QVariant, QVariant> baseline(QVariant(), 0);
150  FiffEvoked evoked(t_fileEvoked, setno, baseline);
151  if(evoked.isEmpty())
152  return 1;
153 
154  std::cout << "evoked first " << evoked.first << "; last " << evoked.last << std::endl;
155 
156  MNEForwardSolution t_Fwd(t_fileFwd);
157  if(t_Fwd.isEmpty())
158  return 1;
159 
160 // AnnotationSet t_annotationSet("./MNE-sample-data/subjects/sample/label/lh.aparc.a2009s.annot", "./MNE-sample-data/subjects/sample/label/rh.aparc.a2009s.annot");
161 // AnnotationSet t_annotationSet("/home/chdinh/sl_data/subjects/mind006/label/lh.aparc.a2009s.annot", "/home/chdinh/sl_data/subjects/mind006/label/rh.aparc.a2009s.annot");
162  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");
163 
164 
165  FiffCov noise_cov(t_fileCov);
166 
167  // regularize noise covariance
168  noise_cov = noise_cov.regularize(evoked.info, 0.05, 0.05, 0.1, true);
169 
170  //
171  // Cluster forward solution;
172  //
173  MNEForwardSolution t_clusteredFwd = t_Fwd.cluster_forward_solution(t_annotationSet, 20);//40);
174 
175 // std::cout << "Size " << t_clusteredFwd.sol->data.rows() << " x " << t_clusteredFwd.sol->data.cols() << std::endl;
176 // std::cout << "Clustered Fwd:\n" << t_clusteredFwd.sol->data.row(0) << std::endl;
177 
178  //
179  // make an inverse operators
180  //
181  FiffInfo info = evoked.info;
182 
183  MNEInverseOperator inverse_operator(info, t_clusteredFwd, noise_cov, 0.2f, 0.8f);
184 
185  //
186  // save clustered inverse
187  //
188  if(!t_sFileNameClusteredInv.isEmpty())
189  {
190  QFile t_fileClusteredInverse(t_sFileNameClusteredInv);
191  inverse_operator.write(t_fileClusteredInverse);
192  }
193 
194  //
195  // Compute inverse solution
196  //
197  MinimumNorm minimumNorm(inverse_operator, lambda2, method);
198  MNESourceEstimate sourceEstimate = minimumNorm.calculateInverse(evoked);
199 
200  if(sourceEstimate.isEmpty())
201  return 1;
202 
203  // View activation time-series
204  std::cout << "\nsourceEstimate:\n" << sourceEstimate.data.block(0,0,10,10) << std::endl;
205  std::cout << "time\n" << sourceEstimate.times.block(0,0,1,10) << std::endl;
206  std::cout << "timeMin\n" << sourceEstimate.times[0] << std::endl;
207  std::cout << "timeMax\n" << sourceEstimate.times[sourceEstimate.times.size()-1] << std::endl;
208  std::cout << "time step\n" << sourceEstimate.tstep << std::endl;
209 
210  //Condition Numbers
211 // MatrixXd mags(102, t_Fwd.sol->data.cols());
212 // qint32 count = 0;
213 // for(qint32 i = 2; i < 306; i += 3)
214 // {
215 // mags.row(count) = t_Fwd.sol->data.row(i);
216 // ++count;
217 // }
218 // MatrixXd magsClustered(102, t_clusteredFwd.sol->data.cols());
219 // count = 0;
220 // for(qint32 i = 2; i < 306; i += 3)
221 // {
222 // magsClustered.row(count) = t_clusteredFwd.sol->data.row(i);
223 // ++count;
224 // }
225 
226 // MatrixXd grads(204, t_Fwd.sol->data.cols());
227 // count = 0;
228 // for(qint32 i = 0; i < 306; i += 3)
229 // {
230 // grads.row(count) = t_Fwd.sol->data.row(i);
231 // ++count;
232 // grads.row(count) = t_Fwd.sol->data.row(i+1);
233 // ++count;
234 // }
235 // MatrixXd gradsClustered(204, t_clusteredFwd.sol->data.cols());
236 // count = 0;
237 // for(qint32 i = 0; i < 306; i += 3)
238 // {
239 // gradsClustered.row(count) = t_clusteredFwd.sol->data.row(i);
240 // ++count;
241 // gradsClustered.row(count) = t_clusteredFwd.sol->data.row(i+1);
242 // ++count;
243 // }
244 
245  VectorXd s;
246 
247  double t_dConditionNumber = MNEMath::getConditionNumber(t_Fwd.sol->data, s);
248  double t_dConditionNumberClustered = MNEMath::getConditionNumber(t_clusteredFwd.sol->data, s);
249 
250 
251  std::cout << "Condition Number:\n" << t_dConditionNumber << std::endl;
252  std::cout << "Clustered Condition Number:\n" << t_dConditionNumberClustered << std::endl;
253 
254  std::cout << "ForwardSolution" << t_Fwd.sol->data.block(0,0,10,10) << std::endl;
255 
256  std::cout << "Clustered ForwardSolution" << t_clusteredFwd.sol->data.block(0,0,10,10) << std::endl;
257 
258 
259 // double t_dConditionNumberMags = MNEMath::getConditionNumber(mags, s);
260 // double t_dConditionNumberMagsClustered = MNEMath::getConditionNumber(magsClustered, s);
261 
262 // std::cout << "Condition Number Magnetometers:\n" << t_dConditionNumberMags << std::endl;
263 // std::cout << "Clustered Condition Number Magnetometers:\n" << t_dConditionNumberMagsClustered << std::endl;
264 
265 // double t_dConditionNumberGrads = MNEMath::getConditionNumber(grads, s);
266 // double t_dConditionNumberGradsClustered = MNEMath::getConditionNumber(gradsClustered, s);
267 
268 // std::cout << "Condition Number Gradiometers:\n" << t_dConditionNumberGrads << std::endl;
269 // std::cout << "Clustered Condition Number Gradiometers:\n" << t_dConditionNumberGradsClustered << std::endl;
270 
271 
272  //Source Estimate end
273  //########################################################################################
274 
275 // SurfaceSet t_surfSet("./MNE-sample-data/subjects/sample/surf/lh.white", "./MNE-sample-data/subjects/sample/surf/rh.white");
276 // SurfaceSet t_surfSet("/home/chdinh/sl_data/subjects/mind006/surf/lh.white", "/home/chdinh/sl_data/subjects/mind006/surf/rh.white");
277  SurfaceSet t_surfSet("E:/Data/sl_data/subjects/mind006/surf/lh.white", "E:/Data/sl_data/subjects/mind006/surf/rh.white");
278 
279 // //only one time point - P100
280 // qint32 sample = 0;
281 // for(qint32 i = 0; i < sourceEstimate.times.size(); ++i)
282 // {
283 // if(sourceEstimate.times(i) >= 0)
284 // {
285 // sample = i;
286 // break;
287 // }
288 // }
289 // sample += (qint32)ceil(0.106/sourceEstimate.tstep); //100ms
290 // sourceEstimate = sourceEstimate.reduce(sample, 1);
291 
292  QList<Label> t_qListLabels;
293  QList<RowVector4i> t_qListRGBAs;
294 
295  //ToDo overload toLabels using instead of t_surfSet rr of MNESourceSpace
296  t_annotationSet.toLabels(t_surfSet, t_qListLabels, t_qListRGBAs);
297 
298 // InverseView view(minimumNorm.getSourceSpace(), t_qListLabels, t_qListRGBAs);
299 
300 // if (view.stereoType() != QGLView::RedCyanAnaglyph)
301 // view.camera()->setEyeSeparation(0.3f);
302 // QStringList args = QCoreApplication::arguments();
303 // int w_pos = args.indexOf("-width");
304 // int h_pos = args.indexOf("-height");
305 // if (w_pos >= 0 && h_pos >= 0)
306 // {
307 // bool ok = true;
308 // int w = args.at(w_pos + 1).toInt(&ok);
309 // if (!ok)
310 // {
311 // qWarning() << "Could not parse width argument:" << args;
312 // return 1;
313 // }
314 // int h = args.at(h_pos + 1).toInt(&ok);
315 // if (!ok)
316 // {
317 // qWarning() << "Could not parse height argument:" << args;
318 // return 1;
319 // }
320 // view.resize(w, h);
321 // }
322 // else
323 // {
324 // view.resize(800, 600);
325 // }
326 // view.show();
327 
328 // //Push Estimate
329 // view.pushSourceEstimate(sourceEstimate);
330 
331  if(!t_sFileNameStc.isEmpty())
332  {
333  QFile t_fileClusteredStc(t_sFileNameStc);
334  sourceEstimate.write(t_fileClusteredStc);
335  }
336 
337 //*/
338  return a.exec();//1;//a.exec();
339 }
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
evoked data
Definition: fiff_evoked.h:91
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.
Minimum norm estimation.
Definition: minimumnorm.h:82
AnnotationSet class declaration.
InverseView class declaration.
Label class declaration.
MNEMath class declaration.
A hemisphere set of surfaces.
Definition: surfaceset.h:83