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