MNE-CPP  0.1.9
A Framework for Electrophysiology
fiff_evoked.cpp
Go to the documentation of this file.
1 //=============================================================================================================
37 //=============================================================================================================
38 // INCLUDES
39 //=============================================================================================================
40 
41 #include "fiff_evoked.h"
42 #include "fiff_stream.h"
43 #include "fiff_tag.h"
44 #include "fiff_dir_node.h"
45 
46 #include <utils/mnemath.h>
47 
48 //=============================================================================================================
49 // USED NAMESPACES
50 //=============================================================================================================
51 
52 using namespace FIFFLIB;
53 using namespace UTILSLIB;
54 using namespace Eigen;
55 
56 //=============================================================================================================
57 // DEFINE MEMBER METHODS
58 //=============================================================================================================
59 
61 : nave(-1)
62 , aspect_kind(-1)
63 , first(-1)
64 , last(-1)
65 , baseline(qMakePair(-1.0f, -1.0f))
66 {
67 
68 }
69 
70 //=============================================================================================================
71 
72 FiffEvoked::FiffEvoked(QIODevice& p_IODevice,
73  QVariant setno,
74  QPair<float,float> t_baseline,
75  bool proj,
76  fiff_int_t p_aspect_kind)
77 {
78  if(!FiffEvoked::read(p_IODevice, *this, setno, t_baseline, proj, p_aspect_kind))
79  {
80  baseline = t_baseline;
81 
82  printf("\tFiff evoked data not found.\n");//ToDo Throw here
83  return;
84  }
85 }
86 
87 //=============================================================================================================
88 
89 FiffEvoked::FiffEvoked(const FiffEvoked& p_FiffEvoked)
90 : info(p_FiffEvoked.info)
91 , nave(p_FiffEvoked.nave)
92 , aspect_kind(p_FiffEvoked.aspect_kind)
93 , first(p_FiffEvoked.first)
94 , last(p_FiffEvoked.last)
95 , comment(p_FiffEvoked.comment)
96 , times(p_FiffEvoked.times)
97 , data(p_FiffEvoked.data)
98 , proj(p_FiffEvoked.proj)
99 , baseline(p_FiffEvoked.baseline)
100 {
101 
102 }
103 
104 //=============================================================================================================
105 
107 {
108 
109 }
110 
111 //=============================================================================================================
112 
114 {
115  info.clear();
116  nave = -1;
117  aspect_kind = -1;
118  first = -1;
119  last = -1;
120  comment = QString("");
121  times = RowVectorXf();
122  data = MatrixXd();
123  proj = MatrixXd();
124 }
125 
126 //=============================================================================================================
127 
128 FiffEvoked FiffEvoked::pick_channels(const QStringList& include,
129  const QStringList& exclude) const
130 {
131  if(include.size() == 0 && exclude.size() == 0)
132  return FiffEvoked(*this);
133 
134  RowVectorXi sel = FiffInfo::pick_channels(this->info.ch_names, include, exclude);
135  if (sel.cols() == 0)
136  {
137  qWarning("Warning : No channels match the selection.\n");
138  return FiffEvoked(*this);
139  }
140 
141  FiffEvoked res(*this);
142  //
143  // Modify the measurement info
144  //
145  res.info = FiffInfo(res.info.pick_info(sel));
146  //
147  // Create the reduced data set
148  //
149  MatrixXd selBlock(1,1);
150 
151  if(selBlock.rows() != sel.cols() || selBlock.cols() != res.data.cols())
152  selBlock.resize(sel.cols(), res.data.cols());
153  for(qint32 l = 0; l < sel.cols(); ++l)
154  {
155  if(sel(0,l) <= res.data.rows()) {
156  selBlock.block(l,0,1,selBlock.cols()) = res.data.block(sel(0,l),0,1,selBlock.cols());
157  } else {
158  qWarning("FiffEvoked::pick_channels - Warning : Selected channel index out of bound.\n");
159  }
160  }
161  res.data.resize(sel.cols(), res.data.cols());
162  res.data = selBlock;
163 
164  return res;
165 }
166 
167 //=============================================================================================================
168 
169 bool FiffEvoked::read(QIODevice& p_IODevice,
170  FiffEvoked& p_FiffEvoked,
171  QVariant setno,
172  QPair<float,float> t_baseline,
173  bool proj,
174  fiff_int_t p_aspect_kind)
175 {
176  p_FiffEvoked.clear();
177 
178  //
179  // Open the file
180  //
181  FiffStream::SPtr t_pStream(new FiffStream(&p_IODevice));
182  QString t_sFileName = t_pStream->streamName();
183 
184  printf("Reading %s ...\n",t_sFileName.toUtf8().constData());
185 
186  if(!t_pStream->open())
187  return false;
188  //
189  // Read the measurement info
190  //
191  FiffInfo info;
192  FiffDirNode::SPtr meas;
193  if(!t_pStream->read_meas_info(t_pStream->dirtree(), info, meas))
194  return false;
195  info.filename = t_sFileName; //move fname storage to read_meas_info member function
196  //
197  // Locate the data of interest
198  //
199  QList<FiffDirNode::SPtr> processed = meas->dir_tree_find(FIFFB_PROCESSED_DATA);
200  if (processed.size() == 0)
201  {
202  qWarning("Could not find processed data");
203  return false;
204  }
205  //
206  QList<FiffDirNode::SPtr> evoked_node = meas->dir_tree_find(FIFFB_EVOKED);
207  if (evoked_node.size() == 0)
208  {
209  qWarning("Could not find evoked data");
210  return false;
211  }
212 
213  // convert setno to an integer
214  if(!setno.isValid())
215  {
216  if (evoked_node.size() > 1)
217  {
218  QStringList comments;
219  QList<fiff_int_t> aspect_kinds;
220  QString t;
221  if(!t_pStream->get_evoked_entries(evoked_node, comments, aspect_kinds, t))
222  t = QString("None found, must use integer");
223  qWarning("%d datasets present, setno parameter must be set. Candidate setno names:\n%s", evoked_node.size(), t.toUtf8().constData());
224  return false;
225  }
226  else
227  setno = 0;
228  }
229  else
230  {
231  // find string-based entry
232  bool t_bIsInteger = true;
233  setno.toInt(&t_bIsInteger);
234  if(!t_bIsInteger)
235  {
236  if(p_aspect_kind != FIFFV_ASPECT_AVERAGE && p_aspect_kind != FIFFV_ASPECT_STD_ERR)
237  {
238  qWarning("kindStat must be \"FIFFV_ASPECT_AVERAGE\" or \"FIFFV_ASPECT_STD_ERR\"");
239  return false;
240  }
241 
242  QStringList comments;
243  QList<fiff_int_t> aspect_kinds;
244  QString t;
245  t_pStream->get_evoked_entries(evoked_node, comments, aspect_kinds, t);
246 
247  bool found = false;
248  for(qint32 i = 0; i < comments.size(); ++i)
249  {
250  if(comments[i].compare(setno.toString()) == 0 && p_aspect_kind == aspect_kinds[i])
251  {
252  setno = i;
253  found = true;
254  break;
255  }
256  }
257  if(!found)
258  {
259  qWarning() << "setno " << setno << " (" << p_aspect_kind << ") not found, out of found datasets:\n " << t;
260  return false;
261  }
262  }
263  }
264 
265  if (setno.toInt() >= evoked_node.size() || setno.toInt() < 0)
266  {
267  qWarning("Data set selector out of range");
268  return false;
269  }
270 
271  FiffDirNode::SPtr my_evoked = evoked_node[setno.toInt()];
272 
273  //
274  // Identify the aspects
275  //
276  QList<FiffDirNode::SPtr> aspects = my_evoked->dir_tree_find(FIFFB_ASPECT);
277 
278  if(aspects.size() > 1)
279  printf("\tMultiple (%d) aspects found. Taking first one.\n", aspects.size());
280 
281  FiffDirNode::SPtr my_aspect = aspects[0];
282 
283  //
284  // Now find the data in the evoked block
285  //
286  fiff_int_t nchan = 0;
287  float sfreq = -1.0f;
288  QList<FiffChInfo> chs;
289  fiff_int_t kind, pos, first=0, last=0;
290  FiffTag::SPtr t_pTag;
291  QString comment("");
292  qint32 k;
293  for (k = 0; k < my_evoked->nent(); ++k)
294  {
295  kind = my_evoked->dir[k]->kind;
296  pos = my_evoked->dir[k]->pos;
297  switch (kind)
298  {
299  case FIFF_COMMENT:
300  t_pStream->read_tag(t_pTag,pos);
301  comment = t_pTag->toString();
302  break;
303  case FIFF_FIRST_SAMPLE:
304  t_pStream->read_tag(t_pTag,pos);
305  first = *t_pTag->toInt();
306  break;
307  case FIFF_LAST_SAMPLE:
308  t_pStream->read_tag(t_pTag,pos);
309  last = *t_pTag->toInt();
310  break;
311  case FIFF_NCHAN:
312  t_pStream->read_tag(t_pTag,pos);
313  nchan = *t_pTag->toInt();
314  break;
315  case FIFF_SFREQ:
316  t_pStream->read_tag(t_pTag,pos);
317  sfreq = *t_pTag->toFloat();
318  break;
319  case FIFF_CH_INFO:
320  t_pStream->read_tag(t_pTag, pos);
321  chs.append( t_pTag->toChInfo() );
322  break;
323  }
324  }
325  if (comment.isEmpty())
326  comment = QString("No comment");
327 
328  //
329  // Local channel information?
330  //
331  if (nchan > 0)
332  {
333  if (chs.size() == 0)
334  {
335  qWarning("Local channel information was not found when it was expected.");
336  return false;
337  }
338  if (chs.size() != nchan)
339  {
340  qWarning("Number of channels and number of channel definitions are different.");
341  return false;
342  }
343  info.chs = chs;
344  info.nchan = nchan;
345  printf("\tFound channel information in evoked data. nchan = %d\n",nchan);
346  if (sfreq > 0.0f)
347  info.sfreq = sfreq;
348  }
349  qint32 nsamp = last-first+1;
350  printf("\tFound the data of interest:\n");
351  printf("\t\tt = %10.2f ... %10.2f ms (%s)\n", 1000*(float)first/info.sfreq, 1000*(float)last/info.sfreq,comment.toUtf8().constData());
352  if (info.comps.size() > 0)
353  printf("\t\t%d CTF compensation matrices available\n", info.comps.size());
354 
355  //
356  // Read the data in the aspect block
357  //
358  fiff_int_t aspect_kind = -1;
359  fiff_int_t nave = -1;
360  QList<FiffTag> epoch;
361  for (k = 0; k < my_aspect->nent(); ++k)
362  {
363  kind = my_aspect->dir[k]->kind;
364  pos = my_aspect->dir[k]->pos;
365 
366  switch (kind)
367  {
368  case FIFF_COMMENT:
369  t_pStream->read_tag(t_pTag, pos);
370  comment = t_pTag->toString();
371  break;
372  case FIFF_ASPECT_KIND:
373  t_pStream->read_tag(t_pTag, pos);
374  aspect_kind = *t_pTag->toInt();
375  break;
376  case FIFF_NAVE:
377  t_pStream->read_tag(t_pTag, pos);
378  nave = *t_pTag->toInt();
379  break;
380  case FIFF_EPOCH:
381  t_pStream->read_tag(t_pTag, pos);
382  epoch.append(FiffTag(t_pTag.data()));
383  break;
384  }
385  }
386  if (nave == -1)
387  nave = 1;
388  printf("\t\tnave = %d - aspect type = %d\n", nave, aspect_kind);
389 
390  qint32 nepoch = epoch.size();
391  MatrixXd all_data;
392  if (nepoch == 1)
393  {
394  //
395  // Only one epoch
396  //
397  all_data = epoch[0].toFloatMatrix().cast<double>();
398  all_data.transposeInPlace();
399  //
400  // May need a transpose if the number of channels is one
401  //
402  if (all_data.cols() == 1 && info.nchan == 1)
403  all_data.transposeInPlace();
404  }
405  else
406  {
407  //
408  // Put the old style epochs together
409  //
410  all_data = epoch[0].toFloatMatrix().cast<double>();
411  all_data.transposeInPlace();
412  qint32 oldsize;
413  for (k = 1; k < nepoch; ++k)
414  {
415  oldsize = all_data.rows();
416  MatrixXd tmp = epoch[k].toFloatMatrix().cast<double>();
417  tmp.transposeInPlace();
418  all_data.conservativeResize(oldsize+tmp.rows(), all_data.cols());
419  all_data.block(oldsize, 0, tmp.rows(), tmp.cols()) = tmp;
420  }
421  }
422  if (all_data.cols() != nsamp)
423  {
424  qWarning("Incorrect number of samples (%d instead of %d)", (int) all_data.cols(), nsamp);
425  return false;
426  }
427 
428  //
429  // Calibrate
430  //
431  printf("\n\tPreprocessing...\n");
432  printf("\t%d channels remain after picking\n",info.nchan);
433 
434  typedef Eigen::Triplet<double> T;
435  std::vector<T> tripletList;
436  tripletList.reserve(info.nchan);
437  for(k = 0; k < info.nchan; ++k)
438  tripletList.push_back(T(k, k, info.chs[k].cal));
439  SparseMatrix<double> cals(info.nchan, info.nchan);
440  cals.setFromTriplets(tripletList.begin(), tripletList.end());
441 
442  all_data = cals * all_data;
443 
444  RowVectorXf times = RowVectorXf(last-first+1);
445  for (k = 0; k < times.size(); ++k)
446  times[k] = ((float)(first+k)) / info.sfreq;
447 
448  //
449  // Set up projection
450  //
451  if(info.projs.size() == 0 || !proj)
452  {
453  printf("\tNo projector specified for these data.\n");
454  p_FiffEvoked.proj = MatrixXd();
455  }
456  else
457  {
458  // Create the projector
459  MatrixXd projection;
460  qint32 nproj = info.make_projector(projection);
461  if(nproj == 0)
462  {
463  printf("\tThe projection vectors do not apply to these channels\n");
464  p_FiffEvoked.proj = MatrixXd();
465  }
466  else
467  {
468  printf("\tCreated an SSP operator (subspace dimension = %d)\n", nproj);
469  p_FiffEvoked.proj = projection;
470  }
471 
472  // The projection items have been activated
474  }
475 
476  if(p_FiffEvoked.proj.rows() > 0)
477  {
478  all_data = p_FiffEvoked.proj * all_data;
479  printf("\tSSP projectors applied to the evoked data\n");
480  }
481 
482  // Put it all together
483  p_FiffEvoked.info = info;
484  p_FiffEvoked.nave = nave;
485  p_FiffEvoked.aspect_kind = aspect_kind;
486  p_FiffEvoked.first = first;
487  p_FiffEvoked.last = last;
488  p_FiffEvoked.comment = comment;
489  p_FiffEvoked.times = times;
490  p_FiffEvoked.data = all_data;
491 
492  // Run baseline correction
493  p_FiffEvoked.applyBaselineCorrection(t_baseline);
494 
495  return true;
496 }
497 
498 //=============================================================================================================
499 
500 void FiffEvoked::setInfo(const FiffInfo &p_info,
501  bool proj)
502 {
503  info = p_info;
504  //
505  // Set up projection
506  //
507  if(info.projs.size() == 0 || !proj)
508  {
509  printf("\tNo projector specified for these data.\n");
510  this->proj = MatrixXd();
511  }
512  else
513  {
514  // Create the projector
515  MatrixXd projection;
516  qint32 nproj = info.make_projector(projection);
517  if(nproj == 0)
518  {
519  printf("\tThe projection vectors do not apply to these channels\n");
520  this->proj = MatrixXd();
521  }
522  else
523  {
524  printf("\tCreated an SSP operator (subspace dimension = %d)\n", nproj);
525  this->proj = projection;
526  }
527 
528  // The projection items have been activated
530  }
531 }
532 
533 //=============================================================================================================
534 
535 FiffEvoked & FiffEvoked::operator+=(const MatrixXd &newData)
536 {
537  //Init matrix if necessary
538  if(nave == -1 || nave == 0) {
539  data = MatrixXd::Zero(newData.rows(),newData.cols());
540  }
541 
542  if(data.cols() == newData.cols() && data.rows() == newData.rows()) {
543  //Revert old averaging
544  data = data*nave;
545 
546  //Do new averaging
547  data += newData;
548  if(nave <= 0) {
549  nave = 1;
550  } else {
551  nave++;
552  }
553 
554  data /= nave;
555  }
556 
557  return *this;
558 }
559 
560 //=============================================================================================================
561 
562 void FiffEvoked::applyBaselineCorrection(QPair<float, float>& p_baseline)
563 {
564  // Run baseline correction
565  printf("Applying baseline correction ... (mode: mean)\n");
566  this->data = MNEMath::rescale(this->data, this->times, p_baseline, QString("mean"));
567  this->baseline = p_baseline;
568 }
#define FIFF_FIRST_SAMPLE
Definition: fiff_file.h:461
static Eigen::RowVectorXi pick_channels(const QStringList &ch_names, const QStringList &include=defaultQStringList, const QStringList &exclude=defaultQStringList)
#define FIFF_EPOCH
Definition: fiff_file.h:560
FiffEvoked & operator+=(const Eigen::MatrixXd &newData)
QList< FiffCtfComp > comps
Definition: fiff_info.h:269
Eigen::RowVectorXf times
Definition: fiff_evoked.h:230
static void activate_projs(QList< FiffProj > &p_qListFiffProj)
Definition: fiff_proj.cpp:104
FiffStream class declaration.
#define FIFF_COMMENT
Definition: fiff_file.h:459
FIFF measurement file information.
Definition: fiff_info.h:84
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)
#define FIFF_NAVE
Definition: fiff_file.h:460
QSharedPointer< FiffDirNode > SPtr
Definition: fiff_dir_node.h:77
#define FIFFV_ASPECT_STD_ERR
Definition: fiff_file.h:439
FiffEvoked class declaration.
QPair< float, float > baseline
Definition: fiff_evoked.h:233
FIFF data tag.
Definition: fiff_tag.h:148
#define FIFF_CH_INFO
Definition: fiff_file.h:456
#define FIFF_NCHAN
Definition: fiff_file.h:453
fiff_int_t aspect_kind
Definition: fiff_evoked.h:226
FiffInfo pick_info(const Eigen::RowVectorXi &sel=defaultVectorXi) const
Definition: fiff_info.cpp:310
QList< FiffProj > projs
Definition: fiff_info.h:268
QSharedPointer< FiffTag > SPtr
Definition: fiff_tag.h:152
evoked data
Definition: fiff_evoked.h:77
FIFF File I/O routines.
Definition: fiff_stream.h:104
#define FIFF_LAST_SAMPLE
Definition: fiff_file.h:462
#define FIFF_SFREQ
Definition: fiff_file.h:454
FiffDirNode class declaration, which provides fiff dir tree processing methods.
qint32 make_projector(Eigen::MatrixXd &proj) const
Definition: fiff_info.h:278
MNEMath class declaration.
#define FIFF_ASPECT_KIND
Definition: fiff_file.h:463
#define FIFFV_ASPECT_AVERAGE
Definition: fiff_file.h:438
QSharedPointer< FiffStream > SPtr
Definition: fiff_stream.h:107
FiffTag class declaration, which provides fiff tag I/O and processing methods.
FiffEvoked pick_channels(const QStringList &include=defaultQStringList, const QStringList &exclude=defaultQStringList) const
void setInfo(const FiffInfo &p_info, bool proj=true)
Eigen::MatrixXd proj
Definition: fiff_evoked.h:232
Eigen::MatrixXd data
Definition: fiff_evoked.h:231
void applyBaselineCorrection(QPair< float, float > &p_baseline)
QList< FiffChInfo > chs