MNE-CPP  beta 1.0
main.cpp
1 //=============================================================================================================
37 //*************************************************************************************************************
38 //=============================================================================================================
39 // INCLUDES
40 //=============================================================================================================
41 
42 #include <iostream>
43 #include <vector>
44 #include <math.h>
45 
46 
47 //*************************************************************************************************************
48 //=============================================================================================================
49 // MNE INCLUDES
50 //=============================================================================================================
51 
52 #include <fiff/fiff_evoked_set.h>
53 #include <mne/mne.h>
54 
55 
56 //*************************************************************************************************************
57 //=============================================================================================================
58 // QT INCLUDES
59 //=============================================================================================================
60 
61 #include <QtCore/QCoreApplication>
62 
63 
64 //*************************************************************************************************************
65 //=============================================================================================================
66 // USED NAMESPACES
67 //=============================================================================================================
68 
69 using namespace FIFFLIB;
70 using namespace MNELIB;
71 
72 
73 //*************************************************************************************************************
74 //=============================================================================================================
75 // MAIN
76 //=============================================================================================================
77 
78 //=============================================================================================================
87 int main(int argc, char *argv[])
88 {
89  QCoreApplication a(argc, argv);
90 
91 // fname_data - Name of the data file
92 // setno - Data set number
93 // fname_inv - Inverse operator file name
94 // nave - Number of averages (scales the noise covariance)
95 // If negative, the number of averages in the data will be
96 // used
97 // lambda2 - The regularization factor
98 // dSPM - do dSPM?
99 // sLORETA - do sLORETA?
100  QFile t_fileEvoked("./MNE-sample-data/MEG/sample/sample_audvis-ave.fif");
101  QFile t_fileInv("./MNE-sample-data/MEG/sample/sample_audvis-meg-eeg-oct-6-meg-eeg-inv.fif");
102  qint32 nave = 40;
103  float snr = 3.0f;
104  float lambda2 = pow(1.0f / snr, 2.0f);
105  bool dSPM = false;
106  bool sLORETA = true;
107 
108  //
109  // Read the data first
110  //
111  FiffEvokedSet evokedSet(t_fileEvoked);
112 
113  //
114  // Then the inverse operator
115  //
116  MNEInverseOperator inv_raw(t_fileInv);
117 
118  //
119  // Iterate over found data sets
120  //
121  for(qint32 setno = 0; setno < evokedSet.evoked.size(); ++setno)
122  {
123  printf(">> Computing inverse for %s data set <<\n", evokedSet.evoked[setno].comment.toLatin1().constData());
124  //
125  // Set up the inverse according to the parameters
126  //
127  if (nave < 0)
128  nave = evokedSet.evoked[setno].nave;
129 
130  MNEInverseOperator inv = inv_raw.prepare_inverse_operator(nave,lambda2,dSPM,sLORETA);
131  //
132  // Pick the correct channels from the data
133  //
134  FiffEvokedSet newEvokedSet = evokedSet.pick_channels(inv.noise_cov->names);
135 
136  evokedSet = newEvokedSet;
137 
138  printf("Picked %d channels from the data\n",evokedSet.info.nchan);
139  printf("Computing inverse...");
140  //
141  // Simple matrix multiplication followed by combination of the
142  // three current components
143  //
144  // This does all the data transformations to compute the weights for the
145  // eigenleads
146  //
147  SparseMatrix<double> reginv(inv.reginv.rows(),inv.reginv.rows());
148  // put this in the MNE algorithm class derived from inverse algorithm
149  //ToDo put this into a function of inv data
150  qint32 i;
151  for(i = 0; i < inv.reginv.rows(); ++i)
152  reginv.insert(i,i) = inv.reginv(i,0);
153 
154  MatrixXd trans = reginv*inv.eigen_fields->data*inv.whitener*inv.proj*evokedSet.evoked[setno].data;
155  //
156  // Transformation into current distributions by weighting the eigenleads
157  // with the weights computed above
158  //
159  MatrixXd sol;
160  if (inv.eigen_leads_weighted)
161  {
162  //
163  // R^0.5 has been already factored in
164  //
165  printf("(eigenleads already weighted)...");
166  sol = inv.eigen_leads->data*trans;
167  }
168  else
169  {
170  //
171  // R^0.5 has to factored in
172  //
173  printf("(eigenleads need to be weighted)...");
174 
175  SparseMatrix<double> sourceCov(inv.source_cov->data.rows(),inv.source_cov->data.rows());
176  for(i = 0; i < inv.source_cov->data.rows(); ++i)
177  sourceCov.insert(i,i) = sqrt(inv.source_cov->data(i,0));
178 
179  sol = sourceCov*inv.eigen_leads->data*trans;
180  }
181 
182  if (inv.source_ori == FIFFV_MNE_FREE_ORI)
183  {
184  printf("combining the current components...");
185  MatrixXd sol1(sol.rows()/3,sol.cols());
186  for(i = 0; i < sol.cols(); ++i)
187  {
188  VectorXd* tmp = MNE::combine_xyz(sol.block(0,i,sol.rows(),1));
189  sol1.block(0,i,sol.rows()/3,1) = tmp->cwiseSqrt();
190  delete tmp;
191  }
192  sol.resize(sol1.rows(),sol1.cols());
193  sol = sol1;
194  }
195  if (dSPM)
196  {
197  printf("(dSPM)...");
198  sol = inv.noisenorm*sol;
199  }
200  else if (sLORETA)
201  {
202  printf("(sLORETA)...");
203  sol = inv.noisenorm*sol;
204  }
205  printf("[done]\n");
206 
207  //Results
208  float tmin = ((float)evokedSet.evoked[setno].first) / evokedSet.info.sfreq;
209  float tstep = 1/evokedSet.info.sfreq;
210 
211  std::cout << "\npart ( block( 0, 0, 10, 10) ) of the inverse solution:\n" << sol.block(0,0,10,10) << std::endl;
212  printf("tmin = %f s\n", tmin);
213  printf("tstep = %f s\n", tstep);
214  }
215 
216  return a.exec();
217 }
FiffEvokedSet pick_channels(const QStringList &include=defaultQStringList, const QStringList &exclude=defaultQStringList) const
FiffNamedMatrix::SDPtr eigen_fields
SparseMatrix< double > noisenorm
FiffNamedMatrix::SDPtr eigen_leads
MNEInverseOperator prepare_inverse_operator(qint32 nave, float lambda2, bool dSPM, bool sLORETA=false) const
Definition: fiff.h:98
FiffEvokedSet class declaration.
evoked data set