MNE-CPP 0.1.9
A Framework for Electrophysiology
Loading...
Searching...
No Matches
fiff_evoked_set.cpp
Go to the documentation of this file.
1//=============================================================================================================
37//=============================================================================================================
38// INCLUDES
39//=============================================================================================================
40
41#include "fiff_evoked_set.h"
42#include "fiff_tag.h"
43#include "fiff_dir_node.h"
44#include "fiff_stream.h"
45
46//=============================================================================================================
47// EIGEN INCLUDES
48//=============================================================================================================
49
50#include <Eigen/SparseCore>
51
52//=============================================================================================================
53// QT INCLUDES
54//=============================================================================================================
55
56#include <QFile>
57
58//=============================================================================================================
59// USED NAMESPACES
60//=============================================================================================================
61
62using namespace FIFFLIB;
63using namespace Eigen;
64
65//=============================================================================================================
66// DEFINE MEMBER METHODS
67//=============================================================================================================
68
70{
71 qRegisterMetaType<FIFFLIB::FiffEvokedSet>("FIFFLIB::FiffEvokedSet");
72 qRegisterMetaType<FIFFLIB::FiffEvokedSet::SPtr>("FIFFLIB::FiffEvokedSet::SPtr");
73}
74
75//=============================================================================================================
76
77FiffEvokedSet::FiffEvokedSet(QIODevice& p_IODevice)
78{
79 qRegisterMetaType<FIFFLIB::FiffEvokedSet>("FIFFLIB::FiffEvokedSet");
80 qRegisterMetaType<FIFFLIB::FiffEvokedSet::SPtr>("FIFFLIB::FiffEvokedSet::SPtr");
81
82 if(!FiffEvokedSet::read(p_IODevice, *this))
83 {
84 printf("\tFiff evoked data set not found.\n");//ToDo Throw here
85 return;
86 }
87}
88
89//=============================================================================================================
90
92: info(p_FiffEvokedSet.info)
93, evoked(p_FiffEvokedSet.evoked)
94{
95}
96
97//=============================================================================================================
98
102
103//=============================================================================================================
104
106{
107 info.clear();
108 evoked.clear();
109}
110
111//=============================================================================================================
112
114 const QStringList& exclude) const
115{
116 FiffEvokedSet res;
117 res.info = this->info;
118 QList<FiffEvoked>::ConstIterator ev;
119 for(ev = evoked.begin(); ev != evoked.end(); ++ev)
120 res.evoked.push_back(ev->pick_channels(include, exclude));
121
122 return res;
123
124 //### OLD MATLAB oriented implementation ###
125
126// if(include.size() == 0 && exclude.size() == 0)
127// return FiffEvokedDataSet(*this);
128
129// RowVectorXi sel = FiffInfo::pick_channels(this->info.ch_names, include, exclude);
130// if (sel.cols() == 0)
131// {
132// qWarning("Warning : No channels match the selection.\n");
133// return FiffEvokedDataSet(*this);
134// }
135
136// FiffEvokedDataSet res(*this);
137// //
138// // Modify the measurement info
139// //
142// res.info = FiffInfo(res.info.pick_info(sel));
143// //
144// // Create the reduced data set
145// //
146// MatrixXd selBlock(1,1);
147// qint32 k, l;
148// for(k = 0; k < res.evoked.size(); ++k)
149// {
150// if(selBlock.rows() != sel.cols() || selBlock.cols() != res.evoked[k]->epochs.cols())
151// selBlock.resize(sel.cols(), res.evoked[k]->epochs.cols());
152// for(l = 0; l < sel.cols(); ++l)
153// {
154// selBlock.block(l,0,1,selBlock.cols()) = res.evoked[k]->epochs.block(sel(0,l),0,1,selBlock.cols());
155// }
156// res.evoked[k]->epochs.resize(sel.cols(), res.evoked[k]->epochs.cols());
157// res.evoked[k]->epochs = selBlock;
158// }
159
160// return res;
161}
162
163//=============================================================================================================
164
166 fiff_int_t to) const
167{
168 qint32 now = p_FiffEvokedSet.info.get_current_comp();
169 FiffCtfComp ctf_comp;
170
171 if(now == to)
172 {
173 printf("Data is already compensated as desired.\n");
174 return false;
175 }
176
177 //Make the compensator and apply it to all data sets
178 p_FiffEvokedSet.info.make_compensator(now,to,ctf_comp);
179
180 for(qint16 i=0; i < p_FiffEvokedSet.evoked.size(); ++i)
181 {
182 p_FiffEvokedSet.evoked[i].data = ctf_comp.data->data*p_FiffEvokedSet.evoked[i].data;
183 }
184
185 //Update the compensation info in the channel descriptors
186 p_FiffEvokedSet.info.set_current_comp(to);
187
188 return true;
189}
190
191//=============================================================================================================
192
193bool FiffEvokedSet::find_evoked(const FiffEvokedSet& p_FiffEvokedSet) const
194{
195 if(!p_FiffEvokedSet.evoked.size()) {
196 printf("No evoked response data sets in %s\n",p_FiffEvokedSet.info.filename.toUtf8().constData());
197 return false;
198 }
199 else
200 printf("\nFound %lld evoked response data sets in %s :\n",p_FiffEvokedSet.evoked.size(),p_FiffEvokedSet.info.filename.toUtf8().constData());
201
202 for(qint32 i = 0; i < p_FiffEvokedSet.evoked.size(); ++i) {
203 printf("%s (%s)\n",p_FiffEvokedSet.evoked.at(i).comment.toUtf8().constData(),p_FiffEvokedSet.evoked.at(i).aspectKindToString().toUtf8().constBegin());
204 }
205
206 return true;
207}
208
209//=============================================================================================================
210
211bool FiffEvokedSet::read(QIODevice& p_IODevice,
212 FiffEvokedSet& p_FiffEvokedSet, QPair<float,float> baseline,
213 bool proj)
214{
215 p_FiffEvokedSet.clear();
216
217 //
218 // Open the file
219 //
220 FiffStream::SPtr t_pStream(new FiffStream(&p_IODevice));
221 QString t_sFileName = t_pStream->streamName();
222
223 printf("Exploring %s ...\n",t_sFileName.toUtf8().constData());
224
225 if(!t_pStream->open())
226 return false;
227 //
228 // Read the measurement info
229 //
231 if(!t_pStream->read_meas_info(t_pStream->dirtree(), p_FiffEvokedSet.info, meas))
232 return false;
233 p_FiffEvokedSet.info.filename = t_sFileName; //move fname storage to read_meas_info member function
234 //
235 // Locate the data of interest
236 //
237 QList<FiffDirNode::SPtr> processed = meas->dir_tree_find(FIFFB_PROCESSED_DATA);
238 if (processed.size() == 0)
239 {
240 qWarning("Could not find processed data");
241 return false;
242 }
243 //
244 QList<FiffDirNode::SPtr> evoked_node = meas->dir_tree_find(FIFFB_EVOKED);
245 if (evoked_node.size() == 0)
246 {
247 qWarning("Could not find evoked data");
248 return false;
249 }
250
251 QStringList comments;
252 QList<fiff_int_t> aspect_kinds;
253 QString t;
254 if(!t_pStream->get_evoked_entries(evoked_node, comments, aspect_kinds, t))
255 t = QString("None found, must use integer");
256 printf("\tFound %lld datasets\n", evoked_node.size());
257
258 for(qint32 i = 0; i < comments.size(); ++i)
259 {
260 QFile t_file(p_FiffEvokedSet.info.filename);
261 printf(">> Processing %s <<\n", comments[i].toUtf8().constData());
262 FiffEvoked t_FiffEvoked;
263 if(FiffEvoked::read(t_file, t_FiffEvoked, i, baseline, proj))
264 p_FiffEvokedSet.evoked.push_back(t_FiffEvoked);
265 }
266
267 return true;
268
269 //### OLD MATLAB oriented implementation ###
270
271// data.clear();
272
273// if (setno < 0)
274// {
275// printf("Data set selector must be positive\n");
276// return false;
277// }
278// //
279// // Open the file
280// //
281// FiffStream::SPtr t_pStream(new FiffStream(&p_IODevice));
282// QString t_sFileName = t_pStream->streamName();
283
284// printf("Reading %s ...\n",t_sFileName.toUtf8().constData());
285// FiffDirNode t_Tree;
286// QList<FiffDirEntry> t_Dir;
287
288// if(!t_pStream->open(t_Tree, t_Dir))
289// {
290// if(t_pStream)
291// delete t_pStream;
292// return false;
293// }
294// //
295// // Read the measurement info
296// //
297// FiffInfo info;// = NULL;
298// FiffDirNode meas;
299// if(!t_pStream->read_meas_info(t_Tree, info, meas))
300// return false;
301// info.filename = t_sFileName; //move fname storage to read_meas_info member function
302// //
303// // Locate the data of interest
304// //
305// QList<FiffDirNode> processed = meas.dir_tree_find(FIFFB_PROCESSED_DATA);
306// if (processed.size() == 0)
307// {
308// if(t_pStream)
309// delete t_pStream;
310// printf("Could not find processed data\n");
311// return false;
312// }
313// //
314// QList<FiffDirNode> evoked = meas.dir_tree_find(FIFFB_EVOKED);
315// if (evoked.size() == 0)
316// {
317// if(t_pStream)
318// delete t_pStream;
319// printf("Could not find evoked data");
320// return false;
321// }
322// //
323// // Identify the aspects
324// //
325// fiff_int_t naspect = 0;
326// fiff_int_t nsaspects = 0;
327// qint32 oldsize = 0;
328// MatrixXi is_smsh(1,0);
329// QList< QList<FiffDirNode> > sets_aspects;
330// QList< qint32 > sets_naspect;
331// QList<FiffDirNode> saspects;
332// qint32 k;
333// for (k = 0; k < evoked.size(); ++k)
334// {
337
338// sets_aspects.append(evoked[k].dir_tree_find(FIFFB_ASPECT));
339// sets_naspect.append(sets_aspects[k].size());
340
341// if (sets_naspect[k] > 0)
342// {
343// oldsize = is_smsh.cols();
344// is_smsh.conservativeResize(1, oldsize + sets_naspect[k]);
345// is_smsh.block(0, oldsize, 1, sets_naspect[k]) = MatrixXi::Zero(1, sets_naspect[k]);
346// naspect += sets_naspect[k];
347// }
348// saspects = evoked[k].dir_tree_find(FIFFB_SMSH_ASPECT);
349// nsaspects = saspects.size();
350// if (nsaspects > 0)
351// {
352// sets_naspect[k] += nsaspects;
353// sets_aspects[k].append(saspects);
354
355// oldsize = is_smsh.cols();
356// is_smsh.conservativeResize(1, oldsize + sets_naspect[k]);
357// is_smsh.block(0, oldsize, 1, sets_naspect[k]) = MatrixXi::Ones(1, sets_naspect[k]);
358// naspect += nsaspects;
359// }
360// }
361// printf("\t%d evoked data sets containing a total of %d data aspects in %s\n",evoked.size(),naspect,t_sFileName.toUtf8().constData());
362// if (setno >= naspect || setno < 0)
363// {
364// if(t_pStream)
365// delete t_pStream;
366// printf("Data set selector out of range\n");
367// return false;
368// }
369// //
370// // Next locate the evoked data set
371// //
372// qint32 p = 0;
373// qint32 a = 0;
374// bool goon = true;
375// FiffDirNode my_evoked;
376// FiffDirNode my_aspect;
377// for(k = 0; k < evoked.size(); ++k)
378// {
379// for (a = 0; a < sets_naspect[k]; ++a)
380// {
381// if(p == setno)
382// {
383// my_evoked = evoked[k];
384// my_aspect = sets_aspects[k][a];
385// goon = false;
386// break;
387// }
388// ++p;
389// }
390// if (!goon)
391// break;
392// }
393// //
394// // The desired data should have been found but better to check
395// //
396// if (my_evoked.isEmpty() || my_aspect.isEmpty())
397// {
398// if(t_pStream)
399// delete t_pStream;
400// printf("Desired data set not found\n");
401// return false;
402// }
403// //
404// // Now find the data in the evoked block
405// //
406// fiff_int_t nchan = 0;
407// float sfreq = -1.0f;
408// QList<FiffChInfo> chs;
409// fiff_int_t kind, pos, first, last;
410// FiffTag* t_pTag = NULL;
411// QString comment("");
412// for (k = 0; k < my_evoked.nent; ++k)
413// {
414// kind = my_evoked.dir[k].kind;
415// pos = my_evoked.dir[k].pos;
416// switch (kind)
417// {
418// case FIFF_COMMENT:
419// FiffTag::read_tag(t_pStream,t_pTag,pos);
420// comment = t_pTag->toString();
421// break;
422// case FIFF_FIRST_SAMPLE:
423// FiffTag::read_tag(t_pStream,t_pTag,pos);
424// first = *t_pTag->toInt();
425// break;
426// case FIFF_LAST_SAMPLE:
427// FiffTag::read_tag(t_pStream,t_pTag,pos);
428// last = *t_pTag->toInt();
429// break;
430// case FIFF_NCHAN:
431// FiffTag::read_tag(t_pStream,t_pTag,pos);
432// nchan = *t_pTag->toInt();
433// break;
434// case FIFF_SFREQ:
435// FiffTag::read_tag(t_pStream,t_pTag,pos);
436// sfreq = *t_pTag->toFloat();
437// break;
438// case FIFF_CH_INFO:
439// FiffTag::read_tag(t_pStream, t_pTag, pos);
440// chs.append( t_pTag->toChInfo() );
441// break;
442// }
443// }
444// if (comment.isEmpty())
445// comment = QString("No comment");
446// //
447// // Local channel information?
448// //
449// if (nchan > 0)
450// {
451// if (chs.size() == 0)
452// {
453// if(t_pStream)
454// delete t_pStream;
455// printf("Local channel information was not found when it was expected.\n");
456// return false;
457// }
458// if (chs.size() != nchan)
459// {
460// if(t_pStream)
461// delete t_pStream;
462// printf("Number of channels and number of channel definitions are different\n");
463// return false;
464// }
465// info.chs = chs;
466// info.nchan = nchan;
467// printf("\tFound channel information in evoked data. nchan = %d\n",nchan);
468// if (sfreq > 0.0f)
469// info.sfreq = sfreq;
470// }
471// qint32 nsamp = last-first+1;
472// printf("\tFound the data of interest:\n");
473// printf("\t\tt = %10.2f ... %10.2f ms (%s)\n", 1000*(float)first/info.sfreq, 1000*(float)last/info.sfreq,comment.toUtf8().constData());
474// if (info.comps.size() > 0)
475// printf("\t\t%d CTF compensation matrices available\n", info.comps.size());
476// //
477// // Read the data in the aspect block
478// //
479// fiff_int_t nepoch = 0;
480// fiff_int_t aspect_kind = -1;
481// fiff_int_t nave = -1;
482// QList<FiffTag*> epoch;
483// for (k = 0; k < my_aspect.nent; ++k)
484// {
485// kind = my_aspect.dir[k].kind;
486// pos = my_aspect.dir[k].pos;
487
488// switch (kind)
489// {
490// case FIFF_COMMENT:
491// FiffTag::read_tag(t_pStream, t_pTag, pos);
492// comment = t_pTag->toString();
493// break;
494// case FIFF_ASPECT_KIND:
495// FiffTag::read_tag(t_pStream, t_pTag, pos);
496// aspect_kind = *t_pTag->toInt();
497// break;
498// case FIFF_NAVE:
499// FiffTag::read_tag(t_pStream, t_pTag, pos);
500// nave = *t_pTag->toInt();
501// break;
502// case FIFF_EPOCH:
503// FiffTag::read_tag(t_pStream, t_pTag, pos);
504// epoch.append(new FiffTag(t_pTag));
505// ++nepoch;
506// break;
507// }
508// }
509// if (nave == -1)
510// nave = 1;
511// printf("\t\tnave = %d aspect type = %d\n", nave, aspect_kind);
512// if (nepoch != 1 && nepoch != info.nchan)
513// {
514// if(t_pStream)
515// delete t_pStream;
516// printf("Number of epoch tags is unreasonable (nepoch = %d nchan = %d)\n", nepoch, info.nchan);
517// return false;
518// }
519// //
520// MatrixXd all;// = NULL;
521// if (nepoch == 1)
522// {
523// //
524// // Only one epoch
525// //
526// all = epoch[0]->toFloatMatrix();
527// all.transposeInPlace();
528// //
529// // May need a transpose if the number of channels is one
530// //
531// if (all.cols() == 1 && info.nchan == 1)
532// all.transposeInPlace();
533// }
534// else
535// {
536// //
537// // Put the old style epochs together
538// //
539// all = epoch[0]->toFloatMatrix();
540// all.transposeInPlace();
541
542// for (k = 2; k < nepoch; ++k)
543// {
544// oldsize = all.rows();
545// MatrixXd tmp = epoch[k]->toFloatMatrix();
546// tmp.transposeInPlace();
547// all.conservativeResize(oldsize+tmp.rows(), all.cols());
548// all.block(oldsize, 0, tmp.rows(), tmp.cols()) = tmp;
549// }
550// }
551// if (all.cols() != nsamp)
552// {
553// if(t_pStream)
554// delete t_pStream;
555// printf("Incorrect number of samples (%d instead of %d)", (int)all.cols(), nsamp);
556// return false;
557// }
558// //
559// // Calibrate
560// //
561// typedef Eigen::Triplet<double> T;
562// std::vector<T> tripletList;
563// tripletList.reserve(info.nchan);
564// for(k = 0; k < info.nchan; ++k)
565// tripletList.push_back(T(k, k, info.chs[k].cal));
566
567// SparseMatrix<double> cals(info.nchan, info.nchan);
568// cals.setFromTriplets(tripletList.begin(), tripletList.end());
569
570// all = cals * all;
571// //
572// // Put it all together
573// //
574// data.info = info;
575
578// data.evoked.append(FiffEvokedData::SDPtr(new FiffEvokedData()));
579// data.evoked[0]->aspect_kind = aspect_kind;
580// data.evoked[0]->is_smsh = is_smsh(0,setno);
581// if (nave != -1)
582// data.evoked[0]->nave = nave;
583// else
584// data.evoked[0]->nave = 1;
585
586// data.evoked[0]->first = first;
587// data.evoked[0]->last = last;
588// if (!comment.isEmpty())
589// data.evoked[0]->comment = comment;
590// //
591// // Times for convenience and the actual epoch data
592// //
593
594// data.evoked[0]->times = MatrixXd(1, last-first+1);
595
596// for (k = 0; k < data.evoked[0]->times.cols(); ++k)
597// data.evoked[0]->times(0, k) = ((float)(first+k)) / info.sfreq;
598
601// data.evoked[0]->epochs = all;
602
603// if(t_pStream)
604// delete t_pStream;
605
606// return true;
607}
FiffTag class declaration, which provides fiff tag I/O and processing methods.
FiffStream class declaration.
FiffDirNode class declaration, which provides fiff dir tree processing methods.
FiffEvokedSet class declaration.
CTF software compensation data.
QSharedPointer< FiffDirNode > SPtr
static bool read(QIODevice &p_IODevice, FiffEvoked &p_FiffEvoked, QVariant setno=0, QPair< float, float > t_baseline=defaultFloatPair, bool proj=true, fiff_int_t p_aspect_kind=FIFFV_ASPECT_AVERAGE)
bool compensate_to(FiffEvokedSet &p_FiffEvokedSet, fiff_int_t to) const
bool find_evoked(const FiffEvokedSet &p_FiffEvokedSet) const
FiffEvokedSet pick_channels(const QStringList &include=defaultQStringList, const QStringList &exclude=defaultQStringList) const
static bool read(QIODevice &p_IODevice, FiffEvokedSet &p_FiffEvokedSet, QPair< float, float > baseline=defaultFloatPair, bool proj=true)
QList< FiffEvoked > evoked
void set_current_comp(fiff_int_t value)
Definition fiff_info.h:292
qint32 get_current_comp()
bool make_compensator(fiff_int_t from, fiff_int_t to, FiffCtfComp &ctf_comp, bool exclude_comp_chs=false) const
FIFF File I/O routines.
QSharedPointer< FiffStream > SPtr