MNE-CPP  0.1.9
A Framework for Electrophysiology
mne_proj_op.cpp
Go to the documentation of this file.
1 //=============================================================================================================
37 //=============================================================================================================
38 // INCLUDES
39 //=============================================================================================================
40 
41 #include <fiff/fiff_constants.h>
42 #include <fiff/fiff_tag.h>
43 
44 #include "mne_proj_op.h"
45 #include "mne_proj_item.h"
46 
47 #include <QFile>
48 
49 #include <Eigen/Core>
50 
51 #ifndef TRUE
52 #define TRUE 1
53 #endif
54 
55 #ifndef FALSE
56 #define FALSE 0
57 #endif
58 
59 #ifndef FAIL
60 #define FAIL -1
61 #endif
62 
63 #ifndef OK
64 #define OK 0
65 #endif
66 
67 #define MALLOC_23(x,t) (t *)malloc((x)*sizeof(t))
68 
69 #define REALLOC_23(x,y,t) (t *)((x == NULL) ? malloc((y)*sizeof(t)) : realloc((x),(y)*sizeof(t)))
70 
71 #define FREE_23(x) if ((char *)(x) != NULL) free((char *)(x))
72 
73 #define FREE_CMATRIX_23(m) mne_free_cmatrix_23((m))
74 
75 void mne_free_cmatrix_23 (float **m)
76 {
77  if (m) {
78  FREE_23(*m);
79  FREE_23(m);
80  }
81 }
82 
83 #define ALLOC_CMATRIX_23(x,y) mne_cmatrix_23((x),(y))
84 
85 static void matrix_error_23(int kind, int nr, int nc)
86 
87 {
88  if (kind == 1)
89  printf("Failed to allocate memory pointers for a %d x %d matrix\n",nr,nc);
90  else if (kind == 2)
91  printf("Failed to allocate memory for a %d x %d matrix\n",nr,nc);
92  else
93  printf("Allocation error for a %d x %d matrix\n",nr,nc);
94  if (sizeof(void *) == 4) {
95  printf("This is probably because you seem to be using a computer with 32-bit architecture.\n");
96  printf("Please consider moving to a 64-bit platform.");
97  }
98  printf("Cannot continue. Sorry.\n");
99  exit(1);
100 }
101 
102 float **mne_cmatrix_23(int nr,int nc)
103 
104 {
105  int i;
106  float **m;
107  float *whole;
108 
109  m = MALLOC_23(nr,float *);
110  if (!m) matrix_error_23(1,nr,nc);
111  whole = MALLOC_23(nr*nc,float);
112  if (!whole) matrix_error_23(2,nr,nc);
113 
114  for(i=0;i<nr;i++)
115  m[i] = whole + i*nc;
116  return m;
117 }
118 
119 float mne_dot_vectors_23 (float *v1,
120  float *v2,
121  int nn)
122 
123 {
124 #ifdef BLAS
125  int one = 1;
126  float res = sdot(&nn,v1,&one,v2,&one);
127  return res;
128 #else
129  float res = 0.0;
130  int k;
131 
132  for (k = 0; k < nn; k++)
133  res = res + v1[k]*v2[k];
134  return res;
135 #endif
136 }
137 
138 //============================= mne_named_matrix.c =============================
139 
140 void mne_string_to_name_list_23(const QString& s, QStringList& listp,int &nlistp)
141 /*
142  * Convert a colon-separated list into a string array
143  */
144 {
145  QStringList list;
146 
147  if (!s.isEmpty() && s.size() > 0) {
149  //list = s.split(":");
150  }
151  listp = list;
152  nlistp = list.size();
153  return;
154 }
155 
156 void fromFloatEigenMatrix_23(const Eigen::MatrixXf& from_mat, float **& to_mat, const int m, const int n)
157 {
158  for ( int i = 0; i < m; ++i)
159  for ( int j = 0; j < n; ++j)
160  to_mat[i][j] = from_mat(i,j);
161 }
162 
163 void fromFloatEigenMatrix_23(const Eigen::MatrixXf& from_mat, float **& to_mat)
164 {
165  fromFloatEigenMatrix_23(from_mat, to_mat, from_mat.rows(), from_mat.cols());
166 }
167 
168 QString mne_name_list_to_string_23(const QStringList& list)
169 /*
170  * Convert a string array to a colon-separated string
171  */
172 {
173  int nlist = list.size();
174  QString res;
175  if (nlist == 0 || list.isEmpty())
176  return res;
177 // res[0] = '\0';
178  for (int k = 0; k < nlist-1; k++) {
179  res += list[k];
180  res += ":";
181  }
182  res += list[nlist-1];
183  return res;
184 }
185 
186 QString mne_channel_names_to_string_23(const QList<FIFFLIB::FiffChInfo>& chs, int nch)
187 /*
188  * Make a colon-separated string out of channel names
189  */
190 {
191  QStringList names;
192  QString res;
193  if (nch <= 0)
194  return res;
195  for (int k = 0; k < nch; k++)
196  names.append(chs.at(k).ch_name);
197  res = mne_name_list_to_string_23(names);
198  return res;
199 }
200 
201 //=============================================================================================================
202 // USED NAMESPACES
203 //=============================================================================================================
204 
205 using namespace Eigen;
206 using namespace FIFFLIB;
207 using namespace MNELIB;
208 
209 //=============================================================================================================
210 // DEFINE MEMBER METHODS
211 //=============================================================================================================
212 
213 MneProjOp::MneProjOp()
214 : nitems (0)
215 , names (NULL)
216 , nch (0)
217 , nvec (0)
218 , proj_data (NULL)
219 {
220 }
221 
222 //=============================================================================================================
223 
225 {
226  // mne_free_proj_op
227  for (int k = 0; k < nitems; k++)
228  if(items[k])
229  delete items[k];
230 
231  // mne_free_proj_op_proj
232 }
233 
234 //=============================================================================================================
235 
237 {
238  if (op == NULL)
239  return;
240 
241  FREE_CMATRIX_23(op->proj_data);
242 
243  op->names.clear();
244  op->nch = 0;
245  op->nvec = 0;
246  op->proj_data = NULL;
247 
248  return;
249 }
250 
251 //=============================================================================================================
252 
253 MneProjOp *MneProjOp::mne_proj_op_combine(MneProjOp *to, MneProjOp *from)
254 /*
255  * Copy items from 'from' operator to 'to' operator
256  */
257 {
258  int k;
259  MneProjItem* it;
260 
261  if (to == NULL)
262  to = new MneProjOp();
263  if (from) {
264  for (k = 0; k < from->nitems; k++) {
265  it = from->items[k];
266  mne_proj_op_add_item(to,it->vecs,it->kind,it->desc);
267  to->items[to->nitems-1]->active_file = it->active_file;
268  }
269  }
270  return to;
271 }
272 
273 //=============================================================================================================
274 
275 void MneProjOp::mne_proj_op_add_item_act(MneProjOp *op, MneNamedMatrix *vecs, int kind, const QString& desc, int is_active)
276 /*
277  * Add a new item to an existing projection operator
278  */
279 {
280  MneProjItem* new_item;
281  int k;
282 
283  // op->items = REALLOC(op->items,op->nitems+1,mneProjItem);
284  // op->items[op->nitems] = new_item = new MneProjItem();
285 
286  new_item = new MneProjItem();
287  op->items.append(new_item);
288 
289  new_item->active = is_active;
290  new_item->vecs = new MneNamedMatrix(*vecs);
291 
292  if (kind == FIFFV_MNE_PROJ_ITEM_EEG_AVREF) {
293  new_item->has_meg = FALSE;
294  new_item->has_eeg = TRUE;
295  }
296  else {
297  for (k = 0; k < vecs->ncol; k++) {
298  if (vecs->collist[k].contains("EEG"))//strstr(vecs->collist[k],"EEG") == vecs->collist[k])
299  new_item->has_eeg = TRUE;
300  if (vecs->collist[k].contains("MEG"))//strstr(vecs->collist[k],"MEG") == vecs->collist[k])
301  new_item->has_meg = TRUE;
302  }
303  if (!new_item->has_meg && !new_item->has_eeg) {
304  new_item->has_meg = TRUE;
305  new_item->has_eeg = FALSE;
306  }
307  else if (new_item->has_meg && new_item->has_eeg) {
308  new_item->has_meg = TRUE;
309  new_item->has_eeg = FALSE;
310  }
311  }
312  if (!desc.isEmpty())
313  new_item->desc = desc;
314  new_item->kind = kind;
315  new_item->nvec = new_item->vecs->nrow;
316 
317  op->nitems++;
318 
319  MneProjOp::mne_free_proj_op_proj(op); /* These data are not valid any more */
320  return;
321 }
322 
323 //=============================================================================================================
324 
325 void MneProjOp::mne_proj_op_add_item(MneProjOp *op, MneNamedMatrix *vecs, int kind, const QString& desc)
326 {
327  mne_proj_op_add_item_act(op, vecs, kind, desc, TRUE);
328 }
329 
330 //=============================================================================================================
331 
332 MneProjOp *MneProjOp::mne_dup_proj_op(MneProjOp *op)
333 /*
334  * Provide a duplicate (item data only)
335  */
336 {
337  MneProjOp* dup = new MneProjOp();
338  MneProjItem* it;
339  int k;
340 
341  if (!op)
342  return NULL;
343 
344  for (k = 0; k < op->nitems; k++) {
345  it = op->items[k];
346  mne_proj_op_add_item_act(dup,it->vecs,it->kind,it->desc,it->active);
347  dup->items[k]->active_file = it->active_file;
348  }
349  return dup;
350 }
351 
352 //=============================================================================================================
353 
354 MneProjOp *MneProjOp::mne_proj_op_average_eeg_ref(const QList<FiffChInfo>& chs, int nch)
355 /*
356  * Make the projection operator for average electrode reference
357  */
358 {
359  int eegcount = 0;
360  int k;
361  float **vec_data;
362  QStringList names;
363  MneNamedMatrix* vecs;
364  MneProjOp* op;
365 
366  for (k = 0; k < nch; k++)
367  if (chs.at(k).kind == FIFFV_EEG_CH)
368  eegcount++;
369  if (eegcount == 0) {
370  qCritical("No EEG channels specified for average reference.");
371  return NULL;
372  }
373 
374  vec_data = ALLOC_CMATRIX_23(1,eegcount);
375 
376  for (k = 0; k < nch; k++)
377  if (chs.at(k).kind == FIFFV_EEG_CH)
378  names.append(chs.at(k).ch_name);
379 
380  for (k = 0; k < eegcount; k++)
381  vec_data[0][k] = 1.0/sqrt((double)eegcount);
382 
383  QStringList emptyList;
384  vecs = MneNamedMatrix::build_named_matrix(1,eegcount,emptyList,names,vec_data);
385 
386  op = new MneProjOp();
387  mne_proj_op_add_item(op,vecs,FIFFV_MNE_PROJ_ITEM_EEG_AVREF,"Average EEG reference");
388 
389  return op;
390 }
391 
392 //=============================================================================================================
393 
394 int MneProjOp::mne_proj_op_affect(MneProjOp *op, const QStringList& list, int nlist)
395 {
396  int k;
397  int naff;
398 
399  if (!op)
400  return 0;
401 
402  for (k = 0, naff = 0; k < op->nitems; k++)
403  if (op->items[k]->active && MneProjItem::mne_proj_item_affect(op->items[k],list,nlist))
404  naff += op->items[k]->nvec;
405 
406  return naff;
407 }
408 
409 //=============================================================================================================
410 
411 int MneProjOp::mne_proj_op_affect_chs(MneProjOp *op, const QList<FiffChInfo>& chs, int nch)
412 {
413  QString ch_string;
414  int res;
415  QStringList list;
416  int nlist;
417 
418  if (nch == 0)
419  return FALSE;
420  ch_string = mne_channel_names_to_string_23(chs,nch);
421  mne_string_to_name_list_23(ch_string,list,nlist);
422  res = mne_proj_op_affect(op,list,nlist);
423  list.clear();
424  return res;
425 }
426 
427 //=============================================================================================================
428 
429 int MneProjOp::mne_proj_op_proj_vector(MneProjOp *op, float *vec, int nvec, int do_complement)
430 /*
431  * Apply projection operator to a vector (floats)
432  * Assume that all dimension checking etc. has been done before
433  */
434 {
435  static float *res = NULL;
436  int res_size = 0;
437  float *pvec;
438  float w;
439  int k,p;
440 
441  if (!op || op->nitems <= 0 || op->nvec <= 0)
442  return OK;
443 
444  if (op->nch != nvec) {
445  printf("Data vector size does not match projection operator");
446  return FAIL;
447  }
448 
449  if (op->nch > res_size) {
450  res = REALLOC_23(res,op->nch,float);
451  res_size = op->nch;
452  }
453 
454  for (k = 0; k < op->nch; k++)
455  res[k] = 0.0;
456 
457  for (p = 0; p < op->nvec; p++) {
458  pvec = op->proj_data[p];
459  w = mne_dot_vectors_23(pvec,vec,op->nch);
460  for (k = 0; k < op->nch; k++)
461  res[k] = res[k] + w*pvec[k];
462  }
463  if (do_complement) {
464  for (k = 0; k < op->nch; k++)
465  vec[k] = vec[k] - res[k];
466  }
467  else {
468  for (k = 0; k < op->nch; k++)
469  vec[k] = res[k];
470  }
471  return OK;
472 }
473 
474 //=============================================================================================================
475 
476 MneProjOp *MneProjOp::mne_read_proj_op_from_node(FiffStream::SPtr &stream, const FiffDirNode::SPtr &start)
477 /*
478  * Load all the linear projection data
479  */
480 {
481  MneProjOp* op = NULL;
482  QList<FiffDirNode::SPtr> proj;
483  FiffDirNode::SPtr start_node;
484  QList<FiffDirNode::SPtr> items;
485  FiffDirNode::SPtr node;
486  int k;
487  QString item_desc,desc_tag;
488  int global_nchan,item_nchan,nlist;
489  QStringList item_names;
490  int item_kind;
491  float **item_vectors = NULL;
492  int item_nvec;
493  int item_active;
494  MneNamedMatrix* item;
495  FiffTag::SPtr t_pTag;
496 
497  if (!stream) {
498  qCritical("File not open mne_read_proj_op_from_node");
499  goto bad;
500  }
501 
502  if (!start || start->isEmpty())
503  start_node = stream->dirtree();
504  else
505  start_node = start;
506 
507  op = new MneProjOp();
508  proj = start_node->dir_tree_find(FIFFB_PROJ);
509  if (proj.size() == 0 || proj[0]->isEmpty()) /* The caller must recognize an empty projection */
510  goto out;
511  /*
512  * Only the first projection block is recognized
513  */
514  items = proj[0]->dir_tree_find(FIFFB_PROJ_ITEM);
515  if (items.size() == 0 || items[0]->isEmpty()) /* The caller must recognize an empty projection */
516  goto out;
517  /*
518  * Get a common number of channels
519  */
520  node = proj[0];
521  if(!node->find_tag(stream, FIFF_NCHAN, t_pTag))
522  global_nchan = 0;
523  else {
524  global_nchan = *t_pTag->toInt();
525  // TAG_FREE(tag);
526  }
527  /*
528  * Proceess each item
529  */
530  for (k = 0; k < items.size(); k++) {
531  node = items[k];
532  /*
533  * Complicated procedure for getting the description
534  */
535  item_desc.clear();
536 
537  if (node->find_tag(stream, FIFF_NAME, t_pTag)) {
538  item_desc += t_pTag->toString();
539  }
540 
541  /*
542  * Take the first line of description if it exists
543  */
544  if (node->find_tag(stream, FIFF_DESCRIPTION, t_pTag)) {
545  desc_tag = t_pTag->toString();
546  int pos;
547  if((pos = desc_tag.indexOf("\n")) >= 0)
548  desc_tag.truncate(pos);
549  if (!item_desc.isEmpty())
550  item_desc += " ";
551  item_desc += desc_tag;
552  }
553  /*
554  * Possibility to override number of channels here
555  */
556  if (!node->find_tag(stream, FIFF_NCHAN, t_pTag)) {
557  item_nchan = global_nchan;
558  }
559  else {
560  item_nchan = *t_pTag->toInt();
561  }
562  if (item_nchan <= 0) {
563  qCritical("Number of channels incorrectly specified for one of the projection items.");
564  goto bad;
565  }
566  /*
567  * Take care of the channel names
568  */
569  if (!node->find_tag(stream, FIFF_PROJ_ITEM_CH_NAME_LIST, t_pTag))
570  goto bad;
571 
572  item_names = FiffStream::split_name_list(t_pTag->toString());
573 
574  if (item_names.size() != item_nchan) {
575  printf("Channel name list incorrectly specified for proj item # %d",k+1);
576  item_names.clear();
577  goto bad;
578  }
579  /*
580  * Kind of item
581  */
582  if (!node->find_tag(stream, FIFF_PROJ_ITEM_KIND, t_pTag))
583  goto bad;
584  item_kind = *t_pTag->toInt();
585  /*
586  * How many vectors
587  */
588  if (!node->find_tag(stream,FIFF_PROJ_ITEM_NVEC, t_pTag))
589  goto bad;
590  item_nvec = *t_pTag->toInt();
591  /*
592  * The projection data
593  */
594  if (!node->find_tag(stream,FIFF_PROJ_ITEM_VECTORS, t_pTag))
595  goto bad;
596 
597  MatrixXf tmp_item_vectors = t_pTag->toFloatMatrix().transpose();
598  item_vectors = ALLOC_CMATRIX_23(tmp_item_vectors.rows(),tmp_item_vectors.cols());
599  fromFloatEigenMatrix_23(tmp_item_vectors, item_vectors);
600 
601  /*
602  * Is this item active?
603  */
604  if (node->find_tag(stream, FIFF_MNE_PROJ_ITEM_ACTIVE, t_pTag)) {
605  item_active = *t_pTag->toInt();
606  }
607  else
608  item_active = FALSE;
609  /*
610  * Ready to add
611  */
612  QStringList emptyList;
613  item = MneNamedMatrix::build_named_matrix(item_nvec,item_nchan,emptyList,item_names,item_vectors);
614  mne_proj_op_add_item_act(op,item,item_kind,item_desc,item_active);
615  delete item;
616  op->items[op->nitems-1]->active_file = item_active;
617  }
618 
619 out :
620  return op;
621 
622 bad : {
623  if(op)
624  delete op;
625  return NULL;
626  }
627 }
628 
629 //=============================================================================================================
630 
631 MneProjOp *MneProjOp::mne_read_proj_op(const QString &name)
632 {
633  QFile file(name);
634  FiffStream::SPtr stream(new FiffStream(&file));
635 
636  if(!stream->open())
637  return NULL;
638 
639  MneProjOp* res = NULL;
640 
641  FiffDirNode::SPtr t_default;
642  res = mne_read_proj_op_from_node(stream,t_default);
643 
644  stream->close();
645 
646  return res;
647 }
648 
649 //=============================================================================================================
650 
651 void MneProjOp::mne_proj_op_report_data(FILE *out, const char *tag, MneProjOp *op, int list_data, char **exclude, int nexclude)
652 /*
653  * Output info about the projection operator
654  */
655 {
656  int j,k,p,q;
657  MneProjItem* it;
658  MneNamedMatrix* vecs;
659  int found;
660 
661  if (out == NULL)
662  return;
663  if (op == NULL)
664  return;
665  if (op->nitems <= 0) {
666  fprintf(out,"Empty operator\n");
667  return;
668  }
669 
670  for (k = 0; k < op->nitems; k++) {
671  it = op->items[k];
672  if (list_data && tag)
673  fprintf(out,"%s\n",tag);
674  if (tag)
675  fprintf(out,"%s",tag);
676  fprintf(out,"# %d : %s : %d vecs : %d chs %s %s\n",
677  k+1,it->desc.toUtf8().constData(),it->nvec,it->vecs->ncol,
678  it->has_meg ? "MEG" : "EEG",
679  it->active ? "active" : "idle");
680  if (list_data && tag)
681  fprintf(out,"%s\n",tag);
682  if (list_data) {
683  vecs = op->items[k]->vecs;
684 
685  for (q = 0; q < vecs->ncol; q++) {
686  fprintf(out,"%-10s",vecs->collist[q].toUtf8().constData());
687  fprintf(out,q < vecs->ncol-1 ? " " : "\n");
688  }
689  for (p = 0; p < vecs->nrow; p++)
690  for (q = 0; q < vecs->ncol; q++) {
691  for (j = 0, found = 0; j < nexclude; j++) {
692  if (QString::compare(exclude[j],vecs->collist[q]) == 0) {
693  found = 1;
694  break;
695  }
696  }
697  fprintf(out,"%10.5g ",found ? 0.0 : vecs->data[p][q]);
698  fprintf(out,q < vecs->ncol-1 ? " " : "\n");
699  }
700  if (list_data && tag)
701  fprintf(out,"%s\n",tag);
702  }
703  }
704  return;
705 }
706 
707 //=============================================================================================================
708 
709 void MneProjOp::mne_proj_op_report(FILE *out, const char *tag, MneProjOp *op)
710 {
711  mne_proj_op_report_data(out,tag,op, FALSE, NULL, 0);
712 }
MneNamedMatrix * vecs
One linear projection item.
Definition: mne_proj_op.h:83
static QStringList split_name_list(QString p_sNameList)
QSharedPointer< FiffDirNode > SPtr
Definition: fiff_dir_node.h:77
MNEProjItem class declaration.
One linear projection item.
Definition: mne_proj_item.h:77
#define FIFF_NCHAN
Definition: fiff_file.h:453
Fiff constants.
#define FIFF_MNE_PROJ_ITEM_ACTIVE
QSharedPointer< FiffTag > SPtr
Definition: fiff_tag.h:152
static void mne_free_proj_op_proj(MneProjOp *op)
#define FIFF_NAME
Definition: fiff_file.h:485
static MneNamedMatrix * build_named_matrix(int nrow, int ncol, const QStringList &rowlist, const QStringList &collist, float **data)
FIFF File I/O routines.
Definition: fiff_stream.h:104
Matrix specification with a channel list.
MNEProjOp class declaration.
QSharedPointer< FiffStream > SPtr
Definition: fiff_stream.h:107
FiffTag class declaration, which provides fiff tag I/O and processing methods.
#define FIFF_DESCRIPTION
Definition: fiff_file.h:486