MNE-CPP  0.1.9
A Framework for Electrophysiology
fwd_comp_data.cpp
Go to the documentation of this file.
1 //=============================================================================================================
37 //=============================================================================================================
38 // INCLUDES
39 //=============================================================================================================
40 
41 #include "fwd_comp_data.h"
42 
44 #include <fiff/fiff_types.h>
45 
46 #include <iostream>
47 
48 #ifndef TRUE
49 #define TRUE 1
50 #endif
51 
52 #ifndef FALSE
53 #define FALSE 0
54 #endif
55 
56 #ifndef FAIL
57 #define FAIL -1
58 #endif
59 
60 #ifndef OK
61 #define OK 0
62 #endif
63 
64 #define X_60 0
65 #define Y_60 1
66 #define Z_60 2
67 
68 #define MALLOC_60(x,t) (t *)malloc((x)*sizeof(t))
69 
70 #define ALLOC_CMATRIX_60(x,y) mne_cmatrix_60((x),(y))
71 
72 #define FREE_60(x) if ((char *)(x) != NULL) free((char *)(x))
73 
74 #define FREE_CMATRIX_60(m) mne_free_cmatrix_60((m))
75 
76 void mne_free_cmatrix_60 (float **m)
77 {
78  if (m) {
79  FREE_60(*m);
80  FREE_60(m);
81  }
82 }
83 
84 static void matrix_error_60 (int kind, int nr, int nc)
85 {
86  if (kind == 1)
87  printf("Failed to allocate memory pointers for a %d x %d matrix\n",nr,nc);
88  else if (kind == 2)
89  printf("Failed to allocate memory for a %d x %d matrix\n",nr,nc);
90  else
91  printf("Allocation error for a %d x %d matrix\n",nr,nc);
92  if (sizeof(void *) == 4) {
93  printf("This is probably because you seem to be using a computer with 32-bit architecture.\n");
94  printf("Please consider moving to a 64-bit platform.");
95  }
96  printf("Cannot continue. Sorry.\n");
97  exit(1);
98 }
99 
100 float **mne_cmatrix_60 (int nr,int nc)
101 {
102  int i;
103  float **m;
104  float *whole;
105 
106  m = MALLOC_60(nr,float *);
107  if (!m) matrix_error_60(1,nr,nc);
108  whole = MALLOC_60(nr*nc,float);
109  if (!whole) matrix_error_60(2,nr,nc);
110 
111  for(i=0;i<nr;i++)
112  m[i] = whole + i*nc;
113  return m;
114 }
115 
116 //=============================================================================================================
117 // USED NAMESPACES
118 //=============================================================================================================
119 
120 using namespace Eigen;
121 using namespace FIFFLIB;
122 using namespace MNELIB;
123 using namespace FWDLIB;
124 
125 //=============================================================================================================
126 // DEFINE MEMBER METHODS
127 //=============================================================================================================
128 
129 FwdCompData::FwdCompData()
130 :comp_coils (NULL)
131 ,field (NULL)
132 ,vec_field (NULL)
133 ,field_grad (NULL)
134 ,client (NULL)
135 ,client_free(NULL)
136 ,set (NULL)
137 ,work (NULL)
138 ,vec_work (NULL)
139 {
140 }
141 
142 //=============================================================================================================
143 
145 {
146 // fwd_free_comp_data((void *)this);
147  if(this->comp_coils)
148  delete this->comp_coils;
149  if(this->set)
150  delete this->set;
151  FREE_60(this->work);
152  FREE_CMATRIX_60(this->vec_work);
153 
154  if (this->client_free && this->client)
155  this->client_free(this->client);
156 }
157 
158 //=============================================================================================================
159 
160 int FwdCompData::fwd_comp_field(float *rd, float *Q, FwdCoilSet *coils, float *res, void *client)
161 /*
162  * Calculate the compensated field (one dipole component)
163  */
164 {
165  FwdCompData* comp = (FwdCompData*)client;
166 
167  if (!comp->field) {
168  printf("Field computation function is missing in fwd_comp_field_vec");
169  return FAIL;
170  }
171  /*
172  * First compute the field in the primary set of coils
173  */
174  if (comp->field(rd,Q,coils,res,comp->client) == FAIL)
175  return FAIL;
176  /*
177  * Compensation needed?
178  */
179  if (!comp->comp_coils || comp->comp_coils->ncoil <= 0 || !comp->set || !comp->set->current)
180  return OK;
181  /*
182  * Workspace needed?
183  */
184  if (!comp->work)
185  comp->work = MALLOC_60(comp->comp_coils->ncoil,float);
186  /*
187  * Compute the field in the compensation coils
188  */
189  if (comp->field(rd,Q,comp->comp_coils,comp->work,comp->client) == FAIL)
190  return FAIL;
191  /*
192  * Compute the compensated field
193  */
194  return MneCTFCompDataSet::mne_apply_ctf_comp(comp->set,TRUE,res,coils->ncoil,comp->work,comp->comp_coils->ncoil);
195 }
196 
197 //=============================================================================================================
198 
199 void FwdCompData::fwd_free_comp_data(void *d)
200 {
201  FwdCompData* comp = (FwdCompData*)d;
202 
203  if (!comp)
204  return;
205 
206  if (comp->client_free && comp->client)
207  comp->client_free(comp->client);
208 
209  if(comp)
210  delete(comp);
211  return;
212 }
213 
214 //=============================================================================================================
215 
216 int FwdCompData::fwd_make_ctf_comp_coils(MneCTFCompDataSet *set,
217  FwdCoilSet *coils,
218  FwdCoilSet *comp_coils) /* The compensation coil set */
219 /*
220  * Call mne_make_ctf_comp using the information in the coil sets
221  */
222 {
223  QList<FiffChInfo> chs;
224  QList<FiffChInfo> compchs;
225  int nchan = 0;
226  int ncomp = 0;
227  FwdCoil* coil;
228  int k,res;
229 
230  if (!coils || coils->ncoil <= 0) {
231  printf("Coil data missing in fwd_make_ctf_comp_coils");
232  return FAIL;
233  }
234  /*
235  * Create the fake channel info which contain just enough information
236  * for mne_make_ctf_comp
237  */
238  for (k = 0; k < coils->ncoil; k++) {
239  chs.append(FiffChInfo());
240  coil = coils->coils[k];
241  chs[k].ch_name = coil->chname;
242  chs[k].chpos.coil_type = coil->type;
243  chs[k].kind = (coil->coil_class == FWD_COILC_EEG) ? FIFFV_EEG_CH : FIFFV_MEG_CH;
244  }
245  nchan = coils->ncoil;
246  if (comp_coils && comp_coils->ncoil > 0) {
247  for (k = 0; k < comp_coils->ncoil; k++) {
248  compchs.append(FiffChInfo());
249  coil = comp_coils->coils[k];
250  compchs[k].ch_name = coil->chname;
251  compchs[k].chpos.coil_type = coil->type;
252  compchs[k].kind = (coil->coil_class == FWD_COILC_EEG) ? FIFFV_EEG_CH : FIFFV_MEG_CH;
253  }
254  ncomp = comp_coils->ncoil;
255  }
256  res = MneCTFCompDataSet::mne_make_ctf_comp(set,chs,nchan,compchs,ncomp);
257 
258  return res;
259 }
260 
261 //=============================================================================================================
262 
263 FwdCompData *FwdCompData::fwd_make_comp_data(MneCTFCompDataSet *set,
264  FwdCoilSet *coils,
265  FwdCoilSet *comp_coils,
266  fwdFieldFunc field,
267  fwdVecFieldFunc vec_field,
268  fwdFieldGradFunc field_grad,
269  void *client,
270  fwdUserFreeFunc client_free)
271 /*
272  * Compose a compensation data set
273  */
274 {
275  FwdCompData* comp = new FwdCompData();
276 
277  if(set)
278  comp->set = new MneCTFCompDataSet(*set);
279  else
280  comp->set = NULL;
281 
282  if (comp_coils) {
283  comp->comp_coils = comp_coils->dup_coil_set(NULL);
284  }
285  else {
286  qWarning("No coils to duplicate");
287  comp->comp_coils = NULL;
288  }
289  comp->field = field;
290  comp->vec_field = vec_field;
291  comp->field_grad = field_grad;
292  comp->client = client;
293  comp->client_free = client_free;
294 
295  if (fwd_make_ctf_comp_coils(comp->set,
296  coils,
297  comp->comp_coils) != OK) {
298  fwd_free_comp_data(comp);
299  return NULL;
300  }
301  else {
302  return comp;
303  }
304 }
305 
306 //=============================================================================================================
307 
308 int FwdCompData::fwd_comp_field_vec(float *rd, FwdCoilSet *coils, float **res, void *client)
309 /*
310  * Calculate the compensated field (all dipole components)
311  */
312 {
313  FwdCompData* comp = (FwdCompData*)client;
314  int k;
315 
316  if (!comp->vec_field) {
317  printf("Field computation function is missing in fwd_comp_field_vec");
318  return FAIL;
319  }
320  /*
321  * First compute the field in the primary set of coils
322  */
323  if (comp->vec_field(rd,coils,res,comp->client) == FAIL)
324  return FAIL;
325  /*
326  * Compensation needed?
327  */
328  if (!comp->comp_coils || comp->comp_coils->ncoil <= 0 || !comp->set || !comp->set->current)
329  return OK;
330  /*
331  * Need workspace?
332  */
333  if (!comp->vec_work)
334  comp->vec_work = ALLOC_CMATRIX_60(3,comp->comp_coils->ncoil);
335  /*
336  * Compute the field at the compensation sensors
337  */
338  if (comp->vec_field(rd,comp->comp_coils,comp->vec_work,comp->client) == FAIL)
339  return FAIL;
340  /*
341  * Compute the compensated field of three orthogonal dipoles
342  */
343  for (k = 0; k < 3; k++) {
344  if (MneCTFCompDataSet::mne_apply_ctf_comp(comp->set,TRUE,res[k],coils->ncoil,comp->vec_work[k],comp->comp_coils->ncoil) == FAIL)
345  return FAIL;
346  }
347  return OK;
348 }
349 
350 //=============================================================================================================
351 
352 int FwdCompData::fwd_comp_field_grad(float *rd, float *Q, FwdCoilSet* coils, float *res, float *xgrad, float *ygrad, float *zgrad, void *client)
353 /*
354  * Calculate the compensated field (one dipole component)
355  */
356 {
357  FwdCompData* comp = (FwdCompData*)client;
358 
359  if (!comp->field_grad) {
360  qCritical("Field and gradient computation function is missing in fwd_comp_field_grad");
361  return FAIL;
362  }
363  /*
364  * First compute the field in the primary set of coils
365  */
366  if (comp->field_grad(rd,Q,coils,res,xgrad,ygrad,zgrad,comp->client) == FAIL)
367  return FAIL;
368  /*
369  * Compensation needed?
370  */
371  if (!comp->comp_coils || comp->comp_coils->ncoil <= 0 || !comp->set || !comp->set->current)
372  return OK;
373  /*
374  * Workspace needed?
375  */
376  if (!comp->work)
377  comp->work = MALLOC_60(comp->comp_coils->ncoil,float);
378  if (!comp->vec_work)
379  comp->vec_work = ALLOC_CMATRIX_60(3,comp->comp_coils->ncoil);
380  /*
381  * Compute the field in the compensation coils
382  */
383  if (comp->field_grad(rd,Q,comp->comp_coils,comp->work,comp->vec_work[0],comp->vec_work[1],comp->vec_work[2],comp->client) == FAIL)
384  return FAIL;
385  /*
386  * Compute the compensated field
387  */
388  if (MneCTFCompDataSet::mne_apply_ctf_comp(comp->set,TRUE,res,coils->ncoil,comp->work,comp->comp_coils->ncoil) != OK)
389  return FAIL;
390  if (MneCTFCompDataSet::mne_apply_ctf_comp(comp->set,TRUE,xgrad,coils->ncoil,comp->vec_work[0],comp->comp_coils->ncoil) != OK)
391  return FAIL;
392  if (MneCTFCompDataSet::mne_apply_ctf_comp(comp->set,TRUE,ygrad,coils->ncoil,comp->vec_work[1],comp->comp_coils->ncoil) != OK)
393  return FAIL;
394  if (MneCTFCompDataSet::mne_apply_ctf_comp(comp->set,TRUE,zgrad,coils->ncoil,comp->vec_work[2],comp->comp_coils->ncoil) != OK)
395  return FAIL;
396  return OK;
397 }
Old fiff_type declarations - replace them.
Channel info descriptor.
Definition: fiff_ch_info.h:74
FwdCoilSet description.
Definition: fwd_coil_set.h:75
FwdCompData class declaration.
One MNE CTF Compensation Data Set description.
FwdCoil description.
Definition: fwd_coil.h:88
QString chname
Definition: fwd_coil.h:167
MneCTFCompDataSet class declaration.
This structure is used in the compensated field calculations.
Definition: fwd_comp_data.h:87
FwdCoilSet * dup_coil_set(const FIFFLIB::FiffCoordTransOld *t) const