MNE-CPP  0.1.9
A Framework for Electrophysiology
fwd_coil_set.cpp
Go to the documentation of this file.
1 //=============================================================================================================
37 //=============================================================================================================
38 // INCLUDES
39 //=============================================================================================================
40 
41 #include "fwd_coil_set.h"
42 #include "fwd_coil.h"
43 
44 #include <fiff/fiff_ch_info.h>
45 
46 //=============================================================================================================
47 // QT INCLUDES
48 //=============================================================================================================
49 
50 #include <QDebug>
51 
52 //=============================================================================================================
53 // USED NAMESPACES
54 //=============================================================================================================
55 
56 using namespace Eigen;
57 using namespace FIFFLIB;
58 using namespace FWDLIB;
59 
60 #define MAXWORD 1000
61 #define BIG 0.5
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 #define MALLOC_6(x,t) (t *)malloc((x)*sizeof(t))
80 #define REALLOC_6(x,y,t) (t *)((x == NULL) ? malloc((y)*sizeof(t)) : realloc((x),(y)*sizeof(t)))
81 #define FREE_6(x) if ((char *)(x) != NULL) free((char *)(x))
82 
83 #define FIFFV_COORD_UNKNOWN 0
84 
85 #define X_6 0
86 #define Y_6 1
87 #define Z_6 2
88 
89 #define VEC_DOT_6(x,y) ((x)[X_6]*(y)[X_6] + (x)[Y_6]*(y)[Y_6] + (x)[Z_6]*(y)[Z_6])
90 #define VEC_LEN_6(x) sqrt(VEC_DOT_6(x,x))
91 
92 #define VEC_COPY_6(to,from) {\
93  (to)[X_6] = (from)[X_6];\
94  (to)[Y_6] = (from)[Y_6];\
95  (to)[Z_6] = (from)[Z_6];\
96  }
97 
98 static void skip_comments(FILE *in)
99 
100 {
101  int c;
102 
103  while (1) {
104  c = fgetc(in);
105  if (c == '#') {
106  for (c = fgetc(in); c != EOF && c != '\n'; c = fgetc(in))
107  ;
108  }
109  else {
110  ungetc(c,in);
111  return;
112  }
113  }
114 }
115 
116 static int whitespace(int c)
117 
118 {
119  if (c == '\t' || c == '\n' || c == ' ')
120  return TRUE;
121  else
122  return FALSE;
123 }
124 
125 static int whitespace_quote(int c, int inquote)
126 
127 {
128  if (inquote)
129  return (c == '"');
130  else
131  return (c == '\t' || c == '\n' || c == ' ');
132 }
133 
134 static char *next_word(FILE *in)
135 
136 {
137  char *next = MALLOC_6(MAXWORD,char);
138  int c;
139  int p,k;
140  int inquote;
141 
142  skip_comments(in);
143 
144  inquote = FALSE;
145  for (k = 0, p = 0, c = fgetc(in); c != EOF && !whitespace_quote(c,inquote) ; c = fgetc(in), k++) {
146  if (k == 0 && c == '"')
147  inquote = TRUE;
148  else
149  next[p++] = c;
150  }
151  if (c == EOF && k == 0) {
152  FREE_6(next);
153  return NULL;
154  }
155  else
156  next[p] = '\0';
157  if (c != EOF) {
158  for (k = 0, c = fgetc(in); whitespace(c) ; c = fgetc(in), k++)
159  ;
160  if (c != EOF)
161  ungetc(c,in);
162  }
163 #ifdef DEBUG
164  if (next)
165  printf("<%s>\n",next);
166 #endif
167  return next;
168 }
169 
170 static int get_ival(FILE *in, int *ival)
171 
172 {
173  char *next = next_word(in);
174  if (next == NULL) {
175  qWarning("missing integer");
176  return FAIL;
177  }
178  else if (sscanf(next,"%d",ival) != 1) {
179  qWarning("bad integer : %s",next);
180  FREE_6(next);
181  return FAIL;
182  }
183  FREE_6(next);
184  return OK;
185 }
186 
187 static int get_fval(FILE *in, float *fval)
188 
189 {
190  char *next = next_word(in);
191  setlocale(LC_NUMERIC, "C");
192  if (next == NULL) {
193  qWarning("bad integer");
194  return FAIL;
195  }
196  else if (sscanf(next,"%g",fval) != 1) {
197  qWarning("bad floating point number : %s",next);
198  FREE_6(next);
199  return FAIL;
200  }
201  FREE_6(next);
202  return OK;
203 }
204 
205 static void normalize(float *rr)
206 /*
207  * Scale vector to unit length
208  */
209 {
210  float ll = VEC_LEN_6(rr);
211  int k;
212  if (ll > 0) {
213  for (k = 0; k < 3; k++)
214  rr[k] = rr[k]/ll;
215  }
216  return;
217 }
218 
219 static FwdCoil* fwd_add_coil_to_set(FwdCoilSet* set,
220  int type, int coil_class, int acc, int np, float size, float base, const QString& desc)
221 
222 {
223  FwdCoil* def;
224 
225  if (set == NULL) {
226  qWarning ("No coil definition set to augment.");
227  return NULL;
228  }
229  if (np <= 0) {
230  qWarning("Number of integration points should be positive (type = %d acc = %d)",type,acc);
231  return NULL;
232  }
233  if (! (acc == FWD_COIL_ACCURACY_POINT ||
234  acc == FWD_COIL_ACCURACY_NORMAL ||
235  acc == FWD_COIL_ACCURACY_ACCURATE) ) {
236  qWarning("Illegal accuracy (type = %d acc = %d)",type,acc);
237  return NULL;
238  }
239  if (! (coil_class == FWD_COILC_MAG ||
240  coil_class == FWD_COILC_AXIAL_GRAD ||
241  coil_class == FWD_COILC_PLANAR_GRAD ||
242  coil_class == FWD_COILC_AXIAL_GRAD2) ) {
243  qWarning("Illegal coil class (type = %d acc = %d class = %d)",type,acc,coil_class);
244  return NULL;
245  }
246 
247  set->coils = REALLOC_6(set->coils,set->ncoil+1,FwdCoil*);
248  def = set->coils[set->ncoil++] = new FwdCoil(np);
249 
250  def->type = type;
251  def->coil_class = coil_class;
252  def->accuracy = acc;
253  def->np = np;
254  def->base = size;
255  def->base = base;
256  if (!desc.isEmpty())
257  def->desc = desc;
258  return def;
259 }
260 
261 //=============================================================================================================
262 // DEFINE MEMBER METHODS
263 //=============================================================================================================
264 
265 FwdCoilSet::FwdCoilSet()
266 {
267  coils = NULL;
268  ncoil = 0;
269  coord_frame = FIFFV_COORD_UNKNOWN;
270  user_data = NULL;
271  user_data_free = NULL;
272 }
273 
274 //=============================================================================================================
275 
276 //FwdCoilSet::FwdCoilSet(const FwdCoilSet& p_FwdCoilSet)
277 //{
278 //}
279 
280 //=============================================================================================================
281 
282 FwdCoilSet::~FwdCoilSet()
283 {
284  for (int k = 0; k < ncoil; k++)
285  delete coils[k];
286  FREE_6(coils);
287 
288  this->fwd_free_coil_set_user_data();
289 }
290 
291 //=============================================================================================================
292 
293 FwdCoil *FwdCoilSet::create_meg_coil(const FiffChInfo& ch, int acc, const FiffCoordTransOld* t)
294 {
295  int k,p,c;
296  FwdCoil* def;
297  FwdCoil* res = NULL;
298 
299  if (ch.kind != FIFFV_MEG_CH && ch.kind != FIFFV_REF_MEG_CH) {
300  qWarning() << ch.ch_name << "is not a MEG channel. Cannot create a coil definition.";
301  goto bad;
302  }
303  /*
304  * Simple linear search from the coil definitions
305  */
306  for (k = 0, def = NULL; k < this->ncoil; k++) {
307  if ((this->coils[k]->type == (ch.chpos.coil_type & 0xFFFF)) &&
308  this->coils[k]->accuracy == acc) {
309  def = this->coils[k];
310  }
311  }
312  if (!def) {
313  printf("Desired coil definition not found (type = %d acc = %d)",ch.chpos.coil_type,acc);
314  goto bad;
315  }
316  /*
317  * Create the result
318  */
319  res = new FwdCoil(def->np);
320 
321  res->chname = ch.ch_name;
322  if (!def->desc.isEmpty())
323  res->desc = def->desc;
324  res->coil_class = def->coil_class;
325  res->accuracy = def->accuracy;
326  res->base = def->base;
327  res->size = def->size;
328  res->type = ch.chpos.coil_type;
329 
330  VEC_COPY_6(res->r0,ch.chpos.r0);
331  VEC_COPY_6(res->ex,ch.chpos.ex);
332  VEC_COPY_6(res->ey,ch.chpos.ey);
333  VEC_COPY_6(res->ez,ch.chpos.ez);
334  /*
335  * Apply a coordinate transformation if so desired
336  */
337  if (t) {
338  FiffCoordTransOld::fiff_coord_trans(res->r0,t,FIFFV_MOVE);
339  FiffCoordTransOld::fiff_coord_trans(res->ex,t,FIFFV_NO_MOVE);
340  FiffCoordTransOld::fiff_coord_trans(res->ey,t,FIFFV_NO_MOVE);
341  FiffCoordTransOld::fiff_coord_trans(res->ez,t,FIFFV_NO_MOVE);
342  res->coord_frame = t->to;
343  }
344  else
345  res->coord_frame = FIFFV_COORD_DEVICE;
346 
347  for (p = 0; p < res->np; p++) {
348  res->w[p] = def->w[p];
349  for (c = 0; c < 3; c++) {
350  res->rmag[p][c] = res->r0[c] + def->rmag[p][X_6]*res->ex[c] + def->rmag[p][Y_6]*res->ey[c] + def->rmag[p][Z_6]*res->ez[c];
351  res->cosmag[p][c] = def->cosmag[p][X_6]*res->ex[c] + def->cosmag[p][Y_6]*res->ey[c] + def->cosmag[p][Z_6]*res->ez[c];
352  }
353  }
354  return res;
355 
356 bad : {
357  return NULL;
358  }
359 }
360 
361 //=============================================================================================================
362 
363 FwdCoilSet *FwdCoilSet::create_meg_coils(const QList<FIFFLIB::FiffChInfo>& chs,
364  int nch,
365  int acc,
366  const FiffCoordTransOld* t)
367 {
368  FwdCoilSet* res = new FwdCoilSet();
369  FwdCoil* next;
370  int k;
371 
372  for (k = 0; k < nch; k++) {
373  if ((next = this->create_meg_coil(chs.at(k),acc,t)) == Q_NULLPTR)
374  goto bad;
375  res->coils = REALLOC_6(res->coils,res->ncoil+1,FwdCoil*);
376  res->coils[res->ncoil++] = next;
377  }
378  if (t)
379  res->coord_frame = t->to;
380  return res;
381 
382 bad : {
383  delete res;
384  return NULL;
385  }
386 }
387 
388 //=============================================================================================================
389 
390 FwdCoilSet *FwdCoilSet::create_eeg_els(const QList<FIFFLIB::FiffChInfo>& chs,
391  int nch,
392  const FiffCoordTransOld* t)
393 {
394  FwdCoilSet* res = new FwdCoilSet();
395  FwdCoil* next;
396  int k;
397 
398  for (k = 0; k < nch; k++) {
399  if ((next = FwdCoil::create_eeg_el(chs.at(k),t)) == Q_NULLPTR)
400  goto bad;
401  res->coils = REALLOC_6(res->coils,res->ncoil+1,FwdCoil*);
402  res->coils[res->ncoil++] = next;
403  }
404  if (t)
405  res->coord_frame = t->to;
406  return res;
407 
408 bad : {
409  delete res;
410  return NULL;
411  }
412 }
413 
414 //=============================================================================================================
415 
416 FwdCoilSet *FwdCoilSet::read_coil_defs(const QString &name)
417 /*
418  * Read a coil definition file
419  */
420 {
421  FILE *in = fopen(name.toUtf8().constData(),"r");
422  char *desc = NULL;
423  int type,coil_class,acc,np;
424  int p;
425  float size,base;
426  FwdCoilSet* res = NULL;
427  FwdCoil* def;
428 
429  if (in == NULL) {
430  qWarning() << "FwdCoilSet::read_coil_defs - File is NULL" << name;
431  goto bad;
432  }
433 
434  res = new FwdCoilSet();
435  while (1) {
436  /*
437  * Read basic info
438  */
439  if (get_ival(in,&coil_class) != OK)
440  break;
441  if (get_ival(in,&type) != OK)
442  goto bad;
443  if (get_ival(in,&acc) != OK)
444  goto bad;
445  if (get_ival(in,&np) != OK)
446  goto bad;
447  if (get_fval(in,&size) != OK)
448  goto bad;
449  if (get_fval(in,&base) != OK)
450  goto bad;
451  desc = next_word(in);
452  if (!desc)
453  goto bad;
454 
455  def = fwd_add_coil_to_set(res,type,coil_class,acc,np,size,base,desc);
456  if (!def)
457  goto bad;
458  FREE_6(desc); desc = NULL;
459 
460  for (p = 0; p < def->np; p++) {
461  /*
462  * Read and verify data for each integration point
463  */
464  if (get_fval(in,def->w+p) != OK)
465  goto bad;
466  if (get_fval(in,def->rmag[p]+X_6) != OK)
467  goto bad;
468  if (get_fval(in,def->rmag[p]+Y_6) != OK)
469  goto bad;
470  if (get_fval(in,def->rmag[p]+Z_6) != OK)
471  goto bad;
472  if (get_fval(in,def->cosmag[p]+X_6) != OK)
473  goto bad;
474  if (get_fval(in,def->cosmag[p]+Y_6) != OK)
475  goto bad;
476  if (get_fval(in,def->cosmag[p]+Z_6) != OK)
477  goto bad;
478 
479  if (VEC_LEN_6(def->rmag[p]) > BIG) {
480  qWarning("Unreasonable integration point: %f %f %f mm (coil type = %d acc = %d)", 1000*def->rmag[p][X_6],1000*def->rmag[p][Y_6],1000*def->rmag[p][Z_6], def->type,def->accuracy);
481  goto bad;
482  }
483  size = VEC_LEN_6(def->cosmag[p]);
484  if (size <= 0) {
485  qWarning("Unreasonable normal: %f %f %f (coil type = %d acc = %d)", def->cosmag[p][X_6],def->cosmag[p][Y_6],def->cosmag[p][Z_6], def->type,def->accuracy);
486  goto bad;
487  }
488  normalize(def->cosmag[p]);
489  }
490  }
491 
492  fclose(in);
493 
494  printf("%d coil definitions read\n",res->ncoil);
495  return res;
496 
497 bad : {
498  delete res;
499  FREE_6(desc);
500  return NULL;
501  }
502 }
503 
504 //=============================================================================================================
505 
506 FwdCoilSet* FwdCoilSet::dup_coil_set(const FiffCoordTransOld* t) const
507 {
508  FwdCoilSet* res;
509  FwdCoil* coil;
510 
511  if (t) {
512  if (this->coord_frame != t->from) {
513  qWarning("Coordinate frame of the transformation does not match the coil set in fwd_dup_coil_set");
514  return NULL;
515  }
516  }
517  res = new FwdCoilSet();
518  if (t)
519  res->coord_frame = t->to;
520  else
521  res->coord_frame = this->coord_frame;
522 
523  res->coils = MALLOC_6(this->ncoil,FwdCoil*);
524  res->ncoil = this->ncoil;
525 
526  for (int k = 0; k < this->ncoil; k++) {
527  coil = res->coils[k] = new FwdCoil(*(this->coils[k]));
528  /*
529  * Optional coordinate transformation
530  */
531  if (t) {
532  FiffCoordTransOld::fiff_coord_trans(coil->r0,t,FIFFV_MOVE);
533  FiffCoordTransOld::fiff_coord_trans(coil->ex,t,FIFFV_NO_MOVE);
534  FiffCoordTransOld::fiff_coord_trans(coil->ey,t,FIFFV_NO_MOVE);
535  FiffCoordTransOld::fiff_coord_trans(coil->ez,t,FIFFV_NO_MOVE);
536 
537  for (int p = 0; p < coil->np; p++) {
538  FiffCoordTransOld::fiff_coord_trans(coil->rmag[p],t,FIFFV_MOVE);
539  FiffCoordTransOld::fiff_coord_trans(coil->cosmag[p],t,FIFFV_NO_MOVE);
540  }
541  coil->coord_frame = t->to;
542  }
543  }
544  return res;
545 }
546 
547 //=============================================================================================================
548 
549 bool FwdCoilSet::is_planar_coil_type(int type) const
550 {
551  if (type == FIFFV_COIL_EEG)
552  return false;
553  for (int k = 0; k < this->ncoil; k++)
554  if (this->coils[k]->type == type)
555  return this->coils[k]->coil_class == FWD_COILC_PLANAR_GRAD;
556  return false;
557 }
558 
559 //=============================================================================================================
560 
561 bool FwdCoilSet::is_axial_coil_type(int type) const
562 {
563  if (type == FIFFV_COIL_EEG)
564  return false;
565  for (int k = 0; k < this->ncoil; k++)
566  if (this->coils[k]->type == type)
567  return (this->coils[k]->coil_class == FWD_COILC_MAG ||
568  this->coils[k]->coil_class == FWD_COILC_AXIAL_GRAD ||
569  this->coils[k]->coil_class == FWD_COILC_AXIAL_GRAD2);
570  return false;
571 }
572 
573 //=============================================================================================================
574 
575 bool FwdCoilSet::is_magnetometer_coil_type(int type) const
576 {
577  if (type == FIFFV_COIL_EEG)
578  return false;
579  for (int k = 0; k < this->ncoil; k++)
580  if (this->coils[k]->type == type)
581  return this->coils[k]->coil_class == FWD_COILC_MAG;
582  return false;
583 }
584 
585 //=============================================================================================================
586 
587 bool FwdCoilSet::is_eeg_electrode_type(int type) const
588 {
589  return type == FIFFV_COIL_EEG;
590 }
591 
Eigen::Vector3f ey
Definition: fiff_ch_pos.h:114
Channel info descriptor.
Definition: fiff_ch_info.h:74
float ** cosmag
Definition: fwd_coil.h:181
float ** rmag
Definition: fwd_coil.h:180
FwdCoilSet description.
Definition: fwd_coil_set.h:75
float ex[3]
Definition: fwd_coil.h:176
float r0[3]
Definition: fwd_coil.h:175
QString desc
Definition: fwd_coil.h:169
float ey[3]
Definition: fwd_coil.h:177
Eigen::Vector3f ez
Definition: fiff_ch_pos.h:115
Coordinate transformation descriptor.
#define FIFFV_COIL_EEG
Eigen::Vector3f ex
Definition: fiff_ch_pos.h:113
FwdCoil description.
Definition: fwd_coil.h:88
QString chname
Definition: fwd_coil.h:167
FiffChInfo class declaration.
FwdCoilSet class declaration.
float ez[3]
Definition: fwd_coil.h:178
#define FIFFV_REF_MEG_CH
fiff_int_t coil_type
Definition: fiff_ch_pos.h:111
Eigen::Vector3f r0
Definition: fiff_ch_pos.h:112
FwdCoil class declaration.