46 #include <fiff/fiff_evoked.h>
67 #include <QGuiApplication>
69 #include <QElapsedTimer>
80 using namespace FSLIB;
101 int main(
int argc,
char *argv[])
103 QGuiApplication a(argc, argv);
112 QFile t_fileFwd(
"D:/Data/MEG/mind006/mind006_051209_auditory01_raw-oct-6p-fwd.fif");
113 QFile t_fileCov(
"D:/Data/MEG/mind006/mind006_051209_auditory01_raw-cov.fif");
114 QFile t_fileRaw(
"D:/Data/MEG/mind006/mind006_051209_auditory01_raw.fif");
115 QString t_sEventName =
"D:/Data/MEG/mind006/mind006_051209_auditory01_raw-eve.fif";
116 AnnotationSet t_annotationSet(
"mind006", 2,
"aparc.a2009s",
"D:/Data/subjects");
117 SurfaceSet t_surfSet(
"mind006", 2,
"white",
"D:/Data/subjects");
125 bool keep_comp =
false;
126 fiff_int_t dest_comp = 0;
127 bool pick_all =
true;
142 picks.resize(raw.info.nchan);
144 for(k = 0; k < raw.info.nchan; ++k)
151 include <<
"STI 014";
152 bool want_meg =
true;
153 bool want_eeg =
false;
154 bool want_stim =
false;
157 picks = raw.info.pick_types(want_meg, want_eeg, want_stim, include, raw.info.bads);
160 QStringList ch_names;
161 for(k = 0; k < picks.cols(); ++k)
162 ch_names << raw.info.ch_names[picks(0,k)];
167 if (raw.info.projs.size() == 0)
168 printf(
"No projector specified for these data\n");
174 for (k = 0; k < raw.info.projs.size(); ++k)
175 raw.info.projs[k].active =
true;
177 printf(
"%d projection items activated\n",raw.info.projs.size());
182 fiff_int_t nproj = raw.info.make_projector(raw.proj);
186 printf(
"The projection vectors do not apply to these channels\n");
190 printf(
"Created an SSP operator (subspace dimension = %d)\n",nproj);
198 qint32 current_comp = raw.info.get_current_comp();
199 if (current_comp > 0)
200 printf(
"Current compensation grade : %d\n",current_comp);
203 dest_comp = current_comp;
205 if (current_comp != dest_comp)
207 qDebug() <<
"This part needs to be debugged";
211 raw.info.set_current_comp(dest_comp);
212 printf(
"Appropriate compensator added to change to grade %d.\n",dest_comp);
216 printf(
"Could not make the compensator\n");
225 if (t_sEventName.size() == 0)
227 p = t_fileRaw.fileName().indexOf(
".fif");
230 t_sEventName = t_fileRaw.fileName().replace(p, 4,
"-eve.fif");
234 printf(
"Raw file name does not end properly\n");
239 t_EventFile.setFileName(t_sEventName);
241 printf(
"Events read from %s\n",t_sEventName.toUtf8().constData());
248 p = t_fileRaw.fileName().indexOf(
".fif");
251 t_EventFile.setFileName(t_sEventName);
254 printf(
"Error while read events.\n");
257 printf(
"Binary event file %s read\n",t_sEventName.toUtf8().constData());
264 printf(
"Text file %s is not supported jet.\n",t_sEventName.toUtf8().constData());
307 MatrixXi selected = MatrixXi::Zero(1, events.rows());
308 for (p = 0; p < events.rows(); ++p)
310 if (events(p,1) == 0 && events(p,2) == event)
312 selected(0,count) = p;
316 selected.conservativeResize(1, count);
318 printf(
"%d matching events found\n",count);
321 printf(
"No desired events found.\n");
326 fiff_int_t event_samp, from, to;
335 for (p = 0; p < count; ++p)
340 event_samp = events(selected(p),0);
341 from = event_samp + tmin*raw.info.sfreq;
342 to = event_samp + floor(tmax*raw.info.sfreq + 0.5);
346 if(raw.read_raw_segment(epoch->
epoch, timesDummy, from, to, picks))
350 times.resize(1, to-from+1);
351 for (qint32 i = 0; i < times.cols(); ++i)
352 times(0, i) = ((float)(from-event_samp+i)) / raw.info.sfreq;
355 epoch->
event = event;
356 epoch->
tmin = ((float)(from)-(float)(raw.first_samp))/raw.info.sfreq;
357 epoch->
tmax = ((
float)(to)-(
float)(raw.first_samp))/raw.info.sfreq;
363 printf(
"Can't read the event data segments");
370 printf(
"Read %d epochs, %d samples each.\n",data.size(),(qint32)data[0]->epoch.cols());
373 std::cout << data[0]->epoch.block(0,0,10,10) << std::endl;
374 qDebug() << data[0]->epoch.rows() <<
" x " << data[0]->epoch.cols();
376 std::cout << times.block(0,0,1,10) << std::endl;
377 qDebug() << times.rows() <<
" x " << times.cols();
404 vecSel << 0, 96, 80, 55, 66, 25, 26, 2, 55, 58, 6, 88;
407 std::cout <<
"Select following epochs to average:\n" << vecSel << std::endl;
409 FiffEvoked evoked = data.
average(raw.info, tmin*raw.info.sfreq, floor(tmax*raw.info.sfreq + 0.5), vecSel);
417 QString method(
"dSPM");
419 QString t_sFileNameClusteredInv(
"");
420 QString t_sFileNameStc(
"test_mind006_051209_auditory01.stc");
423 for(qint32 i = 0; i < argc; ++i)
425 if(strcmp(argv[i],
"-snr") == 0 || strcmp(argv[i],
"--snr") == 0)
428 snr = atof(argv[i+1]);
430 else if(strcmp(argv[i],
"-method") == 0 || strcmp(argv[i],
"--method") == 0)
433 method = QString::fromUtf8(argv[i+1]);
435 else if(strcmp(argv[i],
"-inv") == 0 || strcmp(argv[i],
"--inv") == 0)
438 t_sFileNameClusteredInv = QString::fromUtf8(argv[i+1]);
440 else if(strcmp(argv[i],
"-stc") == 0 || strcmp(argv[i],
"--stc") == 0)
443 t_sFileNameStc = QString::fromUtf8(argv[i+1]);
447 double lambda2 = 1.0 / pow(snr, 2);
448 qDebug() <<
"Start calculation with: SNR" << snr <<
"; Lambda" << lambda2 <<
"; Method" << method <<
"; stc:" << t_sFileNameStc;
465 noise_cov = noise_cov.regularize(evoked.
info, 0.05, 0.05, 0.1,
true);
483 if(!t_sFileNameClusteredInv.isEmpty())
485 QFile t_fileClusteredInverse(t_sFileNameClusteredInv);
486 inverse_operator.write(t_fileClusteredInverse);
492 MinimumNorm minimumNorm(inverse_operator, lambda2, method);
498 minimumNorm.doInverseSetup(vecSel.size(),
false);
501 QList<qint64> qVecElapsedTime;
502 for(qint32 i = 0; i < 100; ++i)
507 sourceEstimate = minimumNorm.calculateInverse(evoked.
data, evoked.
times(0), evoked.
times(1)-evoked.
times(0));
508 qVecElapsedTime.append(timer.elapsed());
511 double meanTime = 0.0;
514 for(qint32 i = offset; i < qVecElapsedTime.size(); ++i)
516 meanTime += qVecElapsedTime[i];
520 meanTime /= (double)c;
523 for(qint32 i = offset; i < qVecElapsedTime.size(); ++i)
524 varTime += pow(qVecElapsedTime[i] - meanTime,2);
526 varTime /= (double)c - 1.0f;
527 varTime = sqrt(varTime);
529 qDebug() <<
"MNE calculation took" << meanTime <<
"+-" << varTime <<
"ms in average";
539 std::cout <<
"\nsourceEstimate:\n" << sourceEstimate.
data.block(0,0,10,10) << std::endl;
540 std::cout <<
"time\n" << sourceEstimate.
times.block(0,0,1,10) << std::endl;
541 std::cout <<
"timeMin\n" << sourceEstimate.
times[0] << std::endl;
542 std::cout <<
"timeMax\n" << sourceEstimate.
times[sourceEstimate.
times.size()-1] << std::endl;
543 std::cout <<
"time step\n" << sourceEstimate.
tstep << std::endl;
586 std::cout <<
"Condition Number:\n" << t_dConditionNumber << std::endl;
587 std::cout <<
"Clustered Condition Number:\n" << t_dConditionNumberClustered << std::endl;
589 std::cout <<
"ForwardSolution" << t_Fwd.sol->data.block(0,0,10,10) << std::endl;
591 std::cout <<
"Clustered ForwardSolution" << t_clusteredFwd.
sol->data.block(0,0,10,10) << std::endl;
623 QList<Label> t_qListLabels;
624 QList<RowVector4i> t_qListRGBAs;
627 t_annotationSet.toLabels(t_surfSet, t_qListLabels, t_qListRGBAs);
629 InverseView view(minimumNorm.getSourceSpace(), t_qListLabels, t_qListRGBAs, 24,
true,
false,
true);
631 if (view.stereoType() != QGLView::RedCyanAnaglyph)
632 view.camera()->setEyeSeparation(0.3f);
633 QStringList args = QCoreApplication::arguments();
634 int w_pos = args.indexOf(
"-width");
635 int h_pos = args.indexOf(
"-height");
636 if (w_pos >= 0 && h_pos >= 0)
639 int w = args.at(w_pos + 1).toInt(&ok);
642 qWarning() <<
"Could not parse width argument:" << args;
645 int h = args.at(h_pos + 1).toInt(&ok);
648 qWarning() <<
"Could not parse height argument:" << args;
655 view.resize(800, 600);
660 view.pushSourceEstimate(sourceEstimate);
662 if(!t_sFileNameStc.isEmpty())
664 QFile t_fileClusteredStc(t_sFileNameStc);
665 sourceEstimate.
write(t_fileClusteredStc);
static bool make_compensator(const FiffInfo &info, fiff_int_t from, fiff_int_t to, FiffCtfComp &ctf_comp, bool exclude_comp_chs=false)
FIFF measurement file information.
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.
static bool read_events(QIODevice &p_IODevice, MatrixXi &eventlist)
FIFF raw measurement data.
static double getConditionNumber(const MatrixXd &A, VectorXd &s)
Minimum norm class declaration.
bool write(QIODevice &p_IODevice)
FiffNamedMatrix::SDPtr sol
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
AnnotationSet class declaration.
InverseView 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.