MNE-CPP  0.1.9
A Framework for Electrophysiology
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 
62 using namespace FIFFLIB;
63 using 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 
77 FiffEvokedSet::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 
100 {
101 }
102 
103 //=============================================================================================================
104 
106 {
107  info.clear();
108  evoked.clear();
109 }
110 
111 //=============================================================================================================
112 
113 FiffEvokedSet FiffEvokedSet::pick_channels(const QStringList& include,
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 
193 bool 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 %d 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 
211 bool 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  //
230  FiffDirNode::SPtr meas;
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 %d 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 }
static bool read(QIODevice &p_IODevice, FiffEvokedSet &p_FiffEvokedSet, QPair< float, float > baseline=defaultFloatPair, bool proj=true)
FiffEvokedSet pick_channels(const QStringList &include=defaultQStringList, const QStringList &exclude=defaultQStringList) const
FiffStream class declaration.
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)
QSharedPointer< FiffDirNode > SPtr
Definition: fiff_dir_node.h:77
evoked data set
QList< FiffEvoked > evoked
CTF software compensation data.
Definition: fiff_ctf_comp.h:73
FiffEvokedSet class declaration.
evoked data
Definition: fiff_evoked.h:77
FiffNamedMatrix::SDPtr data
FIFF File I/O routines.
Definition: fiff_stream.h:104
FiffDirNode class declaration, which provides fiff dir tree processing methods.
qint32 get_current_comp()
Definition: fiff_info.cpp:137
QSharedPointer< FiffStream > SPtr
Definition: fiff_stream.h:107
FiffTag class declaration, which provides fiff tag I/O and processing methods.
bool compensate_to(FiffEvokedSet &p_FiffEvokedSet, fiff_int_t to) const
void set_current_comp(fiff_int_t value)
Definition: fiff_info.h:292
bool find_evoked(const FiffEvokedSet &p_FiffEvokedSet) const
bool make_compensator(fiff_int_t from, fiff_int_t to, FiffCtfComp &ctf_comp, bool exclude_comp_chs=false) const
Definition: fiff_info.cpp:159