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);
119 QFile t_fileRaw(
"D:/Users/Christoph/SkyDrive/Thesis_Data/MIND/mind006_051209_auditory01_raw.fif");
120 QString t_sEventName =
"D:/Users/Christoph/SkyDrive/Thesis_Data/MIND/mind006_051209_auditory01_raw-eve.fif";
121 QFile t_fileFwd(
"D:/Users/Christoph/SkyDrive/Thesis_Data/MIND/mind006_051209_auditory01_raw-oct-6-fwd.fif");
122 AnnotationSet t_annotationSet(
"mind006", 2,
"aparc.a2009s",
"D:/Users/Christoph/SkyDrive/Thesis_Data/subjects");
123 SurfaceSet t_surfSet(
"mind006", 2,
"white",
"D:/Users/Christoph/SkyDrive/Thesis_Data/subjects");
131 QString t_sFileNameStc(
"");
134 bool doMovie =
false;
136 qint32 numDipolePairs = 1;
149 bool keep_comp =
false;
150 fiff_int_t dest_comp = 0;
151 bool pick_all =
true;
157 for(qint32 i = 0; i < argc; ++i)
159 if(strcmp(argv[i],
"-stc") == 0 || strcmp(argv[i],
"--stc") == 0)
162 t_sFileNameStc = QString::fromUtf8(argv[i+1]);
184 picks.resize(raw.info.nchan);
186 for(k = 0; k < raw.info.nchan; ++k)
193 include <<
"STI 014";
194 bool want_meg =
true;
195 bool want_eeg =
false;
196 bool want_stim =
false;
198 picks = raw.info.pick_types(want_meg, want_eeg, want_stim, include, raw.info.bads);
201 QStringList ch_names;
202 for(k = 0; k < picks.cols(); ++k)
203 ch_names << raw.info.ch_names[picks(0,k)];
208 if (raw.info.projs.size() == 0)
209 printf(
"No projector specified for these data\n");
215 for (k = 0; k < raw.info.projs.size(); ++k)
216 raw.info.projs[k].active =
true;
218 printf(
"%d projection items activated\n",raw.info.projs.size());
223 fiff_int_t nproj = raw.info.make_projector(raw.proj);
227 printf(
"The projection vectors do not apply to these channels\n");
231 printf(
"Created an SSP operator (subspace dimension = %d)\n",nproj);
238 qint32 current_comp = raw.info.get_current_comp();
239 if (current_comp > 0)
240 printf(
"Current compensation grade : %d\n",current_comp);
243 dest_comp = current_comp;
245 if (current_comp != dest_comp)
247 qDebug() <<
"This part needs to be debugged";
250 raw.info.set_current_comp(dest_comp);
251 printf(
"Appropriate compensator added to change to grade %d.\n",dest_comp);
255 printf(
"Could not make the compensator\n");
264 if (t_sEventName.size() == 0)
266 p = t_fileRaw.fileName().indexOf(
".fif");
269 t_sEventName = t_fileRaw.fileName().replace(p, 4,
"-eve.fif");
273 printf(
"Raw file name does not end properly\n");
278 t_EventFile.setFileName(t_sEventName);
280 printf(
"Events read from %s\n",t_sEventName.toUtf8().constData());
287 p = t_fileRaw.fileName().indexOf(
".fif");
290 t_EventFile.setFileName(t_sEventName);
293 printf(
"Error while read events.\n");
296 printf(
"Binary event file %s read\n",t_sEventName.toUtf8().constData());
303 printf(
"Text file %s is not supported jet.\n",t_sEventName.toUtf8().constData());
346 MatrixXi selected = MatrixXi::Zero(1, events.rows());
347 for (p = 0; p < events.rows(); ++p)
349 if (events(p,1) == 0 && events(p,2) == event)
351 selected(0,count) = p;
355 selected.conservativeResize(1, count);
357 printf(
"%d matching events found\n",count);
360 printf(
"No desired events found.\n");
365 fiff_int_t event_samp, from, to;
374 for (p = 0; p < count; ++p)
379 event_samp = events(selected(p),0);
380 from = event_samp + tmin*raw.info.sfreq;
381 to = event_samp + floor(tmax*raw.info.sfreq + 0.5);
385 if(raw.read_raw_segment(epoch->
epoch, timesDummy, from, to, picks))
389 times.resize(1, to-from+1);
390 for (qint32 i = 0; i < times.cols(); ++i)
391 times(0, i) = ((float)(from-event_samp+i)) / raw.info.sfreq;
394 epoch->
event = event;
395 epoch->
tmin = ((float)(from)-(float)(raw.first_samp))/raw.info.sfreq;
396 epoch->
tmax = ((
float)(to)-(
float)(raw.first_samp))/raw.info.sfreq;
402 printf(
"Can't read the event data segments");
409 printf(
"Read %d epochs, %d samples each.\n",data.size(),(qint32)data[0]->epoch.cols());
464 for(qint32 i = 0; i < vecSel.size(); ++i)
466 qint32 val = rand() % count;
471 std::cout <<
"Select following epochs to average:\n" << vecSel << std::endl;
473 FiffEvoked evoked = data.
average(raw.info, tmin*raw.info.sfreq, floor(tmax*raw.info.sfreq + 0.5), vecSel);
491 PwlRapMusic t_pwlRapMusic(t_clusteredFwd,
false, numDipolePairs);
494 t_pwlRapMusic.setStcAttr(200,0.5);
499 QList<qint64> qVecElapsedTime;
500 for(qint32 i = 0; i < 100; ++i)
505 sourceEstimate = t_pwlRapMusic.calculateInverse(pickedEvoked);
506 qVecElapsedTime.append(timer.elapsed());
509 double meanTime = 0.0;
512 for(qint32 i = offset; i < qVecElapsedTime.size(); ++i)
514 meanTime += qVecElapsedTime[i];
518 meanTime /= (double)c;
521 for(qint32 i = offset; i < qVecElapsedTime.size(); ++i)
522 varTime += pow(qVecElapsedTime[i] - meanTime,2);
524 varTime /= (double)c - 1.0f;
525 varTime = sqrt(varTime);
527 qDebug() <<
"RAP-MUSIC calculation took" << meanTime <<
"+-" << varTime <<
"ms in average";
559 QList<Label> t_qListLabels;
560 QList<RowVector4i> t_qListRGBAs;
563 t_annotationSet.toLabels(t_surfSet, t_qListLabels, t_qListRGBAs);
565 InverseView view(t_pwlRapMusic.getSourceSpace(), t_qListLabels, t_qListRGBAs, 24,
true,
false,
false);
567 if (view.stereoType() != QGLView::RedCyanAnaglyph)
568 view.camera()->setEyeSeparation(0.3f);
569 QStringList args = QCoreApplication::arguments();
570 int w_pos = args.indexOf(
"-width");
571 int h_pos = args.indexOf(
"-height");
572 if (w_pos >= 0 && h_pos >= 0)
575 int w = args.at(w_pos + 1).toInt(&ok);
578 qWarning() <<
"Could not parse width argument:" << args;
581 int h = args.at(h_pos + 1).toInt(&ok);
584 qWarning() <<
"Could not parse height argument:" << args;
591 view.resize(800, 600);
596 view.pushSourceEstimate(sourceEstimate);
598 if(!t_sFileNameStc.isEmpty())
600 QFile t_fileClusteredStc(t_sFileNameStc);
601 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)
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.
PwlRapMusic algorithm class declaration.
bool write(QIODevice &p_IODevice)
FiffEvoked pick_channels(const QStringList &include=defaultQStringList, const QStringList &exclude=defaultQStringList) const
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.
The PwlRapMusic class provides the POWELL RAP MUSIC Algorithm CPU implementation. ToDo: Paper referen...
MNEMath class declaration.
FIFF class declaration, which provides static wrapper functions to stay consistent with mne matlab to...
A hemisphere set of surfaces.