MNE-CPP  0.1.9
A Framework for Electrophysiology
mne_meas_data.cpp
Go to the documentation of this file.
1 //=============================================================================================================
37 //=============================================================================================================
38 // INCLUDES
39 //=============================================================================================================
40 
41 #include "mne_meas_data.h"
42 #include "mne_meas_data_set.h"
43 #include "mne_inverse_operator.h"
44 
45 #include <mne/c/mne_types.h>
46 #include <mne/c/mne_named_matrix.h>
47 
48 #include <fiff/fiff_types.h>
49 
50 #include <time.h>
51 
52 #include <QFile>
53 
54 //=============================================================================================================
55 // USED NAMESPACES
56 //=============================================================================================================
57 
58 using namespace Eigen;
59 using namespace FIFFLIB;
60 using namespace MNELIB;
61 using namespace INVERSELIB;
62 
63 #ifndef TRUE
64 #define TRUE 1
65 #endif
66 
67 #ifndef FALSE
68 #define FALSE 0
69 #endif
70 
71 #ifndef FAIL
72 #define FAIL -1
73 #endif
74 
75 #ifndef OK
76 #define OK 0
77 #endif
78 
79 #if defined(_WIN32) || defined(_WIN64)
80 #define snprintf _snprintf
81 #define vsnprintf _vsnprintf
82 #define strcasecmp _stricmp
83 #define strncasecmp _strnicmp
84 #endif
85 
86 #define MIN(a,b) ((a) < (b) ? (a) : (b))
87 #define MAX(a,b) ((a) > (b) ? (a) : (b))
88 
89 #define MALLOC_9(x,t) (t *)malloc((x)*sizeof(t))
90 #define REALLOC_9(x,y,t) (t *)((x == NULL) ? malloc((y)*sizeof(t)) : realloc((x),(y)*sizeof(t)))
91 
92 #define FREE_9(x) if ((char *)(x) != NULL) free((char *)(x))
93 
94 #define ALLOC_CMATRIX_9(x,y) mne_cmatrix_9((x),(y))
95 
96 #define FREE_CMATRIX_9(m) mne_free_cmatrix_9((m))
97 
98 void fromFloatEigenMatrix_9(const Eigen::MatrixXf& from_mat, float **& to_mat, const int m, const int n)
99 {
100  for ( int i = 0; i < m; ++i)
101  for ( int j = 0; j < n; ++j)
102  to_mat[i][j] = from_mat(i,j);
103 }
104 
105 void fromFloatEigenMatrix_9(const Eigen::MatrixXf& from_mat, float **& to_mat)
106 {
107  fromFloatEigenMatrix_9(from_mat, to_mat, from_mat.rows(), from_mat.cols());
108 }
109 
110 void mne_free_cmatrix_9 (float **m)
111 {
112  if (m) {
113  FREE_9(*m);
114  FREE_9(m);
115  }
116 }
117 
118 static void matrix_error_9(int kind, int nr, int nc)
119 
120 {
121  if (kind == 1)
122  printf("Failed to allocate memory pointers for a %d x %d matrix\n",nr,nc);
123  else if (kind == 2)
124  printf("Failed to allocate memory for a %d x %d matrix\n",nr,nc);
125  else
126  printf("Allocation error for a %d x %d matrix\n",nr,nc);
127  if (sizeof(void *) == 4) {
128  printf("This is probably because you seem to be using a computer with 32-bit architecture.\n");
129  printf("Please consider moving to a 64-bit platform.");
130  }
131  printf("Cannot continue. Sorry.\n");
132  exit(1);
133 }
134 
135 float **mne_cmatrix_9(int nr,int nc)
136 
137 {
138  int i;
139  float **m;
140  float *whole;
141 
142  m = MALLOC_9(nr,float *);
143  if (!m) matrix_error_9(1,nr,nc);
144  whole = MALLOC_9(nr*nc,float);
145  if (!whole) matrix_error_9(2,nr,nc);
146 
147  for(i=0;i<nr;i++)
148  m[i] = whole + i*nc;
149  return m;
150 }
151 
152 namespace INVERSELIB
153 {
154 
155 typedef struct {
156  int size; /* Size of this buffer in floats */
157  float *data; /* The allocated buffer */
158  float ***datap; /* The matrix which uses this */
159 } *ringBufBuf_9,ringBufBufRec_9;
160 
161 typedef struct {
162  ringBufBuf_9 *bufs;
163  int nbuf;
164  int next;
165 } *ringBuf_9,ringBufRec_9;
166 
167 }
168 
169 void mne_free_ring_buffer_9(void *thisp)
170 
171 {
172  int k;
173  ringBuf_9 this_buf = (ringBuf_9)thisp;
174 
175  if (!this_buf)
176  return;
177 
178  for (k = 0; k < this_buf->nbuf; k++)
179  FREE_9(this_buf->bufs[k]->data);
180  FREE_9(this_buf->bufs);
181  FREE_9(this_buf);
182  return;
183 }
184 
185 //============================= mne_events.c =============================
186 
187 void mne_free_event_9(mneEvent e)
188 {
189  if (!e)
190  return;
191  FREE_9(e->comment);
192  FREE_9(e);
193  return;
194 }
195 
196 void mne_free_event_list_9(mneEventList list)
197 
198 {
199  int k;
200  if (!list)
201  return;
202  for (k = 0; k < list->nevent; k++)
203  mne_free_event_9(list->events[k]);
204  FREE_9(list->events);
205  FREE_9(list);
206  return;
207 }
208 
209 void mne_free_sparse_named_matrix_9(mneSparseNamedMatrix mat)
210 /*
211  * Free the matrix and all the data from within
212  */
213 {
214  if (!mat)
215  return;
216  mat->rowlist.clear();
217  mat->collist.clear();
218  if(mat->data)
219  delete mat->data;
220  FREE_9(mat);
221  return;
222 }
223 
224 //============================= mne_raw_routines.c =============================
225 
226 void mne_ch_selection_free_9(mneChSelection s)
227 
228 {
229  if (!s)
230  return;
231  s->name.clear();
232  FREE_9(s->pick);
233  FREE_9(s->pick_deriv);
234  FREE_9(s->ch_kind);
235  s->chspick.clear();
236  s->chspick_nospace.clear();
237  s->chdef.clear();
238  FREE_9(s);
239  return;
240 }
241 
242 void mne_string_to_name_list_9(const QString& s, QStringList& listp,int &nlistp)
243 /*
244  * Convert a colon-separated list into a string array
245  */
246 {
247  QStringList list;
248 
249  if (!s.isEmpty() && s.size() > 0) {
250  list = FIFFLIB::FiffStream::split_name_list(s);
251  //list = s.split(":");
252  }
253  listp = list;
254  nlistp = list.size();
255  return;
256 }
257 
258 QString mne_name_list_to_string_9(const QStringList& list)
259 /*
260  * Convert a string array to a colon-separated string
261  */
262 {
263  int nlist = list.size();
264  QString res;
265  if (nlist == 0 || list.isEmpty())
266  return res;
267  // res[0] = '\0';
268  for (int k = 0; k < nlist-1; k++) {
269  res += list[k];
270  res += ":";
271  }
272  res += list[nlist-1];
273  return res;
274 }
275 
276 QString mne_channel_names_to_string_9(const QList<FIFFLIB::FiffChInfo>& chs,
277  int nch)
278 /*
279  * Make a colon-separated string out of channel names
280  */
281 {
282  QStringList names;
283  QString res;
284  if (nch <= 0)
285  return res;
286  for (int k = 0; k < nch; k++)
287  names.append(chs[k].ch_name);
288  res = mne_name_list_to_string_9(names);
289  return res;
290 }
291 
292 void mne_channel_names_to_name_list_9(const QList<FIFFLIB::FiffChInfo>& chs,
293  int nch,
294  QStringList& listp,
295  int &nlistp)
296 
297 {
298  QString s = mne_channel_names_to_string_9(chs,nch);
299 
300  mne_string_to_name_list_9(s,listp,nlistp);
301  return;
302 }
303 
304 //============================= read_ch_info.c =============================
305 
306 static FiffDirNode::SPtr find_meas_9 (const FiffDirNode::SPtr& node)
307 /*
308  * Find corresponding meas node
309  */
310 {
311  FiffDirNode::SPtr empty_node;
312  FiffDirNode::SPtr tmp_node = node;
313 
314  while (tmp_node->type != FIFFB_MEAS) {
315  if (tmp_node->parent == NULL)
316  return empty_node;//(NULL);
317  tmp_node = tmp_node->parent;
318  }
319  return (tmp_node);
320 }
321 
322 static FiffDirNode::SPtr find_meas_info_9 (const FiffDirNode::SPtr& node)
323 /*
324  * Find corresponding meas info node
325  */
326 {
327  int k;
328  FiffDirNode::SPtr empty_node;
329  FiffDirNode::SPtr tmp_node = node;
330 
331  while (tmp_node->type != FIFFB_MEAS) {
332  if (tmp_node->parent == NULL)
333  return empty_node;
334  tmp_node = tmp_node->parent;
335  }
336  for (k = 0; k < tmp_node->nchild(); k++)
337  if (tmp_node->children[k]->type == FIFFB_MEAS_INFO)
338  return (tmp_node->children[k]);
339  return empty_node;
340 }
341 
342 //============================= mne_read_evoked.c =============================
343 
344 #define MAXDATE 100
345 
346 static FiffDirNode::SPtr find_evoked (const FiffDirNode::SPtr& node)
347 /*
348  * Find corresponding FIFFB_EVOKED node
349  */
350 {
351  FiffDirNode::SPtr empty_node;
352  FiffDirNode::SPtr tmp_node = node;
353  while (tmp_node->type != FIFFB_EVOKED) {
354  if (tmp_node->parent == NULL)
355  return empty_node;
356  tmp_node = tmp_node->parent;
357  }
358  return (tmp_node);
359 }
360 
361 static QString get_comment ( FiffStream::SPtr& stream,
362  const FiffDirNode::SPtr& start)
363 
364 {
365  int k;
366  FiffTag::SPtr t_pTag;
367  QList<FiffDirEntry::SPtr> ent = start->dir;
368  for (k = 0; k < start->nent(); k++)
369  if (ent[k]->kind == FIFF_COMMENT) {
370  if (stream->read_tag(t_pTag,ent[k]->pos)) {
371  return t_pTag->toString();
372  }
373  }
374  return QString("No comment");
375 }
376 
377 static void get_aspect_name_type( FiffStream::SPtr& stream,
378  const FiffDirNode::SPtr& start,
379  QString& namep, int *typep)
380 
381 {
382  int k;
383  FiffTag::SPtr t_pTag;
384  QList<FiffDirEntry::SPtr> ent = start->dir;
385  QString res = "unknown";
386  int type = -1;
387 
388  for (k = 0; k < start->nent(); k++)
389  if (ent[k]->kind == FIFF_ASPECT_KIND) {
390  if (stream->read_tag(t_pTag,ent[k]->pos)) {
391  type = *t_pTag->toInt();
392  switch (type) {
393  case FIFFV_ASPECT_AVERAGE :
394  res = "average";
395  break;
396  case FIFFV_ASPECT_STD_ERR :
397  res = "std.error";
398  break;
399  case FIFFV_ASPECT_SINGLE :
400  res = "single trace";
401  break;
402  case FIFFV_ASPECT_SAMPLE :
403  res = "sample";
404  break;
405  case FIFFV_ASPECT_SUBAVERAGE :
406  res = "subaverage";
407  break;
408  case FIFFV_ASPECT_ALTAVERAGE :
409  res = "alt. average";
410  break;
411  case FIFFV_ASPECT_POWER_DENSITY :
412  res = "power density spectrum";
413  break;
414  case FIFFV_ASPECT_DIPOLE_WAVE :
415  res = "dipole amplitudes";
416  break;
417  }
418  }
419  break;
420  }
421  // if (namep)
422  namep = res;
423  if (typep)
424  *typep = type;
425  return;
426 }
427 
428 static char *get_meas_date ( FiffStream::SPtr& stream,const FiffDirNode::SPtr& node )
429 {
430  int k;
431  FiffTag::SPtr t_pTag;
432  char *res = NULL;
433  fiff_int_t kind, pos;
434  FiffDirNode::SPtr meas_info;
435 
436  if (!(meas_info = find_meas_info_9(node))) {
437  return res;
438  }
439  for (k = 0; k < meas_info->nent();k++) {
440  kind = meas_info->dir[k]->kind;
441  pos = meas_info->dir[k]->pos;
442  if (kind == FIFF_MEAS_DATE)
443  {
444  if (stream->read_tag(t_pTag,pos)) {
445  fiffTime meas_date = (fiffTime)t_pTag->data();
446  time_t time = meas_date->secs;
447  struct tm *ltime;
448 
449  ltime = localtime(&time);
450  res = MALLOC_9(MAXDATE,char);
451  (void)strftime(res,MAXDATE,"%x %X",ltime);
452  break;
453  }
454  }
455  }
456  return res;
457 }
458 
459 static int get_meas_info ( FiffStream::SPtr& stream, /* The stream we are reading */
460  const FiffDirNode::SPtr& node, /* The directory node containing our data */
461  fiffId *id, /* The block id from the nearest FIFFB_MEAS parent */
462  fiffTime *meas_date, /* Measurement date */
463  int *nchan, /* Number of channels */
464  float *sfreq, /* Sampling frequency */
465  float *highpass, /* Highpass filter setting */
466  float *lowpass, /* Lowpass filter setting */
467  QList<FIFFLIB::FiffChInfo>& chp, /* Channel descriptions */
468  FiffCoordTransOld* *trans) /* Coordinate transformation (head <-> device) */
469 /*
470  * Find channel information from
471  * nearest FIFFB_MEAS_INFO parent of
472  * node.
473  */
474 {
475  QList<FIFFLIB::FiffChInfo> ch;
476  FIFFLIB::FiffChInfo this_ch;
477  FiffCoordTransOld* t = Q_NULLPTR;
478  int j,k;
479  int to_find = 4;
480  QList<FiffDirNode::SPtr> hpi;
481  FiffDirNode::SPtr meas;
482  FiffDirNode::SPtr meas_info;
483  fiff_int_t kind, pos;
484  FiffTag::SPtr t_pTag;
485 
486  *trans = NULL;
487  *id = NULL;
488  /*
489  * Find desired parents
490  */
491  if (!(meas = find_meas_9(node))) {
492  printf ("Meas. block not found!");
493  goto bad;
494  }
495  if (!(meas_info = find_meas_info_9(node))) {
496  printf ("Meas. info not found!");
497  goto bad;
498  }
499  /*
500  * Is there a block id is in the FIFFB_MEAS node?
501  */
502  if (!meas->id.isEmpty()) {
503  *id = MALLOC_9(1,fiffIdRec);
504  (*id)->version = meas->id.version;
505  (*id)->machid[0] = meas->id.machid[0];
506  (*id)->machid[1] = meas->id.machid[1];
507  (*id)->time = meas->id.time;
508  }
509  /*
510  * Others from FIFFB_MEAS_INFO
511  */
512  *lowpass = -1;
513  *highpass = -1;
514  for (k = 0; k < meas_info->nent(); k++) {
515  kind = meas_info->dir[k]->kind;
516  pos = meas_info->dir[k]->pos;
517  switch (kind) {
518 
519  case FIFF_NCHAN :
520  if (!stream->read_tag(t_pTag,pos))
521  goto bad;
522  *nchan = *t_pTag->toInt();
523 
524  for (j = 0; j < *nchan; j++) {
525  ch.append(FiffChInfo());
526  ch[j].scanNo = -1;
527  }
528  to_find = to_find + *nchan - 1;
529  break;
530 
531  case FIFF_SFREQ :
532  if (!stream->read_tag(t_pTag,pos))
533  goto bad;
534  *sfreq = *t_pTag->toFloat();
535  to_find--;
536  break;
537 
538  case FIFF_MEAS_DATE :
539  if (!stream->read_tag(t_pTag,pos))
540  goto bad;
541  if (*meas_date)
542  FREE_9(*meas_date);
543  *meas_date = MALLOC_9(1,fiffTimeRec);
544  **meas_date = *(fiffTime)t_pTag->data();
545  break;
546 
547  case FIFF_LOWPASS :
548  if (!stream->read_tag(t_pTag,pos))
549  goto bad;
550  *lowpass = *t_pTag->toFloat();
551  to_find--;
552  break;
553 
554  case FIFF_HIGHPASS :
555  if (!stream->read_tag(t_pTag,pos))
556  goto bad;
557  *highpass = *t_pTag->toFloat();
558  to_find--;
559  break;
560 
561  case FIFF_CH_INFO : /* Information about one channel */
562 
563  if (!stream->read_tag(t_pTag,pos))
564  goto bad;
565 
566  this_ch = t_pTag->toChInfo();
567  if (this_ch.scanNo <= 0 || this_ch.scanNo > *nchan) {
568  qCritical ("FIFF_CH_INFO : scan # out of range!");
569  goto bad;
570  }
571  else
572  ch[this_ch.scanNo-1] = this_ch;
573  to_find--;
574  break;
575 
576  case FIFF_COORD_TRANS :
577  if (!stream->read_tag(t_pTag,pos))
578  goto bad;
579  if(t)
580  delete t;
581  t = FiffCoordTransOld::read_helper( t_pTag );
582  /*
583  * Require this particular transform!
584  */
585  if (t->from == FIFFV_COORD_DEVICE && t->to == FIFFV_COORD_HEAD) {
586  *trans = t;
587  break;
588  }
589  }
590  }
591  /*
592  * Search for the coordinate transformation from
593  * HPI_RESULT block if it was not previously found
594  */
595 
596  hpi = meas_info->dir_tree_find(FIFFB_HPI_RESULT);
597 
598  if (hpi.size() > 0 && *trans == NULL)
599  for (k = 0; k < hpi[0]->nent(); k++)
600  if (hpi[0]->dir[k]->kind == FIFF_COORD_TRANS) {
601  if (!stream->read_tag(t_pTag,hpi[0]->dir[k]->pos))
602  goto bad;
603  t = FiffCoordTransOld::read_helper( t_pTag );
604 
605  /*
606  * Require this particular transform!
607  */
608  if (t->from == FIFFV_COORD_DEVICE && t->to == FIFFV_COORD_HEAD) {
609  *trans = t;
610  break;
611  }
612  }
613  if (*lowpass < 0) {
614  *lowpass = *sfreq/2.0;
615  to_find--;
616  }
617  if (*highpass < 0) {
618  *highpass = 0.0;
619  to_find--;
620  }
621  if (to_find != 0) {
622  printf ("Not all essential tags were found!");
623  goto bad;
624  }
625 
626  chp = ch;
627  return (0);
628 
629 bad : {
630  return (-1);
631  }
632 }
633 
634 static int find_between ( FiffStream::SPtr& stream,
635  const FiffDirNode::SPtr& low_node,
636  const FiffDirNode::SPtr& high_node,
637  int kind,
638  fiff_byte_t **data)
639 {
640  FiffTag::SPtr t_pTag;
641  FiffDirNode::SPtr node;
642  fiffDirEntry dir;
643  fiff_int_t kind_1, pos;
644  int k;
645 
646  *data = NULL;
647  node = low_node;
648  while (node != NULL) {
649  for (k = 0; k < node->nent(); k++)
650  {
651  kind_1 = node->dir[k]->kind;
652  pos = node->dir[k]->pos;
653  if (kind_1 == kind) {
654  FREE_9(*data);
655  if (!stream->read_tag(t_pTag,pos)) {
656  return (FIFF_FAIL);
657  }
658  else {
659  fiff_byte_t* tmp = MALLOC_9(t_pTag->size(),fiff_byte_t);
660  fiff_byte_t* tmp_current = (fiff_byte_t *)t_pTag->data();
661 
662  for( int k = 0; k < t_pTag->size(); ++k )
663  tmp[k] = tmp_current[k];
664 
665  *data = tmp;
666  return (FIFF_OK);
667  }
668  }
669  }
670  if (node == high_node)
671  break;
672  node = node->parent;
673  }
674  return (FIFF_OK);
675 }
676 
677 static int get_evoked_essentials (FiffStream::SPtr& stream, /* This is our file */
678  const FiffDirNode::SPtr& node, /* The interesting node */
679  float& sfreq, /* Sampling frequency
680  * The value pointed by this is not
681  * modified if individual sampling
682  * frequency is found */
683  float& tmin, /* Time scale minimum */
684  int& nsamp, /* Number of samples */
685  int& nave, /* Number of averaged responses */
686  int& akind, /* Aspect type */
687  int *& artefs, /* Artefact removal parameters */
688  int& nartef)
689 /*
690  * Get the essential info for
691  * given evoked response data
692  */
693 {
694  FiffTag::SPtr t_pTag;
695  int k;
696  int to_find = 2;
697  int first = -1;
698  int last = -1;
699  int my_nsamp = -1;
700  float my_tmin = -1;
701  int res = -1;
702  fiff_int_t kind, pos;
703 
704  fiff_byte_t *tempb;
705 
706  FiffDirNode::SPtr tmp_node = node;
707 
708  /*
709  * This is rather difficult...
710  */
711  if (find_between (stream,tmp_node,tmp_node->parent,FIFF_NAVE,&tempb) == FIFF_FAIL)
712  return res;
713  if (tempb)
714  nave = *(int *)tempb;
715  FREE_9(tempb);
716  if (find_between (stream,tmp_node,tmp_node->parent,
717  FIFF_SFREQ,&tempb) == FIFF_FAIL)
718  return res;
719  if (tempb)
720  sfreq = *(float *)tempb;
721  FREE_9(tempb);
722 
723  if (find_between (stream,tmp_node,tmp_node->parent,
724  FIFF_ASPECT_KIND,&tempb) == FIFF_FAIL)
725  return res;
726  if (tempb)
727  akind = *(int *)tempb;
728  else
729  akind = FIFFV_ASPECT_AVERAGE; /* Just a guess */
730  FREE_9(tempb);
731  /*
732  * Find evoked response descriptive data
733  */
734  tmp_node = tmp_node->parent;
735 
736  // tag.data = NULL;
737  for (k = 0; k < tmp_node->dir_tree.size(); k++) {
738  kind = tmp_node->dir_tree[k]->kind;
739  pos = tmp_node->dir_tree[k]->pos;
740  switch (kind) {
741 
742  case FIFF_FIRST_SAMPLE :
743  if (!stream->read_tag(t_pTag,pos))
744  goto out;
745  first = *t_pTag->toInt(); to_find--;
746  break;
747 
748  case FIFF_LAST_SAMPLE :
749  if (!stream->read_tag(t_pTag,pos))
750  goto out;
751  last = *t_pTag->toInt(); to_find--;
752  break;
753 
754  case FIFF_NO_SAMPLES :
755  if (!stream->read_tag(t_pTag,pos))
756  goto out;
757  my_nsamp = *t_pTag->toInt(); to_find--;
758  break;
759 
760  case FIFF_FIRST_TIME :
761  if (!stream->read_tag(t_pTag,pos))
762  goto out;
763  my_tmin = *t_pTag->toFloat(); to_find--;
764  break;
765  case FIFF_ARTEF_REMOVAL :
766  if (!stream->read_tag(t_pTag,pos))
767  goto out;
768  qDebug() << "TODO: check whether artefs contains the right stuff -> use MatrixXi instead";
769  artefs = t_pTag->toInt();
770  nartef = t_pTag->size()/(3*sizeof(int));
771  break;
772  }
773  }
774  if (to_find > 0) {
775  printf ("Not all essential tags were found!");
776  goto out;
777  }
778  if (first != -1 && last != -1) {
779  nsamp = (last)-(first)+1;
780  tmin = (first)/(sfreq);
781  }
782  else if (my_tmin != -1 && my_nsamp != -1) {
783  tmin = my_tmin;
784  nsamp = my_nsamp;
785  }
786  else {
787  printf("Not enough data for time scale definition!");
788  goto out;
789  }
790  res = 0;
791 
792 out : {
793  return res;
794  }
795 }
796 
797 static int get_evoked_optional( FiffStream::SPtr& stream,
798  const FiffDirNode::SPtr& node, /* The directory node containing our data */
799  int *nchan, /* Number of channels */
800  QList<FiffChInfo>& chp) /* Channel descriptions */
801 /*
802  * The channel info may have been modified
803  */
804 {
805  int res = FIFF_FAIL;
806  QList<FiffChInfo> new_ch;
807  int new_nchan = *nchan;
808  int k,to_find;
809  FiffTag::SPtr t_pTag;
810  fiff_int_t kind, pos;
811  FiffChInfo this_ch;
812  FiffDirNode::SPtr evoked_node;
813 
814  if (!(evoked_node = find_evoked(node))) {
815  res = FIFF_OK;
816  goto out;
817  }
818 
819  to_find = 0;
820  if(evoked_node->find_tag(stream, FIFF_NCHAN, t_pTag))
821  new_nchan = *t_pTag->toInt();
822  else
823  new_nchan = *nchan;
824 
825  for (k = 0; k < evoked_node->nent(); k++) {
826  kind = evoked_node->dir[k]->kind;
827  pos = evoked_node->dir[k]->pos;
828  if (kind == FIFF_CH_INFO) { /* Information about one channel */
829  if (new_ch.isEmpty()) {
830  to_find = new_nchan;
831  for (int i = 0; i < to_find; i++) {
832  new_ch.append(FiffChInfo());
833  }
834  }
835  if (!stream->read_tag(t_pTag,pos))
836  goto out;
837 
838  this_ch = t_pTag->toChInfo();
839  if (this_ch.scanNo <= 0 || this_ch.scanNo > new_nchan) {
840  printf ("FIFF_CH_INFO : scan # out of range!");
841  goto out;
842  }
843  else
844  new_ch[this_ch.scanNo-1] = this_ch;
845  to_find--;
846  }
847  }
848  if (to_find != 0) {
849  printf("All channels were not specified "
850  "at the FIFFB_EVOKED level.");
851  goto out;
852  }
853  res = FIFF_OK;
854  goto out;
855 
856 out : {
857  if (res == FIFF_OK) {
858  *nchan = new_nchan;
859  if (!new_ch.isEmpty()) {
860  chp = new_ch;
861  }
862  }
863  return res;
864  }
865 }
866 
867 static void unpack_data(double offset,
868  double scale,
869  short *packed,
870  int nsamp,
871  float *orig)
872 {
873  int k;
874  for (k = 0; k < nsamp; k++)
875  orig[k] = scale * packed[k] + offset;
876  return;
877 }
878 
879 static float **get_epochs ( FiffStream::SPtr& stream, /* This is our file */
880  const FiffDirNode::SPtr& node, /* The interesting node */
881  int nchan, int nsamp) /* Number of channels and number of samples to be expected */
882 /*
883  * Get the evoked response epochs
884  */
885 {
886  fiff_int_t kind, pos;
887  FiffTag::SPtr t_pTag;
888  int k;
889  int ch;
890  float **epochs = NULL;
891  float offset,scale;
892  short *packed;
893 
894  for (k = 0, ch = 0; k < node->nent() && ch < nchan; k++) {
895  kind = node->dir[k]->kind;
896  pos = node->dir[k]->pos;
897  if (kind == FIFF_EPOCH) {
898  if (!stream->read_tag(t_pTag,pos))
899  goto bad;
900  if (t_pTag->type & FIFFT_MATRIX) {
901  if ((t_pTag->type & ~FIFFT_MATRIX) != FIFFT_FLOAT) {
902  printf("Epochs in matrix should be floats!");
903  goto bad;
904  }
905 
906  qint32 ndim;
907  QVector<qint32> dims;
908  t_pTag->getMatrixDimensions(ndim, dims);
909 
910  if (ndim != 2) {
911  printf("Data matrix dimension should be two!");
912  goto bad;
913  }
914  if (dims[0] != nsamp) {
915  printf("Incorrect number of samples in data matrix!");
916  goto bad;
917  }
918  if (dims[1] != nchan) {
919  printf("Incorrect number of channels in data matrix!");
920  goto bad;
921  }
922  MatrixXf tmp_epochs = t_pTag->toFloatMatrix().transpose();
923  epochs = ALLOC_CMATRIX_9(tmp_epochs.rows(),tmp_epochs.cols());
924  fromFloatEigenMatrix_9(tmp_epochs, epochs);
925  ch = nchan;
926  break; /* We have the data */
927  }
928  else { /* Individual epochs */
929  if (epochs == NULL)
930  epochs = ALLOC_CMATRIX_9(nchan,nsamp);
931  if (t_pTag->type == FIFFT_OLD_PACK) {
932  offset = ((float *)t_pTag->data())[0];
933  scale = ((float *)t_pTag->data())[1];
934  packed = (short *)(((float *)t_pTag->data())+2);
935  unpack_data(offset,scale,packed,nsamp,epochs[ch++]);
936  }
937  else if (t_pTag->type == FIFFT_FLOAT)
938  memcpy(epochs[ch++],t_pTag->data(),nsamp*sizeof(float));
939  else {
940  printf ("Unknown data packing type!");
941  FREE_CMATRIX_9(epochs);
942  return (NULL);
943  }
944  }
945  if (ch == nchan)
946  return (epochs);
947  }
948  }
949  if (ch < nchan) {
950  printf ("All epochs were not found!");
951  goto bad;
952  }
953  return (epochs);
954 
955 bad : {
956  FREE_CMATRIX_9(epochs);
957  return (NULL);
958  }
959 }
960 
961 int mne_find_evoked_types_comments ( FiffStream::SPtr& stream,
962  QList<FiffDirNode::SPtr>& nodesp,
963  int **aspect_typesp,
964  QStringList* commentsp)
965 /*
966  * Find all data we are able to process
967  */
968 {
969  QList<FiffDirNode::SPtr> evoked;
970  QList<FiffDirNode::SPtr> meas;
971  QList<FiffDirNode::SPtr> nodes;
972  int evoked_count,count;
973  QString part,type,meas_date;
974  QStringList comments;
975  int *types = NULL;
976  int j,k,p;
977 
978  if (stream == NULL)
979  return 0;
980  /*
981  * First find all measurements
982  */
983  meas = stream->dirtree()->dir_tree_find(FIFFB_MEAS);
984  /*
985  * Process each measurement
986  */
987  for (count = 0,p = 0; p < meas.size(); p++) {
988  evoked = meas[p]->dir_tree_find(FIFFB_EVOKED);
989  /*
990  * Count the entries
991  */
992  for (evoked_count = 0, j = 0; j < evoked.size(); j++) {
993  for (k = 0; k < evoked[j]->nchild(); k++) {
994  if (evoked[j]->children[k]->type == FIFFB_ASPECT) {
995  evoked_count++;
996  }
997  }
998  }
999  /*
1000  * Enlarge tables
1001  */
1002  types = REALLOC_9(types,count+evoked_count+1,int);
1003  /*
1004  * Insert node references and compile associated comments...
1005  */
1006  for (j = 0; j < evoked.size(); j++) /* Evoked data */
1007  for (k = 0; k < evoked[j]->nchild(); k++)
1008  if (evoked[j]->children[k]->type == FIFFB_ASPECT) {
1009  meas_date = get_meas_date(stream,evoked[j]);
1010  part = get_comment(stream,evoked[j]);
1011  get_aspect_name_type(stream,evoked[j]->children[k],type,types+count);
1012  if (!meas_date.isEmpty()) {
1013  comments.append(QString("%1>%2>%3").arg(meas_date).arg(part).arg(type));
1014  }
1015  else {
1016  comments.append(QString("%1>%2").arg(part).arg(type));
1017  }
1018  nodes.append(evoked[j]->children[k]);
1019  count++;
1020  }
1021  }
1022  if (count == 0) { /* Nothing to report */
1023  comments.clear();
1024  nodesp.clear();
1025  if (commentsp)
1026  *commentsp = comments;
1027  if (aspect_typesp)
1028  *aspect_typesp = NULL;
1029  return 0;
1030  }
1031  else { /* Return the appropriate variables */
1032 // comments[count] = NULL;
1033  types[count] = -1;
1034  nodesp = nodes;
1035  if (commentsp)
1036  *commentsp = comments;
1037  else
1038  comments.clear();
1039  if (aspect_typesp)
1040  *aspect_typesp = types;
1041  else
1042  FREE_9(types);
1043  return count;
1044  }
1045 }
1046 
1047 QList<FiffDirNode::SPtr> mne_find_evoked ( FiffStream::SPtr& stream, QStringList* commentsp)
1048 /* Optionally return the compiled comments here */
1049 {
1050  QList<FiffDirNode::SPtr> evoked;
1051  mne_find_evoked_types_comments(stream,evoked,NULL,commentsp);
1052  return evoked;
1053 }
1054 
1055 static void remove_artefacts (float *resp,
1056  int nsamp,
1057  int *artefs,
1058  int nartef)
1059 /*
1060  * Apply the artefact removal
1061  */
1062 {
1063  int start,end;
1064  int j,k;
1065  float a,b;
1066  int remove_jump;
1067 
1068  for (k = 0; k < nartef; k++) {
1069  if (artefs[3*k] == FIFFV_ARTEF_NONE || artefs[3*k] == FIFFV_ARTEF_KEEP)
1070  continue;
1071  remove_jump = (artefs[3*k] == FIFFV_ARTEF_NOJUMP);
1072  /*
1073  * Find out the indices for the start and end times
1074  */
1075  start = artefs[3*k+1];
1076  end = artefs[3*k+2];
1077  start = MAX(0,MIN(start,nsamp));
1078  end = MAX(0,MIN(end,nsamp));
1079  /*
1080  * Replace the artefact region with a straight line
1081  */
1082  if (start < end) {
1083  if (remove_jump) { /* Remove jump... */
1084  a = resp[end] - resp[start];
1085  for (j = 0; j <=start; j++)
1086  resp[j] = resp[j] + a;
1087  for (j = start+1 ; j < end; j++)
1088  resp[j] = resp[end];
1089  }
1090  else { /* Just connect... */
1091  a = (resp[end]-resp[start])/(end-start);
1092  b = (resp[start]*end - resp[end]*start)/(end-start);
1093  for (j = start+1 ; j < end; j++)
1094  resp[j] = a*j+b;
1095  }
1096  }
1097  }
1098  return;
1099 }
1100 
1101 int mne_read_evoked(const QString& name, /* Name of the file */
1102  int setno, /* Which data set */
1103  int *nchanp, /* How many channels */
1104  int *nsampp, /* Number of time points */
1105  float *tminp, /* First time point */
1106  float *sfreqp, /* Sampling frequency */
1107  QList<FiffChInfo>& chsp, /* Channel info (this is now optional as well) */
1108  float ***epochsp, /* Data, channel by channel */
1109  /*
1110  * Optional items follow
1111  */
1112  QString* commentp, /* Comment for these data */
1113  float *highpassp, /* Highpass frequency */
1114  float *lowpassp, /* Lowpass frequency */
1115  int *navep, /* How many averages */
1116  int *aspect_kindp, /* What kind of an evoked data */
1117  FiffCoordTransOld* *transp, /* Coordinate transformation */
1118  fiffId *idp, /* Measurement id */
1119  fiffTime *meas_datep) /* Measurement date */
1120 /*
1121  * Load evoked-response data from a fif file
1122  */
1123 {
1124  QFile file(name);
1125  FiffStream::SPtr stream(new FiffStream(&file));
1126 
1127  QList<FiffDirNode::SPtr> evoked; /* The evoked data nodes */
1128  int nset = 0;
1129  int nchan = 0; /* How many channels */
1130  QStringList comments; /* The associated comments */
1131  float sfreq = 0.0; /* What sampling frequency */
1132  FiffDirNode::SPtr start;
1133  QList<FiffChInfo> chs; /* Channel info */
1134  int *artefs = NULL; /* Artefact limits */
1135  int nartef = 0; /* How many */
1136  float **epochs = NULL; /* The averaged epochs */
1137  FiffCoordTransOld* trans = NULL; /* The coordinate transformation */
1138  fiffId id = NULL; /* Measurement id */
1139  fiffTime meas_date = NULL; /* Measurement date */
1140  int nave = 1; /* Number of averaged responses */
1141  float tmin = 0; /* Time scale minimum */
1142  float lowpass; /* Lowpass filter frequency */
1143  float highpass = 0.0; /* Highpass filter frequency */
1144  int nsamp = 0; /* Samples in epoch */
1145  int aspect_kind; /* What kind of data */
1146  int res = FAIL; /* A little bit of pessimism */
1147 
1148  float *epoch;
1149  int j,k;
1150 
1151  if (setno < 0) {
1152  printf ("Evoked response selector must be positive!");
1153  goto out;
1154  }
1155 
1156  if(!stream->open())
1157  goto out;
1158 
1159  /*
1160  * Select correct data set
1161  */
1162  evoked = mne_find_evoked(stream,(commentp == NULL) ? NULL : &comments);
1163  if (!evoked.size()) {
1164  printf ("No evoked response data available here");
1165  goto out;
1166  }
1167 
1168  nset = evoked.size();
1169 
1170  if (setno < nset) {
1171  start = evoked[setno];
1172  }
1173  else {
1174  printf ("Too few evoked response data sets (how come?)");
1175  goto out;
1176  }
1177  /*
1178  * Get various things...
1179  */
1180  if (get_meas_info (stream,
1181  start,
1182  &id,
1183  &meas_date,
1184  &nchan,
1185  &sfreq,
1186  &highpass,
1187  &lowpass,
1188  chs,
1189  &trans) == -1)
1190  goto out;
1191 
1192  /*
1193  * sfreq is listed here again because
1194  * there might be an individual one in the
1195  * evoked-response data
1196  */
1197  if (get_evoked_essentials(stream,start,sfreq,
1198  tmin,nsamp,nave,aspect_kind,
1199  artefs,nartef) == -1)
1200  goto out;
1201  /*
1202  * Some things may be redefined at a lower level
1203  */
1204  if (get_evoked_optional(stream,
1205  start,
1206  &nchan,
1207  chs) == -1)
1208  goto out;
1209  /*
1210  * Omit nonmagnetic channels
1211  */
1212  if ((epochs = get_epochs(stream,start,nchan,nsamp)) == NULL)
1213  goto out;
1214  /*
1215  * Change artefact limits to start from 0
1216  */
1217  for (k = 0; k < nartef; k++) {
1218  qDebug() << "TODO: Artefact Vectors do not contain the right stuff!";
1219  artefs[2*k+1] = artefs[2*k+1] - sfreq*tmin;
1220  artefs[2*k+2] = artefs[2*k+2] - sfreq*tmin;
1221  }
1222  for (k = 0; k < nchan; k++) {
1223  epoch = epochs[k];
1224  for (j = 0; j < nsamp; j++)
1225  epoch[j] = chs[k].cal*epoch[j];
1226  remove_artefacts(epoch,nsamp,artefs,nartef);
1227  }
1228  /*
1229  * Ready to go
1230  */
1231  chsp = chs;
1232  *tminp = tmin;
1233  *nchanp = nchan;
1234  *nsampp = nsamp;
1235  *sfreqp = sfreq;
1236  *epochsp = epochs; epochs = NULL;
1237  /*
1238  * Fill in the optional data
1239  */
1240  if (commentp) {
1241  *commentp = comments[setno];
1242  comments[setno] = "";
1243  }
1244  if (highpassp)
1245  *highpassp = highpass;
1246  if (lowpassp)
1247  *lowpassp = lowpass;
1248  if (transp) {
1249  *transp = trans;
1250  trans = NULL;
1251  }
1252  if (navep)
1253  *navep = nave;
1254  if (aspect_kindp)
1255  *aspect_kindp = aspect_kind;
1256  if (idp) {
1257  *idp = id;
1258  id = NULL;
1259  }
1260  if (meas_datep) {
1261  *meas_datep = meas_date;
1262  meas_date = NULL;
1263  }
1264  res = OK;
1265  /*
1266  * FREE all allocated data on exit
1267  */
1268 out : {
1269  comments.clear();
1270  FREE_9(artefs);
1271  FREE_9(trans);
1272  FREE_9(id);
1273  FREE_9(meas_date);
1274  FREE_CMATRIX_9(epochs);
1275  stream->close();
1276  return res;
1277  }
1278 }
1279 
1280 //============================= mne_inverse_io.c =============================
1281 
1282 #define MAXBUF 200
1283 
1284 char *mne_format_file_id (fiffId id)
1285 
1286 {
1287  char buf[MAXBUF];
1288  static char s[300];
1289  struct tm *ltime;
1290  time_t secs;
1291 
1292  secs = id->time.secs;
1293  ltime = localtime(&secs);
1294  (void)strftime(buf,MAXBUF,"%c",ltime);
1295 
1296  sprintf(s,"%d.%d 0x%x%x %s",id->version>>16,id->version & 0xFFFF,id->machid[0],id->machid[1],buf);
1297  return s;
1298 }
1299 
1300 int mne_read_meg_comp_eeg_ch_info_9(const QString& name,
1301  QList<FiffChInfo>& megp, /* MEG channels */
1302  int *nmegp,
1303  QList<FiffChInfo>& meg_compp,
1304  int *nmeg_compp,
1305  QList<FiffChInfo>& eegp, /* EEG channels */
1306  int *neegp,
1307  FiffCoordTransOld* *meg_head_t,
1308  fiffId *idp) /* The measurement ID */
1309 /*
1310  * Read the channel information and split it into three arrays,
1311  * one for MEG, one for MEG compensation channels, and one for EEG
1312  */
1313 {
1314  QFile file(name);
1315  FiffStream::SPtr stream(new FiffStream(&file));
1316 
1317  QList<FiffChInfo> chs;
1318  int nchan = 0;
1319  QList<FiffChInfo> meg;
1320  int nmeg = 0;
1321  QList<FiffChInfo> meg_comp;
1322  int nmeg_comp = 0;
1323  QList<FiffChInfo> eeg;
1324  int neeg = 0;
1325  fiffId id = NULL;
1326  QList<FiffDirNode::SPtr> nodes;
1327  FiffDirNode::SPtr info;
1328  FiffTag::SPtr t_pTag;
1329  FiffChInfo this_ch;
1330  FiffCoordTransOld* t = NULL;
1331  fiff_int_t kind, pos;
1332  int j,k,to_find;
1333 
1334  if(!stream->open())
1335  goto bad;
1336 
1337  nodes = stream->dirtree()->dir_tree_find(FIFFB_MNE_PARENT_MEAS_FILE);
1338 
1339  if (nodes.size() == 0) {
1340  nodes = stream->dirtree()->dir_tree_find(FIFFB_MEAS_INFO);
1341  if (nodes.size() == 0) {
1342  qCritical ("Could not find the channel information.");
1343  goto bad;
1344  }
1345  }
1346  info = nodes[0];
1347  to_find = 0;
1348  for (k = 0; k < info->nent(); k++) {
1349  kind = info->dir[k]->kind;
1350  pos = info->dir[k]->pos;
1351  switch (kind) {
1352  case FIFF_NCHAN :
1353  if (!stream->read_tag(t_pTag,pos))
1354  goto bad;
1355  nchan = *t_pTag->toInt();
1356 
1357  for (j = 0; j < nchan; j++) {
1358  chs.append(FiffChInfo());
1359  chs[j].scanNo = -1;
1360  }
1361  to_find = nchan;
1362  break;
1363 
1364  case FIFF_PARENT_BLOCK_ID :
1365  if(!stream->read_tag(t_pTag, pos))
1366  goto bad;
1367  // id = t_pTag->toFiffID();
1368  *id = *(fiffId)t_pTag->data();
1369  break;
1370 
1371  case FIFF_COORD_TRANS :
1372  if(!stream->read_tag(t_pTag, pos))
1373  goto bad;
1374  // t = t_pTag->toCoordTrans();
1375  t = FiffCoordTransOld::read_helper( t_pTag );
1376  if (t->from != FIFFV_COORD_DEVICE || t->to != FIFFV_COORD_HEAD)
1377  t = NULL;
1378  break;
1379 
1380  case FIFF_CH_INFO : /* Information about one channel */
1381  if(!stream->read_tag(t_pTag, pos))
1382  goto bad;
1383 
1384  this_ch = t_pTag->toChInfo();
1385  if (this_ch.scanNo <= 0 || this_ch.scanNo > nchan) {
1386  printf ("FIFF_CH_INFO : scan # out of range %d (%d)!",this_ch.scanNo,nchan);
1387  goto bad;
1388  }
1389  else
1390  chs[this_ch.scanNo-1] = this_ch;
1391  to_find--;
1392  break;
1393  }
1394  }
1395  if (to_find != 0) {
1396  qCritical("Some of the channel information was missing.");
1397  goto bad;
1398  }
1399  if (t == NULL && meg_head_t != NULL) {
1400  /*
1401  * Try again in a more general fashion
1402  */
1403  if ((t = FiffCoordTransOld::mne_read_meas_transform(name)) == NULL) {
1404  qCritical("MEG -> head coordinate transformation not found.");
1405  goto bad;
1406  }
1407  }
1408  /*
1409  * Sort out the channels
1410  */
1411  for (k = 0; k < nchan; k++) {
1412  if (chs[k].kind == FIFFV_MEG_CH) {
1413  meg.append(chs[k]);
1414  nmeg++;
1415  } else if (chs[k].kind == FIFFV_REF_MEG_CH) {
1416  meg_comp.append(chs[k]);
1417  nmeg_comp++;
1418  } else if (chs[k].kind == FIFFV_EEG_CH) {
1419  eeg.append(chs[k]);
1420  neeg++;
1421  }
1422  }
1423  // fiff_close(in);
1424  stream->close();
1425 
1426  megp = meg;
1427  if(nmegp) {
1428  *nmegp = nmeg;
1429  }
1430 
1431  meg_compp = meg_comp;
1432  if(nmeg_compp) {
1433  *nmeg_compp = nmeg_comp;
1434  }
1435  eegp = eeg;
1436  if(neegp) {
1437  *neegp = neeg;
1438  }
1439 
1440  if (idp == NULL) {
1441  FREE_9(id);
1442  }
1443  else
1444  *idp = id;
1445  if (meg_head_t == NULL) {
1446  FREE_9(t);
1447  }
1448  else
1449  *meg_head_t = t;
1450 
1451  return FIFF_OK;
1452 
1453 bad : {
1454  // fiff_close(in);
1455  stream->close();
1456  FREE_9(id);
1457  // FREE(tag.data);
1458  FREE_9(t);
1459  return FIFF_FAIL;
1460  }
1461 }
1462 
1463 int mne_read_bad_channel_list_from_node_9( FiffStream::SPtr& stream,
1464  const FiffDirNode::SPtr& pNode, QStringList& listp, int& nlistp)
1465 {
1466  FiffDirNode::SPtr node,bad;
1467  QList<FiffDirNode::SPtr> temp;
1468  QStringList list;
1469  int nlist = 0;
1470  FiffTag::SPtr t_pTag;
1471  QString names;
1472 
1473  if (pNode->isEmpty())
1474  node = stream->dirtree();
1475  else
1476  node = pNode;
1477 
1478  temp = node->dir_tree_find(FIFFB_MNE_BAD_CHANNELS);
1479  if (temp.size() > 0) {
1480  bad = temp[0];
1481 
1482  bad->find_tag(stream, FIFF_MNE_CH_NAME_LIST, t_pTag);
1483  if (t_pTag) {
1484  names = t_pTag->toString();
1485  mne_string_to_name_list_9(names,list,nlist);
1486  }
1487  }
1488  listp = list;
1489  nlistp = nlist;
1490  return OK;
1491 }
1492 
1493 int mne_read_bad_channel_list_9(const QString& name, QStringList& listp, int& nlistp)
1494 
1495 {
1496  QFile file(name);
1497  FiffStream::SPtr stream(new FiffStream(&file));
1498 
1499  int res;
1500 
1501  if(!stream->open())
1502  return FAIL;
1503 
1504  res = mne_read_bad_channel_list_from_node_9(stream,stream->dirtree(),listp,nlistp);
1505 
1506  stream->close();
1507 
1508  return res;
1509 }
1510 
1511 //=============================================================================================================
1512 // DEFINE MEMBER METHODS
1513 //=============================================================================================================
1514 
1515 MneMeasData::MneMeasData()
1516  :meas_id (NULL)
1517  ,current (NULL)
1518  ,ch_major (FALSE)
1519  ,nset (0)
1520  ,nchan (0)
1521  ,op (NULL)
1522  ,fwd (NULL)
1523  ,meg_head_t(NULL)
1524  ,mri_head_t(NULL)
1525  ,proj (NULL)
1526  ,comp (NULL)
1527  ,raw (NULL)
1528  ,chsel (NULL)
1529  ,bad (NULL)
1530  ,nbad (0)
1531  ,badlist (NULL)
1532 {
1533  meas_date.secs = 0;
1534  meas_date.usecs = 0;
1535 }
1536 
1537 //=============================================================================================================
1538 
1539 MneMeasData::~MneMeasData()
1540 {
1541  int k;
1542 
1543  filename.clear();
1544  FREE_9(meas_id);
1545  if(meg_head_t)
1546  delete meg_head_t;
1547  if(mri_head_t)
1548  delete mri_head_t;
1549  if(proj)
1550  delete proj;
1551  if(comp)
1552  delete comp;
1553  FREE_9(bad);
1554  badlist.clear();
1555 
1556  for (k = 0; k < nset; k++)
1557  delete sets[k];
1558 
1559  if(raw)
1560  delete raw;
1561  mne_ch_selection_free_9(chsel);
1562 
1563  return;
1564 }
1565 
1566 //=============================================================================================================
1567 
1568 void MneMeasData::adjust_baselines(float bmin, float bmax)
1569 {
1570  int b1,b2;
1571  float sfreq,tmin,tmax;
1572  float **data;
1573  float ave;
1574  int s,c;
1575 
1576  if (!this->current)
1577  return;
1578 
1579  sfreq = 1.0/ this->current->tstep;
1580  tmin = this->current->tmin;
1581  tmax = this->current->tmin + ( this->current->np-1)/sfreq;
1582 
1583  if (bmin < tmin)
1584  b1 = 0;
1585  else if (bmin > tmax)
1586  b1 = this->current->np;
1587  else {
1588  for (b1 = 0; b1/sfreq + tmin < bmin; b1++)
1589  ;
1590  if (b1 < 0)
1591  b1 = 0;
1592  else if (b1 > this->current->np)
1593  b1 = this->current->np;
1594  }
1595  if (bmax < tmin)
1596  b2 = 0;
1597  else if (bmax > tmax)
1598  b2 = this->current->np;
1599  else {
1600  for (b2 = this->current->np; b2/sfreq + tmin > bmax; b2--)
1601  ;
1602  if (b2 < 0)
1603  b2 = 0;
1604  else if (b2 > this->current->np)
1605  b2 = this->current->np;
1606  }
1607  data = this->current->data;
1608  if (b2 > b1) {
1609  for (c = 0; c < this->nchan; c++) {
1610  for (s = b1, ave = 0.0; s < b2; s++)
1611  ave += data[s][c];
1612  ave = ave/(b2-b1);
1613  this->current->baselines[c] += ave;
1614  for (s = 0; s < this->current->np; s++)
1615  data[s][c] = data[s][c] - ave;
1616  }
1617  qDebug() << "TODO: Check comments content";
1618  printf("\t%s : using baseline %7.1f ... %7.1f ms\n",
1619  this->current->comment.toUtf8().constData() ? this->current->comment.toUtf8().constData() : "unknown",
1620  1000*(tmin+b1/sfreq),
1621  1000*(tmin+b2/sfreq));
1622  }
1623  return;
1624 }
1625 
1626 //=============================================================================================================
1627 
1628 MneMeasData *MneMeasData::mne_read_meas_data_add(const QString &name,
1629  int set,
1630  MneInverseOperator* op,
1631  MneNamedMatrix *fwd,
1632  const QStringList& namesp,
1633  int nnamesp,
1634  MneMeasData *add_to) /* Add to this */
1635 /*
1636  * Read an evoked-response data file
1637  */
1638 {
1639  /*
1640  * Data read from the file
1641  */
1642  QList<FiffChInfo> chs;
1643  int nchan_file,nsamp;
1644  float dtmin,dtmax,sfreq;
1645  QString comment;
1646  float **data = NULL;
1647  float lowpass,highpass;
1648  int nave;
1649  int aspect_kind;
1650  fiffId id = NULL;
1651  FiffCoordTransOld* t = NULL;
1652  fiffTime meas_date = NULL;
1653  QString stim14_name;
1654  /*
1655  * Desired channels
1656  */
1657  QStringList names;
1658  int nchan = 0;
1659  /*
1660  * Selected channels
1661  */
1662  int *sel = NULL;
1663  int stim14 = -1;
1664  /*
1665  * Other stuff
1666  */
1667  float *source,tmin,tmax;
1668  int k,p,c,np,n1,n2;
1669  MneMeasData* res = NULL;
1670  MneMeasData* new_data = add_to;
1671  MneMeasDataSet* dataset = NULL;
1672 
1673  stim14_name = getenv(MNE_ENV_TRIGGER_CH);
1674  if (stim14_name.isEmpty() || stim14_name.size() == 0)
1675  stim14_name = MNE_DEFAULT_TRIGGER_CH;
1676 
1677  if (add_to)
1678  mne_channel_names_to_name_list_9(add_to->chs,add_to->nchan,names,nchan);
1679  else {
1680  if (op) {
1681  names = op->eigen_fields->collist;
1682  nchan = op->nchan;
1683  }
1684  else if (fwd) {
1685  names = fwd->collist;
1686  nchan = fwd->ncol;
1687  }
1688  else {
1689  names = namesp;
1690  nchan = nnamesp;
1691  }
1692  if (names.isEmpty())
1693  nchan = 0;
1694  }
1695  /*
1696  * Read the evoked data file
1697  */
1698  if (mne_read_evoked(name,set-1,
1699  &nchan_file,
1700  &nsamp,
1701  &dtmin,
1702  &sfreq,
1703  chs,
1704  &data,
1705  &comment,
1706  &highpass,
1707  &lowpass,
1708  &nave,
1709  &aspect_kind,
1710  &t,
1711  &id,
1712  &meas_date) == FAIL)
1713  goto out;
1714 
1715  if (id)
1716  printf("\tMeasurement file id: %s\n",mne_format_file_id(id));
1717 
1718 #ifdef FOO
1719  if (add_to) { /* Should add consistency check here */
1720  printf("\tWarning: data set consistency check is still in the works.\n");
1721  }
1722 #endif
1723  /*
1724  * Pick out the necessary channels
1725  */
1726  if (nchan > 0) {
1727  sel = MALLOC_9(nchan,int);
1728  for (k = 0; k < nchan; k++)
1729  sel[k] = -1;
1730  for (c = 0; c < nchan_file; c++) {
1731  for (k = 0; k < nchan; k++) {
1732  if (sel[k] == -1 && QString::compare(chs[c].ch_name,names[k]) == 0) {
1733  sel[k] = c;
1734  break;
1735  }
1736  }
1737  if (QString::compare(stim14_name,chs[c].ch_name) == 0) {
1738  stim14 = c;
1739  }
1740  }
1741  for (k = 0; k < nchan; k++)
1742  if (sel[k] == -1) {
1743  printf("All channels needed were not in the MEG/EEG data file "
1744  "(first missing: %s).",names[k].toUtf8().constData());
1745  goto out;
1746  }
1747  }
1748  else { /* Load all channels */
1749  sel = MALLOC_9(nchan_file,int);
1750  for (c = 0, nchan = 0; c < nchan_file; c++) {
1751  if (chs[c].kind == FIFFV_MEG_CH || chs[c].kind == FIFFV_EEG_CH) {
1752  sel[nchan] = c;
1753  nchan++;
1754  }
1755  if (QString::compare(stim14_name,chs[c].ch_name) == 0) {
1756  stim14 = c;
1757  }
1758  }
1759  }
1760  /*
1761  * Cut the data to the analysis time range
1762  */
1763  n1 = 0;
1764  n2 = nsamp;
1765  np = n2 - n1;
1766  dtmax = dtmin + (np-1)/sfreq;
1767  /*
1768  * Then the analysis time range
1769  */
1770  tmin = dtmin;
1771  tmax = dtmax;
1772  printf("\tData time range: %8.1f ... %8.1f ms\n",1000*tmin,1000*tmax);
1773  /*
1774  * Just put it together
1775  */
1776  if (!new_data) { /* We need a new meas data structure */
1777  new_data = new MneMeasData;
1778  new_data->filename = name;
1779  new_data->meas_id = id; id = NULL;
1780  /*
1781  * Getting starting time from measurement ID is not too accurate...
1782  */
1783  if (meas_date)
1784  new_data->meas_date = *meas_date;
1785  else {
1786  if (new_data->meas_id)
1787  new_data->meas_date = new_data->meas_id->time;
1788  else {
1789  new_data->meas_date.secs = 0;
1790  new_data->meas_date.usecs = 0;
1791  }
1792  }
1793  new_data->lowpass = lowpass;
1794  new_data->highpass = highpass;
1795  new_data->nchan = nchan;
1796  new_data->sfreq = sfreq;
1797 
1798  if (t) {
1799  new_data->meg_head_t = t;
1800  t = NULL;
1801  printf("\tUsing MEG <-> head transform from the present data set\n");
1802  }
1803  if (op != NULL && op->mri_head_t != NULL) { /* Copy if available */
1804  if (!new_data->mri_head_t)
1805  new_data->mri_head_t = new FiffCoordTransOld;
1806  *(new_data->mri_head_t) = *(op->mri_head_t);
1807  printf("\tPicked MRI <-> head transform from the inverse operator\n");
1808  }
1809  /*
1810  * Channel list
1811  */
1812  for (k = 0; k < nchan; k++) {
1813  new_data->chs.append(FiffChInfo());
1814  new_data->chs[k] = chs[sel[k]];
1815  }
1816 
1817  new_data->op = op; /* Attach inverse operator */
1818  new_data->fwd = fwd; /* ...or a fwd operator */
1819  if (op) /* Attach the projection operator and CTF compensation info to the data, too */
1820  new_data->proj = MneProjOp::mne_dup_proj_op(op->proj);
1821  else {
1822  new_data->proj = MneProjOp::mne_read_proj_op(name);
1823  if (new_data->proj && new_data->proj->nitems > 0) {
1824  printf("\tLoaded projection from %s:\n",name.toUtf8().data());
1825  MneProjOp::mne_proj_op_report(stderr,"\t\t",new_data->proj);
1826  }
1827  new_data->comp = MneCTFCompDataSet::mne_read_ctf_comp_data(name);
1828  if (new_data->comp == NULL)
1829  goto out;
1830  if (new_data->comp->ncomp > 0)
1831  printf("\tRead %d compensation data sets from %s\n",new_data->comp->ncomp,name.toUtf8().data());
1832  }
1833  /*
1834  * Th bad channel stuff
1835  */
1836  {
1837  int b;
1838 
1839  new_data->bad = MALLOC_9(new_data->nchan,int);
1840  for (k = 0; k < new_data->nchan; k++)
1841  new_data->bad[k] = FALSE;
1842 
1843  if (mne_read_bad_channel_list_9(name,new_data->badlist,new_data->nbad) == OK) {
1844  for (b = 0; b < new_data->nbad; b++) {
1845  for (k = 0; k < new_data->nchan; k++) {
1846  if (QString::compare(new_data->chs[k].ch_name,new_data->badlist[b],Qt::CaseInsensitive) == 0) {
1847  new_data->bad[k] = TRUE;
1848  break;
1849  }
1850  }
1851  }
1852  printf("\t%d bad channels read from %s%s",new_data->nbad,name.toUtf8().data(),new_data->nbad > 0 ? ":\n" : "\n");
1853  if (new_data->nbad > 0) {
1854  printf("\t\t");
1855  for (k = 0; k < new_data->nbad; k++)
1856  printf("%s%c",new_data->badlist[k].toUtf8().constData(),k < new_data->nbad-1 ? ' ' : '\n');
1857  }
1858  }
1859  }
1860  }
1861  /*
1862  * New data set is created anyway
1863  */
1864  dataset = new MneMeasDataSet;
1865  dataset->tmin = tmin;
1866  dataset->tstep = 1.0/sfreq;
1867  dataset->first = n1;
1868  dataset->np = np;
1869  dataset->nave = nave;
1870  dataset->kind = aspect_kind;
1871  dataset->data = ALLOC_CMATRIX_9(np,nchan);
1872  dataset->comment = comment; comment.clear();
1873  dataset->baselines = MALLOC_9(nchan,float);
1874  /*
1875  * Pick data from all channels
1876  */
1877  for (k = 0; k < nchan; k++) {
1878  source = data[sel[k]];
1879  /*
1880  * Shift the response
1881  */
1882  for (p = 0; p < np; p++)
1883  dataset->data[p][k] = source[p+n1];
1884  dataset->baselines[k] = 0.0;
1885  }
1886  /*
1887  * Pick the digital trigger channel, too
1888  */
1889  if (stim14 >= 0) {
1890  dataset->stim14 = MALLOC_9(np,float);
1891  source = data[stim14];
1892  for (p = 0; p < np; p++) /* Copy the data and correct for the possible non-unit calibration */
1893  dataset->stim14[p] = source[p+n1]/chs[stim14].cal;
1894  }
1895  new_data->sets.append(dataset); dataset = NULL;
1896  new_data->nset++;
1897  if (!add_to)
1898  new_data->current = new_data->sets[0];
1899  res = new_data;
1900  printf("\t%s dataset %s from %s\n",
1901  add_to ? "Added" : "Loaded",
1902  new_data->sets[new_data->nset-1]->comment.toUtf8().constData() ? new_data->sets[new_data->nset-1]->comment.toUtf8().constData() : "unknown",name.toUtf8().data());
1903 
1904 out : {
1905  FREE_9(sel);
1906  comment.clear();
1907  FREE_CMATRIX_9(data);
1908  FREE_9(t);
1909  FREE_9(id);
1910  if (res == NULL && !add_to)
1911  delete new_data;
1912  if (add_to)
1913  names.clear();
1914  return res;
1915  }
1916 }
1917 
1918 //=============================================================================================================
1919 
1920 MneMeasData *MneMeasData::mne_read_meas_data(const QString &name,
1921  int set,
1922  MneInverseOperator* op,
1923  MneNamedMatrix *fwd,
1924  const QStringList& namesp,
1925  int nnamesp)
1926 
1927 {
1928  return mne_read_meas_data_add(name,set,op,fwd,namesp,nnamesp,NULL);
1929 }
FIFFV_ASPECT_SINGLE
#define FIFFV_ASPECT_SINGLE
Definition: fiff_file.h:440
FIFFLIB::_fiffIdRec
ToDo Old implementation use new fiff_id.h instead.
Definition: fiff_types_mne-c.h:214
FIFFLIB::_fiffIdRec::time
fiffTimeRec time
Definition: fiff_types_mne-c.h:217
MNELIB::mneSparseNamedMatrix
Definition: mne_types.h:316
FIFFV_ASPECT_POWER_DENSITY
#define FIFFV_ASPECT_POWER_DENSITY
Definition: fiff_file.h:444
FIFF_CH_INFO
#define FIFF_CH_INFO
Definition: fiff_file.h:456
FIFFLIB::FiffCoordTransOld::to
FIFFLIB::fiff_int_t to
Definition: fiff_coord_trans_old.h:196
FIFF_EPOCH
#define FIFF_EPOCH
Definition: fiff_file.h:558
FIFFLIB::FiffStream::SPtr
QSharedPointer< FiffStream > SPtr
Definition: fiff_stream.h:107
FIFF_LOWPASS
#define FIFF_LOWPASS
Definition: fiff_file.h:472
FIFF_ASPECT_KIND
#define FIFF_ASPECT_KIND
Definition: fiff_file.h:463
FIFFLIB::_fiffDirEntryRec
ToDo Old implementation use new fiff_dir_entry.h instead.
Definition: fiff_types_mne-c.h:389
FIFFLIB::FiffDirNode::SPtr
QSharedPointer< FiffDirNode > SPtr
Definition: fiff_dir_node.h:76
FIFFV_ASPECT_SAMPLE
#define FIFFV_ASPECT_SAMPLE
Definition: fiff_file.h:443
INVERSELIB::MneMeasDataSet
One data set, used in mneMeasData.
Definition: mne_meas_data_set.h:80
FIFF_LAST_SAMPLE
#define FIFF_LAST_SAMPLE
Definition: fiff_file.h:462
FIFFLIB::_fiffTimeRec::usecs
fiff_int_t usecs
Definition: fiff_types_mne-c.h:184
FIFFLIB::FiffCoordTransOld::from
FIFFLIB::fiff_int_t from
Definition: fiff_coord_trans_old.h:195
FIFF_FIRST_SAMPLE
#define FIFF_FIRST_SAMPLE
Definition: fiff_file.h:461
FIFFLIB::_fiffTimeRec
Definition: fiff_types_mne-c.h:182
mne_inverse_operator.h
MNE Inverse Operator (MneInverseOperator) class declaration.
mne_meas_data.h
MNE Meas Data (MneMeasData) class declaration.
k
int k
Definition: fiff_tag.cpp:322
FIFFV_REF_MEG_CH
#define FIFFV_REF_MEG_CH
Definition: fiff_constants.h:270
FIFFLIB::FiffCoordTransOld
Coordinate transformation descriptor.
Definition: fiff_coord_trans_old.h:80
FIFF_NAVE
#define FIFF_NAVE
Definition: fiff_file.h:460
mne_meas_data_set.h
MNE Meas Data Set (MneMeasDataSet) class declaration.
FIFF_MEAS_DATE
#define FIFF_MEAS_DATE
Definition: fiff_file.h:457
FIFFLIB::_fiffIdRec::machid
fiff_int_t machid[2]
Definition: fiff_types_mne-c.h:216
FIFF_COMMENT
#define FIFF_COMMENT
Definition: fiff_file.h:459
INVERSELIB::ringBufBuf_9
Definition: mne_meas_data.cpp:155
FIFFLIB::FiffChInfo::scanNo
fiff_int_t scanNo
Definition: fiff_ch_info.h:119
FIFF_COORD_TRANS
#define FIFF_COORD_TRANS
Definition: fiff_file.h:475
INVERSELIB::MneMeasData::MneMeasData
MneMeasData()
Definition: mne_meas_data.cpp:1515
FIFFLIB::_fiffIdRec::version
fiff_int_t version
Definition: fiff_types_mne-c.h:215
FIFFLIB::FiffStream::split_name_list
static QStringList split_name_list(QString p_sNameList)
Definition: fiff_stream.cpp:1914
MNELIB::mneChSelection
Definition: mne_types.h:538
FIFFV_ASPECT_ALTAVERAGE
#define FIFFV_ASPECT_ALTAVERAGE
Definition: fiff_file.h:442
FIFFLIB::FiffTag::SPtr
QSharedPointer< FiffTag > SPtr
Definition: fiff_tag.h:152
FIFFV_ASPECT_STD_ERR
#define FIFFV_ASPECT_STD_ERR
Definition: fiff_file.h:439
FIFF_FIRST_TIME
#define FIFF_FIRST_TIME
Definition: fiff_file.h:482
FIFFLIB::FiffChInfo
Channel info descriptor.
Definition: fiff_ch_info.h:74
INVERSELIB::MneMeasData::adjust_baselines
void adjust_baselines(float bmin, float bmax)
Definition: mne_meas_data.cpp:1568
MNELIB::MneNamedMatrix
Matrix specification with a channel list.
Definition: mne_named_matrix.h:84
MNELIB::mneEventList
Definition: mne_types.h:569
FIFFLIB::_fiffTimeRec::secs
fiff_int_t secs
Definition: fiff_types_mne-c.h:183
FIFF_NO_SAMPLES
#define FIFF_NO_SAMPLES
Definition: fiff_file.h:481
INVERSELIB::MneMeasData
easurement data representation in MNE calculations
Definition: mne_meas_data.h:92
FIFF_NCHAN
#define FIFF_NCHAN
Definition: fiff_file.h:453
FIFFV_ASPECT_AVERAGE
#define FIFFV_ASPECT_AVERAGE
Definition: fiff_file.h:438
INVERSELIB::MneInverseOperator
An inverse operator.
Definition: mne_inverse_operator.h:93
FIFFLIB::FiffStream
FIFF File I/O routines.
Definition: fiff_stream.h:104
FIFF_ARTEF_REMOVAL
#define FIFF_ARTEF_REMOVAL
Definition: fiff_file.h:474
mne_named_matrix.h
MNE Named Matrix (MneNamedMatrix) class declaration.
fiff_types.h
Definitions for describing the objects in a FIFF file.
INVERSELIB::MneMeasData::~MneMeasData
~MneMeasData()
Definition: mne_meas_data.cpp:1539
FIFF_SFREQ
#define FIFF_SFREQ
Definition: fiff_file.h:454
FIFF_HIGHPASS
#define FIFF_HIGHPASS
Definition: fiff_file.h:476
INVERSELIB::ringBuf_9
Definition: mne_meas_data.cpp:161
FIFFV_ASPECT_DIPOLE_WAVE
#define FIFFV_ASPECT_DIPOLE_WAVE
Definition: fiff_file.h:445
MNELIB::mneEvent
Definition: mne_types.h:560
Generated on Wed Aug 9 2023 17:55:41 for MNE-CPP by   doxygen 1.8.17