45 #include <fiff/fiff_evoked.h>
64 #include <QGuiApplication>
74 using namespace FSLIB;
95 int main(
int argc,
char *argv[])
97 QGuiApplication a(argc, argv);
126 QFile t_fileRaw(
"D:/Dropbox/Masterarbeit DB/Messdaten/EEG/2014_02_04_Lorenz_Esch_005/Processed/averaged/EEG_data_002_voluntary_right_tapping_filtered_07_40_Averaged_-1_1_184.fif");
127 QFile t_fileFwd(
"D:/Dropbox/Masterarbeit DB/Messdaten/Forward solutions/Lorenz-140128-Duke128-fwd.fif");
128 AnnotationSet t_annotationSet(
"D:/Dropbox/Masterarbeit DB/Messdaten/Forward solutions/atlas/lh.aparc.a2009s.annot",
"D:/Dropbox/Masterarbeit DB/Messdaten/Forward solutions/atlas/rh.aparc.a2009s.annot");
129 SurfaceSet t_surfSet(
"D:/Dropbox/Masterarbeit DB/Messdaten/Forward solutions/surface/lh.white",
"D:/Dropbox/Masterarbeit DB/Messdaten/Forward solutions/surface/rh.white");
130 QString t_sFileNameStc(
"D:/Dropbox/Masterarbeit DB/Messdaten/EEG/2014_02_04_Lorenz_Esch_005/Processed/stc/EEG_data_002_voluntary_right_tapping_filtered_07_40_Averaged_-1_1_184_RAP.stc");
138 qint32 numDipolePairs = 7;
140 qint32 samplesStcWindow = 300;
141 float stcOverlap = 0.0f;
143 qint32 startSample = 0;
144 qint32 numSample = 40000;
146 bool in_samples =
true;
147 bool keep_comp =
true;
150 for(qint32 i = 0; i < argc; ++i)
152 if(strcmp(argv[i],
"-stc") == 0 || strcmp(argv[i],
"--stc") == 0)
155 t_sFileNameStc = QString::fromUtf8(argv[i+1]);
166 QList<Label> t_qListLabels;
167 QList<RowVector4i> t_qListRGBAs;
170 t_annotationSet.toLabels(t_surfSet, t_qListLabels, t_qListRGBAs);
172 QList<Label> t_qListLabelSelection;
175 t_qListLabelSelection << t_qListLabels[28] << t_qListLabels[29] << t_qListLabels[45] << t_qListLabels[66];
177 t_qListLabelSelection << t_qListLabels[103] << t_qListLabels[104] << t_qListLabels[120]<< t_qListLabels[141];
194 float from = raw.first_samp + startSample;
195 float to = from + numSample - 1;
202 bool want_meg =
false;
203 bool want_eeg =
true;
204 bool want_stim =
false;
206 RowVectorXi picks = raw.info.pick_types(want_meg, want_eeg, want_stim, include);
212 if (raw.info.projs.size() == 0)
213 printf(
"No projector specified for these data\n");
219 for (k = 0; k < raw.info.projs.size(); ++k)
220 raw.info.projs[k].active =
true;
222 printf(
"%d projection items activated\n",raw.info.projs.size());
226 fiff_int_t nproj = raw.info.make_projector(raw.proj);
229 printf(
"The projection vectors do not apply to these channels\n");
231 printf(
"Created an SSP operator (subspace dimension = %d)\n",nproj);
237 qint32 current_comp = raw.info.get_current_comp();
238 qint32 dest_comp = -1;
240 if (current_comp > 0)
241 printf(
"Current compensation grade : %d\n",current_comp);
244 dest_comp = current_comp;
246 if (current_comp != dest_comp)
248 qDebug() <<
"This part needs to be debugged";
251 raw.info.set_current_comp(dest_comp);
252 printf(
"Appropriate compensator added to change to grade %d.\n",dest_comp);
256 printf(
"Could not make the compensator\n");
264 bool readSuccessful =
false;
268 readSuccessful = raw.read_raw_segment(data, times, (qint32)from, (qint32)to, picks);
270 readSuccessful = raw.read_raw_segment_times(data, times, from, to, picks);
274 printf(
"Could not read raw segment.\n");
278 printf(
"Read %d samples.\n",(qint32)data.cols());
292 PwlRapMusic t_pwlRapMusic(t_clusteredFwd,
false, numDipolePairs);
296 float tstep = 1.0f/raw.info.sfreq;
307 sourceEstimate.
data = MatrixXd::Zero(t_clusteredFwd.nsource, data.cols());
310 sourceEstimate.
vertices = VectorXi(t_clusteredFwd.src[0].vertno.size() + t_clusteredFwd.src[1].vertno.size());
311 sourceEstimate.
vertices << t_clusteredFwd.src[0].vertno, t_clusteredFwd.src[1].vertno;
313 sourceEstimate.
times = RowVectorXf::Zero(data.cols());
314 sourceEstimate.
times[0] = tmin;
315 for(qint32 i = 1; i < sourceEstimate.
times.size(); ++i)
316 sourceEstimate.
times[i] = sourceEstimate.
times[i-1] + tstep;
317 sourceEstimate.
tmin = tmin;
318 sourceEstimate.
tstep = tstep;
324 qint32 t_iNumSensors = data.rows();
325 qint32 t_iNumSteps = data.cols();
327 qint32 t_iSamplesOverlap = (qint32)floor(((
float)samplesStcWindow)*stcOverlap);
328 qint32 t_iSamplesDiscard = t_iSamplesOverlap/2;
330 MatrixXd measData = MatrixXd::Zero(t_iNumSensors, samplesStcWindow);
332 qint32 curSample = 0;
333 qint32 curResultSample = 0;
334 qint32 stcWindowSize = samplesStcWindow - 2*t_iSamplesDiscard;
339 if(curSample + samplesStcWindow >= t_iNumSteps)
342 measData = data.block(0, data.cols()-samplesStcWindow, t_iNumSensors, samplesStcWindow);
345 measData = data.block(0, curSample, t_iNumSensors, samplesStcWindow);
348 curSample += (samplesStcWindow - t_iSamplesOverlap);
350 curSample -= t_iSamplesDiscard;
353 MNESourceEstimate stcData = t_pwlRapMusic.calculateInverse(measData, 0.0f, tstep);
357 stcWindowSize = stcData.
data.cols() - curResultSample;
359 sourceEstimate.
data.block(0,curResultSample,sourceEstimate.
data.rows(),stcWindowSize) =
360 stcData.
data.block(0,0,stcData.
data.rows(),stcWindowSize);
362 curResultSample += stcWindowSize;
372 if(!t_sFileNameStc.isEmpty())
374 QFile t_fileClusteredStc(t_sFileNameStc);
375 sourceEstimate.
write(t_fileClusteredStc);
381 InverseView view(t_pwlRapMusic.getSourceSpace(), t_qListLabels, t_qListRGBAs, 24,
true,
false,
false);
383 if (view.stereoType() != QGLView::RedCyanAnaglyph)
384 view.camera()->setEyeSeparation(0.3f);
385 QStringList args = QCoreApplication::arguments();
386 int w_pos = args.indexOf(
"-width");
387 int h_pos = args.indexOf(
"-height");
388 if (w_pos >= 0 && h_pos >= 0)
391 int w = args.at(w_pos + 1).toInt(&ok);
394 qWarning() <<
"Could not parse width argument:" << args;
397 int h = args.at(h_pos + 1).toInt(&ok);
400 qWarning() <<
"Could not parse height argument:" << args;
407 view.resize(800, 600);
412 view.pushSourceEstimate(sourceEstimate);
414 if(!t_sFileNameStc.isEmpty())
416 QFile t_fileClusteredStc(t_sFileNameStc);
417 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
SurfaceSet class declaration.
FIFF raw measurement data.
PwlRapMusic algorithm class declaration.
bool write(QIODevice &p_IODevice)
MNEForwardSolution pick_regions(const QList< Label > &p_qListLabels) const
MNESourceEstimate class declaration.
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.