45 #include <fiff/fiff_evoked.h>
66 #include <QGuiApplication>
76 using namespace FSLIB;
97 int main(
int argc,
char *argv[])
99 QGuiApplication a(argc, argv);
111 QFile t_fileRaw(
"E:/Data/sl_data/MEG/mind006/mind006_051209_median01_raw.fif");
112 QFile t_fileFwd(
"E:/Data/sl_data/MEG/mind006/mind006_051209_median01_raw-oct-6-fwd.fif");
113 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");
114 SurfaceSet t_surfSet(
"E:/Data/sl_data/subjects/mind006/surf/lh.white",
"E:/Data/sl_data/subjects/mind006/surf/rh.white");
116 QString t_sFileNameStc(
"");
118 qint32 numDipolePairs = 7;
120 qint32 samplesStcWindow = 100;
121 float stcOverlap = 0.0f;
123 qint32 startSample = 0;
124 qint32 numSample = 100000;
126 bool in_samples =
true;
127 bool keep_comp =
true;
130 for(qint32 i = 0; i < argc; ++i)
132 if(strcmp(argv[i],
"-stc") == 0 || strcmp(argv[i],
"--stc") == 0)
135 t_sFileNameStc = QString::fromUtf8(argv[i+1]);
151 float from = raw.first_samp + startSample;
152 float to = from + numSample - 1;
159 bool want_meg =
true;
160 bool want_eeg =
false;
161 bool want_stim =
false;
163 RowVectorXi picks = raw.info.pick_types(want_meg, want_eeg, want_stim, include);
169 if (raw.info.projs.size() == 0)
170 printf(
"No projector specified for these data\n");
176 for (k = 0; k < raw.info.projs.size(); ++k)
177 raw.info.projs[k].active =
true;
179 printf(
"%d projection items activated\n",raw.info.projs.size());
183 fiff_int_t nproj = raw.info.make_projector(raw.proj);
186 printf(
"The projection vectors do not apply to these channels\n");
188 printf(
"Created an SSP operator (subspace dimension = %d)\n",nproj);
194 qint32 current_comp = raw.info.get_current_comp();
195 qint32 dest_comp = -1;
197 if (current_comp > 0)
198 printf(
"Current compensation grade : %d\n",current_comp);
201 dest_comp = current_comp;
203 if (current_comp != dest_comp)
205 qDebug() <<
"This part needs to be debugged";
208 raw.info.set_current_comp(dest_comp);
209 printf(
"Appropriate compensator added to change to grade %d.\n",dest_comp);
213 printf(
"Could not make the compensator\n");
221 bool readSuccessful =
false;
225 readSuccessful = raw.read_raw_segment(data, times, (qint32)from, (qint32)to, picks);
227 readSuccessful = raw.read_raw_segment_times(data, times, from, to, picks);
231 printf(
"Could not read raw segment.\n");
235 printf(
"Read %d samples.\n",(qint32)data.cols());
248 PwlRapMusic t_pwlRapMusic(t_clusteredFwd,
false, numDipolePairs);
252 float tstep = 1.0f/raw.info.sfreq;
263 sourceEstimate.
data = MatrixXd::Zero(t_clusteredFwd.nsource, data.cols());
266 sourceEstimate.
vertices = VectorXi(t_clusteredFwd.src[0].vertno.size() + t_clusteredFwd.src[1].vertno.size());
267 sourceEstimate.
vertices << t_clusteredFwd.src[0].vertno, t_clusteredFwd.src[1].vertno;
269 sourceEstimate.
times = RowVectorXf::Zero(data.cols());
270 sourceEstimate.
times[0] = tmin;
271 for(qint32 i = 1; i < sourceEstimate.
times.size(); ++i)
272 sourceEstimate.
times[i] = sourceEstimate.
times[i-1] + tstep;
273 sourceEstimate.
tmin = tmin;
274 sourceEstimate.
tstep = tstep;
280 qint32 t_iNumSensors = data.rows();
281 qint32 t_iNumSteps = data.cols();
283 qint32 t_iSamplesOverlap = (qint32)floor(((
float)samplesStcWindow)*stcOverlap);
284 qint32 t_iSamplesDiscard = t_iSamplesOverlap/2;
286 MatrixXd measData = MatrixXd::Zero(t_iNumSensors, samplesStcWindow);
288 qint32 curSample = 0;
289 qint32 curResultSample = 0;
290 qint32 stcWindowSize = samplesStcWindow - 2*t_iSamplesDiscard;
295 if(curSample + samplesStcWindow >= t_iNumSteps)
298 measData = data.block(0, data.cols()-samplesStcWindow, t_iNumSensors, samplesStcWindow);
301 measData = data.block(0, curSample, t_iNumSensors, samplesStcWindow);
304 curSample += (samplesStcWindow - t_iSamplesOverlap);
306 curSample -= t_iSamplesDiscard;
309 MNESourceEstimate stcData = t_pwlRapMusic.calculateInverse(measData, 0.0f, tstep);
313 stcWindowSize = stcData.
data.cols() - curResultSample;
315 sourceEstimate.
data.block(0,curResultSample,sourceEstimate.
data.rows(),stcWindowSize) =
316 stcData.
data.block(0,0,stcData.
data.rows(),stcWindowSize);
318 curResultSample += stcWindowSize;
348 QList<Label> t_qListLabels;
349 QList<RowVector4i> t_qListRGBAs;
352 t_annotationSet.toLabels(t_surfSet, t_qListLabels, t_qListRGBAs);
354 InverseView view(t_pwlRapMusic.getSourceSpace(), t_qListLabels, t_qListRGBAs, 24,
true,
false,
false);
356 if (view.stereoType() != QGLView::RedCyanAnaglyph)
357 view.camera()->setEyeSeparation(0.3f);
358 QStringList args = QCoreApplication::arguments();
359 int w_pos = args.indexOf(
"-width");
360 int h_pos = args.indexOf(
"-height");
361 if (w_pos >= 0 && h_pos >= 0)
364 int w = args.at(w_pos + 1).toInt(&ok);
367 qWarning() <<
"Could not parse width argument:" << args;
370 int h = args.at(h_pos + 1).toInt(&ok);
373 qWarning() <<
"Could not parse height argument:" << args;
380 view.resize(800, 600);
385 view.pushSourceEstimate(sourceEstimate);
387 if(!t_sFileNameStc.isEmpty())
389 QFile t_fileClusteredStc(t_sFileNameStc);
390 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.
FIFF raw measurement data.
PwlRapMusic algorithm class declaration.
bool write(QIODevice &p_IODevice)
Surface class declaration.
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.