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