MNE-CPP  0.1.9
A Framework for Electrophysiology
mne_ctf_comp_data_set.cpp
Go to the documentation of this file.
1 //=============================================================================================================
37 //=============================================================================================================
38 // INCLUDES
39 //=============================================================================================================
40 
41 #include "mne_ctf_comp_data_set.h"
42 #include "mne_ctf_comp_data.h"
43 
44 #include "mne_types.h"
45 
47 
48 #include <fiff/fiff_types.h>
49 
50 #include <Eigen/Core>
51 
52 //=============================================================================================================
53 // QT INCLUDES
54 //=============================================================================================================
55 
56 #include <QFile>
57 
58 #ifndef TRUE
59 #define TRUE 1
60 #endif
61 
62 #ifndef FALSE
63 #define FALSE 0
64 #endif
65 
66 #ifndef FAIL
67 #define FAIL -1
68 #endif
69 
70 #ifndef OK
71 #define OK 0
72 #endif
73 
74 #define MALLOC_32(x,t) (t *)malloc((x)*sizeof(t))
75 #define REALLOC_32(x,y,t) (t *)((x == NULL) ? malloc((y)*sizeof(t)) : realloc((x),(y)*sizeof(t)))
76 
77 #define FREE_32(x) if ((char *)(x) != NULL) free((char *)(x))
78 
79 #define FREE_CMATRIX_32(m) mne_free_cmatrix_32((m))
80 #define ALLOC_CMATRIX_32(x,y) mne_cmatrix_32((x),(y))
81 
82 static void matrix_error_32(int kind, int nr, int nc)
83 
84 {
85  if (kind == 1)
86  printf("Failed to allocate memory pointers for a %d x %d matrix\n",nr,nc);
87  else if (kind == 2)
88  printf("Failed to allocate memory for a %d x %d matrix\n",nr,nc);
89  else
90  printf("Allocation error for a %d x %d matrix\n",nr,nc);
91  if (sizeof(void *) == 4) {
92  printf("This is probably because you seem to be using a computer with 32-bit architecture.\n");
93  printf("Please consider moving to a 64-bit platform.");
94  }
95  printf("Cannot continue. Sorry.\n");
96  exit(1);
97 }
98 
99 float **mne_cmatrix_32(int nr,int nc)
100 
101 {
102  int i;
103  float **m;
104  float *whole;
105 
106  m = MALLOC_32(nr,float *);
107  if (!m) matrix_error_32(1,nr,nc);
108  whole = MALLOC_32(nr*nc,float);
109  if (!whole) matrix_error_32(2,nr,nc);
110 
111  for(i=0;i<nr;i++)
112  m[i] = whole + i*nc;
113  return m;
114 }
115 
116 void mne_free_cmatrix_32 (float **m)
117 {
118  if (m) {
119  FREE_32(*m);
120  FREE_32(m);
121  }
122 }
123 
124 //=============================================================================================================
125 // USED NAMESPACES
126 //=============================================================================================================
127 
128 using namespace Eigen;
129 using namespace FIFFLIB;
130 using namespace MNELIB;
131 
132 //============================= mne_read_forward_solution.c =============================
133 
134 int mne_read_meg_comp_eeg_ch_info_32(const QString& name,
135  QList<FIFFLIB::FiffChInfo>& megp, /* MEG channels */
136  int *nmegp,
137  QList<FIFFLIB::FiffChInfo>& meg_compp,
138  int *nmeg_compp,
139  QList<FIFFLIB::FiffChInfo>& eegp, /* EEG channels */
140  int *neegp,
141  FiffCoordTransOld* *meg_head_t,
142  fiffId *idp) /* The measurement ID */
143 /*
144  * Read the channel information and split it into three arrays,
145  * one for MEG, one for MEG compensation channels, and one for EEG
146  */
147 {
148  QFile file(name);
149  FiffStream::SPtr stream(new FiffStream(&file));
150 
151  QList<FIFFLIB::FiffChInfo> chs;
152  int nchan = 0;
153  QList<FIFFLIB::FiffChInfo> meg;
154  int nmeg = 0;
155  QList<FIFFLIB::FiffChInfo> meg_comp;
156  int nmeg_comp = 0;
157  QList<FIFFLIB::FiffChInfo> eeg;
158  int neeg = 0;
159  fiffId id = NULL;
160  QList<FiffDirNode::SPtr> nodes;
161  FiffDirNode::SPtr info;
162  FiffTag::SPtr t_pTag;
163  FIFFLIB::FiffChInfo this_ch;
164  FiffCoordTransOld* t = NULL;
165  fiff_int_t kind, pos;
166  int j,k,to_find;
167 
168  if(!stream->open())
169  goto bad;
170 
171  nodes = stream->dirtree()->dir_tree_find(FIFFB_MNE_PARENT_MEAS_FILE);
172 
173  if (nodes.size() == 0) {
174  nodes = stream->dirtree()->dir_tree_find(FIFFB_MEAS_INFO);
175  if (nodes.size() == 0) {
176  qCritical ("Could not find the channel information.");
177  goto bad;
178  }
179  }
180  info = nodes[0];
181  to_find = 0;
182  for (k = 0; k < info->nent(); k++) {
183  kind = info->dir[k]->kind;
184  pos = info->dir[k]->pos;
185  switch (kind) {
186  case FIFF_NCHAN :
187  if (!stream->read_tag(t_pTag,pos))
188  goto bad;
189  nchan = *t_pTag->toInt();
190 
191  for (j = 0; j < nchan; j++) {
192  chs.append(FiffChInfo());
193  chs[j].scanNo = -1;
194  }
195  to_find = nchan;
196  break;
197 
198  case FIFF_PARENT_BLOCK_ID :
199  if(!stream->read_tag(t_pTag, pos))
200  goto bad;
201 // id = t_pTag->toFiffID();
202  *id = *(fiffId)t_pTag->data();
203  break;
204 
205  case FIFF_COORD_TRANS :
206  if(!stream->read_tag(t_pTag, pos))
207  goto bad;
208 // t = t_pTag->toCoordTrans();
209  t = FiffCoordTransOld::read_helper( t_pTag );
210  if (t->from != FIFFV_COORD_DEVICE || t->to != FIFFV_COORD_HEAD)
211  t = NULL;
212  break;
213 
214  case FIFF_CH_INFO : /* Information about one channel */
215  if(!stream->read_tag(t_pTag, pos))
216  goto bad;
217 
218  this_ch = t_pTag->toChInfo();
219  if (this_ch.scanNo <= 0 || this_ch.scanNo > nchan) {
220  printf ("FIFF_CH_INFO : scan # out of range %d (%d)!",this_ch.scanNo,nchan);
221  goto bad;
222  }
223  else
224  chs[this_ch.scanNo-1] = this_ch;
225  to_find--;
226  break;
227  }
228  }
229  if (to_find != 0) {
230  qCritical("Some of the channel information was missing.");
231  goto bad;
232  }
233  if (t == NULL && meg_head_t != NULL) {
234  /*
235  * Try again in a more general fashion
236  */
237  if ((t = FiffCoordTransOld::mne_read_meas_transform(name)) == NULL) {
238  qCritical("MEG -> head coordinate transformation not found.");
239  goto bad;
240  }
241  }
242  /*
243  * Sort out the channels
244  */
245  for (k = 0; k < nchan; k++) {
246  if (chs[k].kind == FIFFV_MEG_CH) {
247  meg.append(chs[k]);
248  nmeg++;
249  } else if (chs[k].kind == FIFFV_REF_MEG_CH) {
250  meg_comp.append(chs[k]);
251  nmeg_comp++;
252  } else if (chs[k].kind == FIFFV_EEG_CH) {
253  eeg.append(chs[k]);
254  neeg++;
255  }
256  }
257 // fiff_close(in);
258  stream->close();
259 
260  megp = meg;
261  if(nmegp) {
262  *nmegp = nmeg;
263  }
264 
265  meg_compp = meg_comp;
266  if(nmeg_compp) {
267  *nmeg_compp = nmeg_comp;
268  }
269 
270  eegp = eeg;
271  if(neegp) {
272  *neegp = neeg;
273  }
274 
275  if (idp == NULL) {
276  FREE_32(id);
277  }
278  else
279  *idp = id;
280  if (meg_head_t == NULL) {
281  FREE_32(t);
282  }
283  else
284  *meg_head_t = t;
285 
286  return FIFF_OK;
287 
288 bad : {
289 // fiff_close(in);
290  stream->close();
291  FREE_32(id);
292 // FREE_32(tag.data);
293  FREE_32(t);
294  return FIFF_FAIL;
295  }
296 }
297 
298 #define MNE_CTFV_COMP_UNKNOWN -1
299 #define MNE_CTFV_COMP_NONE 0
300 #define MNE_CTFV_COMP_G1BR 0x47314252
301 #define MNE_CTFV_COMP_G2BR 0x47324252
302 #define MNE_CTFV_COMP_G3BR 0x47334252
303 #define MNE_CTFV_COMP_G2OI 0x47324f49
304 #define MNE_CTFV_COMP_G3OI 0x47334f49
305 
306 static struct {
307  int grad_comp;
308  int ctf_comp;
309 } compMap[] = { { MNE_CTFV_NOGRAD, MNE_CTFV_COMP_NONE },
310 { MNE_CTFV_GRAD1, MNE_CTFV_COMP_G1BR },
311 { MNE_CTFV_GRAD2, MNE_CTFV_COMP_G2BR },
312 { MNE_CTFV_GRAD3, MNE_CTFV_COMP_G3BR },
313 { MNE_4DV_COMP1, MNE_4DV_COMP1 }, /* One-to-one mapping for 4D data */
314 { MNE_CTFV_COMP_UNKNOWN, MNE_CTFV_COMP_UNKNOWN }};
315 
316 int mne_unmap_ctf_comp_kind(int ctf_comp)
317 
318 {
319  int k;
320 
321  for (k = 0; compMap[k].grad_comp >= 0; k++)
322  if (ctf_comp == compMap[k].ctf_comp)
323  return compMap[k].grad_comp;
324  return ctf_comp;
325 }
326 
327 FiffSparseMatrix* mne_convert_to_sparse(float **dense, /* The dense matrix to be converted */
328  int nrow, /* Number of rows in the dense matrix */
329  int ncol, /* Number of columns in the dense matrix */
330  int stor_type, /* Either FIFFTS_MC_CCS or FIFFTS_MC_RCS */
331  float small) /* How small elements should be ignored? */
332 /*
333  * Create the compressed row or column storage sparse matrix representation
334  * including a vector containing the nonzero matrix element values,
335  * the row or column pointer vector and the appropriate index vector(s).
336  */
337 {
338  int j,k;
339  int nz;
340  int ptr;
341  FiffSparseMatrix* sparse = NULL;
342  int size;
343 
344  if (small < 0) { /* Automatic scaling */
345  float maxval = 0.0;
346  float val;
347 
348  for (j = 0; j < nrow; j++)
349  for (k = 0; k < ncol; k++) {
350  val = std::fabs(dense[j][k]);
351  if (val > maxval)
352  maxval = val;
353  }
354  if (maxval > 0)
355  small = maxval*std::fabs(small);
356  else
357  small = std::fabs(small);
358  }
359  for (j = 0, nz = 0; j < nrow; j++)
360  for (k = 0; k < ncol; k++) {
361  if (std::fabs(dense[j][k]) > small)
362  nz++;
363  }
364 
365  if (nz <= 0) {
366  printf("No nonzero elements found.");
367  return NULL;
368  }
369  if (stor_type == FIFFTS_MC_CCS) {
370  size = nz*(sizeof(fiff_float_t) + sizeof(fiff_int_t)) +
371  (ncol+1)*(sizeof(fiff_int_t));
372  }
373  else if (stor_type == FIFFTS_MC_RCS) {
374  size = nz*(sizeof(fiff_float_t) + sizeof(fiff_int_t)) +
375  (nrow+1)*(sizeof(fiff_int_t));
376  }
377  else {
378  printf("Unknown sparse matrix storage type: %d",stor_type);
379  return NULL;
380  }
381  sparse = new FiffSparseMatrix;
382  sparse->coding = stor_type;
383  sparse->m = nrow;
384  sparse->n = ncol;
385  sparse->nz = nz;
386  sparse->data = (float *)malloc(size);
387  sparse->inds = (int *)(sparse->data+nz);
388  sparse->ptrs = sparse->inds+nz;
389 
390  if (stor_type == FIFFTS_MC_RCS) {
391  for (j = 0, nz = 0; j < nrow; j++) {
392  ptr = -1;
393  for (k = 0; k < ncol; k++)
394  if (std::fabs(dense[j][k]) > small) {
395  sparse->data[nz] = dense[j][k];
396  if (ptr < 0)
397  ptr = nz;
398  sparse->inds[nz++] = k;
399  }
400  sparse->ptrs[j] = ptr;
401  }
402  sparse->ptrs[nrow] = nz;
403  for (j = nrow - 1; j >= 0; j--) /* Take care of the empty rows */
404  if (sparse->ptrs[j] < 0)
405  sparse->ptrs[j] = sparse->ptrs[j+1];
406  }
407  else if (stor_type == FIFFTS_MC_CCS) {
408  for (k = 0, nz = 0; k < ncol; k++) {
409  ptr = -1;
410  for (j = 0; j < nrow; j++)
411  if (std::fabs(dense[j][k]) > small) {
412  sparse->data[nz] = dense[j][k];
413  if (ptr < 0)
414  ptr = nz;
415  sparse->inds[nz++] = j;
416  }
417  sparse->ptrs[k] = ptr;
418  }
419  sparse->ptrs[ncol] = nz;
420  for (k = ncol-1; k >= 0; k--) /* Take care of the empty columns */
421  if (sparse->ptrs[k] < 0)
422  sparse->ptrs[k] = sparse->ptrs[k+1];
423  }
424  return sparse;
425 }
426 
427 int mne_sparse_mat_mult2_32(FiffSparseMatrix* mat, /* The sparse matrix */
428  float **mult, /* Matrix to be multiplied */
429  int ncol, /* How many columns in the above */
430  float **res) /* Result of the multiplication */
431 /*
432  * Multiply a dense matrix by a sparse matrix.
433  */
434 {
435  int i,j,k;
436  float val;
437 
438  if (mat->coding == FIFFTS_MC_RCS) {
439  for (i = 0; i < mat->m; i++) {
440  for (k = 0; k < ncol; k++) {
441  val = 0.0;
442  for (j = mat->ptrs[i]; j < mat->ptrs[i+1]; j++)
443  val += mat->data[j]*mult[mat->inds[j]][k];
444  res[i][k] = val;
445  }
446  }
447  }
448  else if (mat->coding == FIFFTS_MC_CCS) {
449  for (k = 0; k < ncol; k++) {
450  for (i = 0; i < mat->m; i++)
451  res[i][k] = 0.0;
452  for (i = 0; i < mat->n; i++)
453  for (j = mat->ptrs[i]; j < mat->ptrs[i+1]; j++)
454  res[mat->inds[j]][k] += mat->data[j]*mult[i][k];
455  }
456  }
457  else {
458  printf("mne_sparse_mat_mult2: unknown sparse matrix storage type: %d",mat->coding);
459  return -1;
460  }
461  return 0;
462 }
463 
464 float **mne_mat_mat_mult_32 (float **m1,float **m2,int d1,int d2,int d3)
465 /* Matrix multiplication
466  * result(d1 x d3) = m1(d1 x d2) * m2(d2 x d3) */
467 
468 {
469 #ifdef BLAS
470  float **result = ALLOC_CMATRIX_32(d1,d3);
471  char *transa = "N";
472  char *transb = "N";
473  float zero = 0.0;
474  float one = 1.0;
475  sgemm (transa,transb,&d3,&d1,&d2,
476  &one,m2[0],&d3,m1[0],&d2,&zero,result[0],&d3);
477  return (result);
478 #else
479  float **result = ALLOC_CMATRIX_32(d1,d3);
480  int j,k,p;
481  float sum;
482 
483  for (j = 0; j < d1; j++)
484  for (k = 0; k < d3; k++) {
485  sum = 0.0;
486  for (p = 0; p < d2; p++)
487  sum = sum + m1[j][p]*m2[p][k];
488  result[j][k] = sum;
489  }
490  return (result);
491 #endif
492 }
493 
494 int mne_sparse_vec_mult2_32(FiffSparseMatrix* mat, /* The sparse matrix */
495  float *vector, /* Vector to be multiplied */
496  float *res) /* Result of the multiplication */
497 /*
498  * Multiply a vector by a sparse matrix.
499  */
500 {
501  int i,j;
502 
503  if (mat->coding == FIFFTS_MC_RCS) {
504  for (i = 0; i < mat->m; i++) {
505  res[i] = 0.0;
506  for (j = mat->ptrs[i]; j < mat->ptrs[i+1]; j++)
507  res[i] += mat->data[j]*vector[mat->inds[j]];
508  }
509  return 0;
510  }
511  else if (mat->coding == FIFFTS_MC_CCS) {
512  for (i = 0; i < mat->m; i++)
513  res[i] = 0.0;
514  for (i = 0; i < mat->n; i++)
515  for (j = mat->ptrs[i]; j < mat->ptrs[i+1]; j++)
516  res[mat->inds[j]] += mat->data[j]*vector[i];
517  return 0;
518  }
519  else {
520  printf("mne_sparse_vec_mult2: unknown sparse matrix storage type: %d",mat->coding);
521  return -1;
522  }
523 }
524 
525 float mne_dot_vectors_32 (float *v1,
526  float *v2,
527  int nn)
528 
529 {
530 #ifdef BLAS
531  int one = 1;
532  float res = sdot(&nn,v1,&one,v2,&one);
533  return res;
534 #else
535  float res = 0.0;
536  int k;
537 
538  for (k = 0; k < nn; k++)
539  res = res + v1[k]*v2[k];
540  return res;
541 #endif
542 }
543 
544 void mne_mat_vec_mult2_32 (float **m,float *v,float *result, int d1,int d2)
545 /*
546  * Matrix multiplication
547  * result(d1) = m(d1 x d2) * v(d2)
548  */
549 
550 {
551  int j;
552 
553  for (j = 0; j < d1; j++)
554  result[j] = mne_dot_vectors_32 (m[j],v,d2);
555  return;
556 }
557 
558 //=============================================================================================================
559 // DEFINE MEMBER METHODS
560 //=============================================================================================================
561 
562 MneCTFCompDataSet::MneCTFCompDataSet()
563 :ncomp(0)
564 ,nch(0)
565 ,current(NULL)
566 ,undo(NULL)
567 {
568 }
569 
570 //=============================================================================================================
571 
573 {
574 // if (!set)
575 // return NULL;
576 
577  if (set.ncomp > 0) {
578 // this->comps = MALLOC_32(set.ncomp,mneCTFcompData);
579  this->ncomp = set.comps.size();
580  for (int k = 0; k < this->ncomp; k++)
581  if(set.comps[k])
582  this->comps.append(new MneCTFCompData(*set.comps[k]));
583  }
584 
585  if(set.current)
586  this->current = new MneCTFCompData(*set.current);
587 }
588 
589 //=============================================================================================================
590 
592 {
593 
594  for (int k = 0; k < comps.size(); k++)
595  if(comps[k])
596  delete comps[k];
597 
598  if(current)
599  delete current;
600 }
601 
602 //=============================================================================================================
603 
604 MneCTFCompDataSet *MneCTFCompDataSet::mne_read_ctf_comp_data(const QString &name)
605 /*
606  * Read all CTF compensation data from a given file
607  */
608 {
609  QFile file(name);
610  FiffStream::SPtr stream(new FiffStream(&file));
611 
612  MneCTFCompDataSet* set = NULL;
613  MneCTFCompData* one;
614  QList<FiffDirNode::SPtr> nodes;
615  QList<FiffDirNode::SPtr> comps;
616  int ncomp;
617  MneNamedMatrix* mat = NULL;
618  int kind,k;
619  FiffTag::SPtr t_pTag;
620  QList<FiffChInfo> chs;
621  int nch = 0;
622  int calibrated;
623  /*
624  * Read the channel information
625  */
626  {
627  QList<FiffChInfo> comp_chs, temp;
628  int ncompch = 0;
629 
630  if (mne_read_meg_comp_eeg_ch_info_32(name,
631  chs,
632  &nch,
633  comp_chs,
634  &ncompch,
635  temp,
636  NULL,
637  NULL,
638  NULL) == FAIL)
639  goto bad;
640  if (ncompch > 0) {
641  for (k = 0; k < ncompch; k++)
642  chs.append(comp_chs[k]);
643  nch = nch + ncompch;
644  }
645  }
646  /*
647  * Read the rest of the stuff
648  */
649  if(!stream->open())
650  goto bad;
651  set = new MneCTFCompDataSet();
652  /*
653  * Locate the compensation data sets
654  */
655  nodes = stream->dirtree()->dir_tree_find(FIFFB_MNE_CTF_COMP);
656  if (nodes.size() == 0)
657  goto good; /* Nothing more to do */
658  comps = nodes[0]->dir_tree_find(FIFFB_MNE_CTF_COMP_DATA);
659  if (comps.size() == 0)
660  goto good;
661  ncomp = comps.size();
662  /*
663  * Set the channel info
664  */
665  set->chs = chs;
666  set->nch = nch;
667  /*
668  * Read each data set
669  */
670  for (k = 0; k < ncomp; k++) {
672  if (!mat)
673  goto bad;
674  comps[k]->find_tag(stream, FIFF_MNE_CTF_COMP_KIND, t_pTag);
675  if (t_pTag) {
676  kind = *t_pTag->toInt();
677  }
678  else
679  goto bad;
680  comps[k]->find_tag(stream, FIFF_MNE_CTF_COMP_CALIBRATED, t_pTag);
681  if (t_pTag) {
682  calibrated = *t_pTag->toInt();
683  }
684  else
685  calibrated = FALSE;
686  /*
687  * Add these data to the set
688  */
689  one = new MneCTFCompData;
690  one->data = mat; mat = NULL;
691  one->kind = kind;
692  one->mne_kind = mne_unmap_ctf_comp_kind(one->kind);
693  one->calibrated = calibrated;
694 
695  if (MneCTFCompData::mne_calibrate_ctf_comp(one,set->chs,set->nch,TRUE) == FAIL) {
696  printf("Warning: Compensation data for '%s' omitted\n", mne_explain_ctf_comp(one->kind));//,err_get_error(),mne_explain_ctf_comp(one->kind));
697  if(one)
698  delete one;
699  }
700  else {
701  // set->comps = REALLOC_9(set->comps,set->ncomp+1,mneCTFcompData);
702  // set->comps[set->ncomp++] = one;
703  set->comps.append(one);
704  set->ncomp++;
705  }
706  }
707 #ifdef DEBUG
708  printf("%d CTF compensation data sets read from %s\n",set->ncomp,name);
709 #endif
710  goto good;
711 
712 bad : {
713  if(mat)
714  delete mat;
715  stream->close();
716  if(set)
717  delete set;
718  return NULL;
719  }
720 
721 good : {
722  stream->close();
723  return set;
724  }
725 }
726 
727 //=============================================================================================================
728 
729 int MneCTFCompDataSet::mne_make_ctf_comp(MneCTFCompDataSet* set,
730  const QList<FiffChInfo>& chs,
731  int nch,
732  QList<FiffChInfo> compchs,
733  int ncomp) /* How many of these */
734 /*
735  * Make compensation data to apply to a set of channels to yield (or uncompensated) compensated data
736  */
737 {
738  int *comps = NULL;
739  int need_comp;
740  int first_comp;
741  MneCTFCompData* this_comp;
742  int *comp_sel = NULL;
743  QStringList names;
744  QString name;
745  int j,k,p;
746 
747  FiffSparseMatrix* presel = NULL;
748  FiffSparseMatrix* postsel = NULL;
749  MneNamedMatrix* data = NULL;
750 
751  QStringList emptyList;
752 
753  if (compchs.isEmpty()) {
754  compchs = chs;
755  ncomp = nch;
756  }
757  printf("Setting up compensation data...\n");
758  if (nch == 0)
759  return OK;
760  if (set) {
761  if(set->current)
762  delete set->current;
763  set->current = NULL;
764  }
765  comps = MALLOC_32(nch,int);
766  for (k = 0, need_comp = 0, first_comp = MNE_CTFV_COMP_NONE; k < nch; k++) {
767  if (chs[k].kind == FIFFV_MEG_CH) {
768  comps[k] = chs[k].chpos.coil_type >> 16;
769  if (comps[k] != MNE_CTFV_COMP_NONE) {
770  if (first_comp == MNE_CTFV_COMP_NONE)
771  first_comp = comps[k];
772  else {
773  if (comps[k] != first_comp) {
774  printf("We do not support nonuniform compensation yet.");
775  goto bad;
776  }
777  }
778  need_comp++;
779  }
780  }
781  else
782  comps[k] = MNE_CTFV_COMP_NONE;
783  }
784  if (need_comp == 0) {
785  printf("\tNo compensation set. Nothing more to do.\n");
786  FREE_32(comps);
787  return OK;
788  }
789  printf("\t%d out of %d channels have the compensation set.\n",need_comp,nch);
790  if (!set) {
791  printf("No compensation data available for the required compensation.");
792  return FAIL;
793  }
794  /*
795  * Find the desired compensation data matrix
796  */
797  for (k = 0, this_comp = NULL; k < set->ncomp; k++) {
798  if (set->comps[k]->mne_kind == first_comp) {
799  this_comp = set->comps[k];
800  break;
801  }
802  }
803  if (!this_comp) {
804  printf("Did not find the desired compensation data : %s",
805  mne_explain_ctf_comp(mne_map_ctf_comp_kind(first_comp)));
806  goto bad;
807  }
808  printf("\tDesired compensation data (%s) found.\n",mne_explain_ctf_comp(mne_map_ctf_comp_kind(first_comp)));
809  /*
810  * Find the compensation channels
811  */
812  comp_sel = MALLOC_32(this_comp->data->ncol,int);
813  for (k = 0; k < this_comp->data->ncol; k++) {
814  comp_sel[k] = -1;
815  name = this_comp->data->collist[k];
816  for (p = 0; p < ncomp; p++)
817  if (QString::compare(name,compchs[p].ch_name) == 0) {
818  comp_sel[k] = p;
819  break;
820  }
821  if (comp_sel[k] < 0) {
822  printf("Compensation channel %s not found",name.toUtf8().constData());
823  goto bad;
824  }
825  }
826  printf("\tAll compensation channels found.\n");
827  /*
828  * Create the preselector
829  */
830  {
831  float **sel = ALLOC_CMATRIX_32(this_comp->data->ncol,ncomp);
832  for (j = 0; j < this_comp->data->ncol; j++) {
833  for (k = 0; k < ncomp; k++)
834  sel[j][k] = 0.0;
835  sel[j][comp_sel[j]] = 1.0;
836  }
837  if ((presel = mne_convert_to_sparse(sel,this_comp->data->ncol,ncomp,FIFFTS_MC_RCS,1e-30)) == NULL) {
838  FREE_CMATRIX_32(sel);
839  goto bad;
840  }
841  FREE_CMATRIX_32(sel);
842  printf("\tPreselector created.\n");
843  }
844  /*
845  * Pick the desired channels
846  */
847  for (k = 0; k < nch; k++) {
848  if (comps[k] != MNE_CTFV_COMP_NONE)
849  names.append(chs[k].ch_name);
850  }
851 
852  if ((data = this_comp->data->pick_from_named_matrix(names,need_comp,emptyList,0)) == NULL)
853  goto bad;
854  printf("\tCompensation data matrix created.\n");
855  /*
856  * Create the postselector
857  */
858  {
859  float **sel = ALLOC_CMATRIX_32(nch,data->nrow);
860  for (j = 0, p = 0; j < nch; j++) {
861  for (k = 0; k < data->nrow; k++)
862  sel[j][k] = 0.0;
863  if (comps[j] != MNE_CTFV_COMP_NONE)
864  sel[j][p++] = 1.0;
865  }
866  if ((postsel = mne_convert_to_sparse(sel,nch,data->nrow,FIFFTS_MC_RCS,1e-30)) == NULL) {
867  FREE_CMATRIX_32(sel);
868  goto bad;
869  }
870  FREE_CMATRIX_32(sel);
871  printf("\tPostselector created.\n");
872  }
873  set->current = new MneCTFCompData();
874  set->current->kind = this_comp->kind;
875  set->current->mne_kind = this_comp->mne_kind;
876  set->current->data = data;
877  set->current->presel = presel;
878  set->current->postsel = postsel;
879 
880  printf("\tCompensation set up.\n");
881 
882  names.clear();
883  FREE_32(comps);
884  FREE_32(comp_sel);
885 
886  return OK;
887 
888 bad : {
889  if(presel)
890  delete presel;
891  if(postsel)
892  delete postsel;
893  if(data)
894  delete data;
895  names.clear();
896  FREE_32(comps);
897  FREE_32(comp_sel);
898  return FAIL;
899  }
900 }
901 
902 //=============================================================================================================
903 
904 int MneCTFCompDataSet::mne_set_ctf_comp(QList<FIFFLIB::FiffChInfo>& chs,
905  int nch,
906  int comp)
907 /*
908  * Set the compensation bits to the desired value
909  */
910 {
911  int k;
912  int nset;
913  for (k = 0, nset = 0; k < nch; k++) {
914  if (chs[k].kind == FIFFV_MEG_CH) {
915  chs[k].chpos.coil_type = (chs[k].chpos.coil_type & 0xFFFF) | (comp << 16);
916  nset++;
917  }
918  }
919  printf("A new compensation value (%s) was assigned to %d MEG channels.\n",
920  mne_explain_ctf_comp(mne_map_ctf_comp_kind(comp)),nset);
921  return nset;
922 }
923 
924 //=============================================================================================================
925 
926 int MneCTFCompDataSet::mne_apply_ctf_comp(MneCTFCompDataSet *set, int do_it, float *data, int ndata, float *compdata, int ncompdata)
927 /*
928  * Apply compensation or revert to uncompensated data
929  */
930 {
931  MneCTFCompData* this_comp;
932  float *presel,*comp;
933  int k;
934 
935  if (compdata == NULL) {
936  compdata = data;
937  ncompdata = ndata;
938  }
939  if (!set || !set->current)
940  return OK;
941  this_comp = set->current;
942  /*
943  * Dimension checks
944  */
945  if (this_comp->presel) {
946  if (this_comp->presel->n != ncompdata) {
947  printf("Compensation data dimension mismatch. Expected %d, got %d channels.",
948  this_comp->presel->n,ncompdata);
949  return FAIL;
950  }
951  }
952  else if (this_comp->data->ncol != ncompdata) {
953  printf("Compensation data dimension mismatch. Expected %d, got %d channels.",
954  this_comp->data->ncol,ncompdata);
955  return FAIL;
956  }
957  if (this_comp->postsel) {
958  if (this_comp->postsel->m != ndata) {
959  printf("Data dimension mismatch. Expected %d, got %d channels.",
960  this_comp->postsel->m,ndata);
961  return FAIL;
962  }
963  }
964  else if (this_comp->data->nrow != ndata) {
965  printf("Data dimension mismatch. Expected %d, got %d channels.",
966  this_comp->data->nrow,ndata);
967  return FAIL;
968  }
969  /*
970  * Preselection is optional
971  */
972  if (this_comp->presel) {
973  if (!this_comp->presel_data)
974  this_comp->presel_data = MALLOC_32(this_comp->presel->m,float);
975  if (mne_sparse_vec_mult2_32(this_comp->presel,compdata,this_comp->presel_data) != OK)
976  return FAIL;
977  presel = this_comp->presel_data;
978  }
979  else
980  presel = compdata;
981  /*
982  * This always happens
983  */
984  if (!this_comp->comp_data)
985  this_comp->comp_data = MALLOC_32(this_comp->data->nrow,float);
986  mne_mat_vec_mult2_32(this_comp->data->data,presel,this_comp->comp_data,this_comp->data->nrow,this_comp->data->ncol);
987  /*
988  * Optional postselection
989  */
990  if (!this_comp->postsel)
991  comp = this_comp->comp_data;
992  else {
993  if (!this_comp->postsel_data) {
994  this_comp->postsel_data = MALLOC_32(this_comp->postsel->m,float);
995  }
996  if (mne_sparse_vec_mult2_32(this_comp->postsel,this_comp->comp_data,this_comp->postsel_data) != OK)
997  return FAIL;
998  comp = this_comp->postsel_data;
999  }
1000  /*
1001  * Compensate or revert compensation?
1002  */
1003  if (do_it) {
1004  for (k = 0; k < ndata; k++)
1005  data[k] = data[k] - comp[k];
1006  }
1007  else {
1008  for (k = 0; k < ndata; k++)
1009  data[k] = data[k] + comp[k];
1010  }
1011  return OK;
1012 }
1013 
1014 //=============================================================================================================
1015 
1016 int MneCTFCompDataSet::mne_apply_ctf_comp_t(MneCTFCompDataSet *set, int do_it, float **data, int ndata, int ns) /* Number of samples */
1017 /*
1018  * Apply compensation or revert to uncompensated data
1019  */
1020 {
1021  MneCTFCompData* this_comp;
1022  float **presel,**comp;
1023  float **compdata = data;
1024  int ncompdata = ndata;
1025  int k,p;
1026 
1027  if (!set || !set->current)
1028  return OK;
1029  this_comp = set->current;
1030  /*
1031  * Dimension checks
1032  */
1033  if (this_comp->presel) {
1034  if (this_comp->presel->n != ncompdata) {
1035  printf("Compensation data dimension mismatch. Expected %d, got %d channels.",
1036  this_comp->presel->n,ncompdata);
1037  return FAIL;
1038  }
1039  }
1040  else if (this_comp->data->ncol != ncompdata) {
1041  printf("Compensation data dimension mismatch. Expected %d, got %d channels.",
1042  this_comp->data->ncol,ncompdata);
1043  return FAIL;
1044  }
1045  if (this_comp->postsel) {
1046  if (this_comp->postsel->m != ndata) {
1047  printf("Data dimension mismatch. Expected %d, got %d channels.",
1048  this_comp->postsel->m,ndata);
1049  return FAIL;
1050  }
1051  }
1052  else if (this_comp->data->nrow != ndata) {
1053  printf("Data dimension mismatch. Expected %d, got %d channels.",
1054  this_comp->data->nrow,ndata);
1055  return FAIL;
1056  }
1057  /*
1058  * Preselection is optional
1059  */
1060  if (this_comp->presel) {
1061  presel = ALLOC_CMATRIX_32(this_comp->presel->m,ns);
1062  if (mne_sparse_mat_mult2_32(this_comp->presel,compdata,ns,presel) != OK) {
1063  FREE_CMATRIX_32(presel);
1064  return FAIL;
1065  }
1066  }
1067  else
1068  presel = data;
1069  /*
1070  * This always happens
1071  */
1072  comp = mne_mat_mat_mult_32(this_comp->data->data,presel,this_comp->data->nrow,this_comp->data->ncol,ns);
1073  if (this_comp->presel)
1074  FREE_CMATRIX_32(presel);
1075  /*
1076  * Optional postselection
1077  */
1078  if (this_comp->postsel) {
1079  float **postsel = ALLOC_CMATRIX_32(this_comp->postsel->m,ns);
1080  if (mne_sparse_mat_mult2_32(this_comp->postsel,comp,ns,postsel) != OK) {
1081  FREE_CMATRIX_32(postsel);
1082  return FAIL;
1083  }
1084  FREE_CMATRIX_32(comp);
1085  comp = postsel;
1086  }
1087  /*
1088  * Compensate or revert compensation?
1089  */
1090  if (do_it) {
1091  for (k = 0; k < ndata; k++)
1092  for (p = 0; p < ns; p++)
1093  data[k][p] = data[k][p] - comp[k][p];
1094  }
1095  else {
1096  for (k = 0; k < ndata; k++)
1097  for (p = 0; p < ns; p++)
1098  data[k][p] = data[k][p] + comp[k][p];
1099  }
1100  FREE_CMATRIX_32(comp);
1101  return OK;
1102 }
1103 
1104 //=============================================================================================================
1105 
1106 int MneCTFCompDataSet::mne_get_ctf_comp(const QList<FIFFLIB::FiffChInfo> &chs, int nch)
1107 {
1108  int res = MNE_CTFV_NOGRAD;
1109  int first_comp,comp;
1110  int k;
1111 
1112  for (k = 0, first_comp = -1; k < nch; k++) {
1113  if (chs[k].kind == FIFFV_MEG_CH) {
1114  comp = chs[k].chpos.coil_type >> 16;
1115  if (first_comp < 0)
1116  first_comp = comp;
1117  else if (first_comp != comp) {
1118  printf("Non uniform compensation not supported.");
1119  return FAIL;
1120  }
1121  }
1122  }
1123  if (first_comp >= 0)
1124  res = first_comp;
1125  return res;
1126 }
1127 
1128 //=============================================================================================================
1129 
1130 int MneCTFCompDataSet::mne_map_ctf_comp_kind(int grad)
1131 /*
1132  * Simple mapping
1133  */
1134 {
1135  int k;
1136 
1137  for (k = 0; compMap[k].grad_comp >= 0; k++)
1138  if (grad == compMap[k].grad_comp)
1139  return compMap[k].ctf_comp;
1140  return grad;
1141 }
1142 
1143 //=============================================================================================================
1144 
1145 const char *MneCTFCompDataSet::mne_explain_ctf_comp(int kind)
1146 {
1147  static struct {
1148  int kind;
1149  const char *expl;
1150  } explain[] = { { MNE_CTFV_COMP_NONE, "uncompensated" },
1151  { MNE_CTFV_COMP_G1BR, "first order gradiometer" },
1152  { MNE_CTFV_COMP_G2BR, "second order gradiometer" },
1153  { MNE_CTFV_COMP_G3BR, "third order gradiometer" },
1154  { MNE_4DV_COMP1, "4D comp 1" },
1155  { MNE_CTFV_COMP_UNKNOWN, "unknown" } };
1156  int k;
1157 
1158  for (k = 0; explain[k].kind != MNE_CTFV_COMP_UNKNOWN; k++)
1159  if (explain[k].kind == kind)
1160  return explain[k].expl;
1161  return explain[k].expl;
1162 }
1163 
1164 //=============================================================================================================
1165 
1167  int compensate_to,
1168  QList<FiffChInfo>& chs,
1169  int nchan,
1170  QList<FiffChInfo> comp_chs,
1171  int ncomp_chan) /* How many */
1172 /*
1173  * Make data which has the third-order gradient compensation applied
1174  */
1175 {
1176  int k;
1177  int have_comp_chs;
1178  int comp_was = MNE_CTFV_COMP_UNKNOWN;
1179 
1180  if (!set) {
1181  if (compensate_to == MNE_CTFV_NOGRAD)
1182  return OK;
1183  else {
1184  printf("Cannot do compensation because compensation data are missing");
1185  return FAIL;
1186  }
1187  }
1188  if (comp_chs.isEmpty()) {
1189  comp_chs = chs;
1190  ncomp_chan = nchan;
1191  }
1192  if (set) {
1193  if(set->undo)
1194  delete set->undo;
1195  set->undo = NULL;
1196  if(set->current)
1197  delete set->current;
1198  set->current = NULL;
1199  }
1200  for (k = 0, have_comp_chs = 0; k < ncomp_chan; k++)
1201  if (comp_chs[k].kind == FIFFV_REF_MEG_CH)
1202  have_comp_chs++;
1203  if (have_comp_chs == 0 && compensate_to != MNE_CTFV_NOGRAD) {
1204  printf("No compensation channels in these data.");
1205  return FAIL;
1206  }
1207  /*
1208  * Update the 'current' field in 'set' to reflect the compensation possibly present in the data now
1209  */
1210  if (mne_make_ctf_comp(set,chs,nchan,comp_chs,ncomp_chan) == FAIL)
1211  goto bad;
1212  /*
1213  * Are we there already?
1214  */
1215  if (set->current && set->current->mne_kind == compensate_to) {
1216  printf("No further compensation necessary (comp = %s)\n",mne_explain_ctf_comp(set->current->kind));
1217  if(set->current)
1218  delete set->current;
1219  set->current = NULL;
1220  return OK;
1221  }
1222  set->undo = set->current;
1223  set->current = NULL;
1224  if (compensate_to == MNE_CTFV_NOGRAD) {
1225  printf("No compensation was requested.\n");
1226  mne_set_ctf_comp(chs,nchan,compensate_to);
1227  return OK;
1228  }
1229  if (mne_set_ctf_comp(chs,nchan,compensate_to) > 0) {
1230  if (set->undo)
1231  comp_was = set->undo->mne_kind;
1232  else
1233  comp_was = MNE_CTFV_NOGRAD;
1234  if (mne_make_ctf_comp(set,chs,nchan,comp_chs,ncomp_chan) == FAIL)
1235  goto bad;
1236  printf("Compensation set up as requested (%s -> %s).\n",
1237  mne_explain_ctf_comp(mne_map_ctf_comp_kind(comp_was)),
1238  mne_explain_ctf_comp(set->current->kind));
1239  }
1240  return OK;
1241 
1242 bad : {
1243  if (comp_was != MNE_CTFV_COMP_UNKNOWN)
1244  mne_set_ctf_comp(chs,nchan,comp_was);
1245  return FAIL;
1246  }
1247 }
Old fiff_type declarations - replace them.
MneNamedMatrix * pick_from_named_matrix(const QStringList &pickrowlist, int picknrow, const QStringList &pickcollist, int pickncol) const
Channel info descriptor.
Definition: fiff_ch_info.h:74
ToDo Old implementation use new fiff_id.h instead.
Definition: fiff_types.h:218
QSharedPointer< FiffDirNode > SPtr
Definition: fiff_dir_node.h:77
static int mne_ctf_set_compensation(MneCTFCompDataSet *set, int compensate_to, QList< FIFFLIB::FiffChInfo > &chs, int nchan, QList< FIFFLIB::FiffChInfo > comp_chs, int ncomp_chan)
FIFFLIB::fiff_int_t coding
#define FIFF_MNE_CTF_COMP_CALIBRATED
#define FIFF_CH_INFO
Definition: fiff_file.h:456
#define FIFF_NCHAN
Definition: fiff_file.h:453
One MNE CTF compensation description.
#define FIFF_MNE_CTF_COMP_KIND
QSharedPointer< FiffTag > SPtr
Definition: fiff_tag.h:152
One MNE CTF Compensation Data Set description.
FIFF File I/O routines.
Definition: fiff_stream.h:104
Coordinate transformation descriptor.
Matrix specification with a channel list.
MneCTFCompDataSet class declaration.
#define FIFF_COORD_TRANS
Definition: fiff_file.h:475
FIFFLIB::fiff_int_t * ptrs
FIFFLIB::fiff_int_t * inds
static MneNamedMatrix * read_named_matrix(QSharedPointer< FIFFLIB::FiffStream > &stream, const QSharedPointer< FIFFLIB::FiffDirNode > &node, int kind)
#define FIFF_MNE_CTF_COMP_DATA
QSharedPointer< FiffStream > SPtr
Definition: fiff_stream.h:107
FIFFLIB::fiff_float_t * data
#define FIFFV_REF_MEG_CH
Data associated with MNE computations for each mneMeasDataSet.
MneCTFCompData class declaration.
FiffCoordTransOld class declaration.