MNE-CPP  0.1.9
A Framework for Electrophysiology
guess_data.cpp
Go to the documentation of this file.
1 //=============================================================================================================
37 //=============================================================================================================
38 // INCLUDES
39 //=============================================================================================================
40 
41 #include "guess_data.h"
42 #include "dipole_fit_data.h"
43 #include "dipole_forward.h"
44 #include <mne/c/mne_surface_old.h>
46 
47 #include <fiff/fiff_stream.h>
48 #include <fiff/fiff_tag.h>
49 
50 #include <QFile>
51 
52 //=============================================================================================================
53 // USED NAMESPACES
54 //=============================================================================================================
55 
56 using namespace Eigen;
57 using namespace FIFFLIB;
58 using namespace MNELIB;
59 using namespace FWDLIB;
60 using namespace INVERSELIB;
61 
62 //ToDo remove later on
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 X_16 0
80 #define Y_16 1
81 #define Z_16 2
82 
83 #define VEC_COPY_16(to,from) {\
84  (to)[X_16] = (from)[X_16];\
85  (to)[Y_16] = (from)[Y_16];\
86  (to)[Z_16] = (from)[Z_16];\
87  }
88 
89 #define MALLOC_16(x,t) (t *)malloc((x)*sizeof(t))
90 
91 #define REALLOC_16(x,y,t) (t *)((x == NULL) ? malloc((y)*sizeof(t)) : realloc((x),(y)*sizeof(t)))
92 
93 #define ALLOC_CMATRIX_16(x,y) mne_cmatrix_16((x),(y))
94 
95 static void matrix_error_16(int kind, int nr, int nc)
96 
97 {
98  if (kind == 1)
99  printf("Failed to allocate memory pointers for a %d x %d matrix\n",nr,nc);
100  else if (kind == 2)
101  printf("Failed to allocate memory for a %d x %d matrix\n",nr,nc);
102  else
103  printf("Allocation error for a %d x %d matrix\n",nr,nc);
104  if (sizeof(void *) == 4) {
105  printf("This is probably because you seem to be using a computer with 32-bit architecture.\n");
106  printf("Please consider moving to a 64-bit platform.");
107  }
108  printf("Cannot continue. Sorry.\n");
109  exit(1);
110 }
111 
112 float **mne_cmatrix_16(int nr,int nc)
113 
114 {
115  int i;
116  float **m;
117  float *whole;
118 
119  m = MALLOC_16(nr,float *);
120  if (!m) matrix_error_16(1,nr,nc);
121  whole = MALLOC_16(nr*nc,float);
122  if (!whole) matrix_error_16(2,nr,nc);
123 
124  for(i=0;i<nr;i++)
125  m[i] = whole + i*nc;
126  return m;
127 }
128 
129 #define FREE_16(x) if ((char *)(x) != NULL) free((char *)(x))
130 #define FREE_CMATRIX_16(m) mne_free_cmatrix_16((m))
131 
132 void mne_free_cmatrix_16 (float **m)
133 {
134  if (m) {
135  FREE_16(*m);
136  FREE_16(m);
137  }
138 }
139 
140 void fromFloatEigenMatrix_16(const Eigen::MatrixXf& from_mat, float **& to_mat, const int m, const int n)
141 {
142  for ( int i = 0; i < m; ++i)
143  for ( int j = 0; j < n; ++j)
144  to_mat[i][j] = from_mat(i,j);
145 }
146 
147 void fromFloatEigenMatrix_16(const Eigen::MatrixXf& from_mat, float **& to_mat)
148 {
149  fromFloatEigenMatrix_16(from_mat, to_mat, from_mat.rows(), from_mat.cols());
150 }
151 
152 //int
153 void fromIntEigenMatrix_16(const Eigen::MatrixXi& from_mat, int **&to_mat, const int m, const int n)
154 {
155  for ( int i = 0; i < m; ++i)
156  for ( int j = 0; j < n; ++j)
157  to_mat[i][j] = from_mat(i,j);
158 }
159 
160 void fromIntEigenMatrix_16(const Eigen::MatrixXi& from_mat, int **&to_mat)
161 {
162  fromIntEigenMatrix_16(from_mat, to_mat, from_mat.rows(), from_mat.cols());
163 }
164 
165 //=============================================================================================================
166 // DEFINE MEMBER METHODS
167 //=============================================================================================================
168 
169 GuessData::GuessData()
170 : rr(NULL)
171 , guess_fwd(NULL)
172 , nguess(0)
173 {
174 }
175 
176 //=============================================================================================================
177 
178 //GuessData::GuessData(const GuessData& p_GuessData)
179 //{
180 //}
181 
182 //=============================================================================================================
183 
184 GuessData::GuessData(const QString &guessname, const QString &guess_surfname, float mindist, float exclude, float grid, DipoleFitData *f)
185 {
186  MneSourceSpaceOld* *sp = NULL;
187  int nsp = 0;
188 // GuessData* res = new GuessData();
189  int k,p;
190  float guessrad = 0.080;
191  MneSourceSpaceOld* guesses = NULL;
192  dipoleFitFuncs orig;
193 
194  if (!guessname.isEmpty()) {
195  /*
196  * Read the guesses and transform to the appropriate coordinate frame
197  */
198  if (MneSurfaceOrVolume::mne_read_source_spaces(guessname,&sp,&nsp) == FAIL)
199  goto bad;
200  if (nsp != 1) {
201  printf("Incorrect number of source spaces in guess file");
202  for (k = 0; k < nsp; k++)
203  delete sp[k];
204  FREE_16(sp);
205  goto bad;
206  }
207  printf("Read guesses from %s\n",guessname.toUtf8().constData());
208  guesses = sp[0]; FREE_16(sp);
209  }
210  else {
211  MneSurfaceOld* inner_skull = NULL;
212  int free_inner_skull = FALSE;
213  float r0[3];
214 
215  VEC_COPY_16(r0,f->r0);
216  FiffCoordTransOld::fiff_coord_trans_inv(r0,f->mri_head_t,TRUE);
217  if (f->bem_model) {
218  printf("Using inner skull surface from the BEM (%s)...\n",f->bemname.toUtf8().constData());
219  if ((inner_skull = f->bem_model->fwd_bem_find_surface(FIFFV_BEM_SURF_ID_BRAIN)) == NULL)
220  goto bad;
221  }
222  else if (!guess_surfname.isEmpty()) {
223  printf("Reading inner skull surface from %s...\n",guess_surfname.toUtf8().data());
224  if ((inner_skull = MneSurfaceOrVolume::read_bem_surface(guess_surfname,FIFFV_BEM_SURF_ID_BRAIN,TRUE,NULL)) == NULL)
225  goto bad;
226  free_inner_skull = TRUE;
227  }
228  if ((guesses = (MneSourceSpaceOld*)FwdBemModel::make_guesses(inner_skull,guessrad,r0,grid,exclude,mindist)) == NULL)
229  goto bad;
230  if (free_inner_skull)
231  delete inner_skull;
232  }
233  if (MneSurfaceOrVolume::mne_transform_source_spaces_to(f->coord_frame,f->mri_head_t,&guesses,1) != OK)
234  goto bad;
235  printf("Guess locations are now in %s coordinates.\n",FiffCoordTransOld::mne_coord_frame_name(f->coord_frame));
236  this->nguess = guesses->nuse;
237  this->rr = ALLOC_CMATRIX_16(guesses->nuse,3);
238  for (k = 0, p = 0; k < guesses->np; k++)
239  if (guesses->inuse[k]) {
240  VEC_COPY_16(this->rr[p],guesses->rr[k]);
241  p++;
242  }
243  delete guesses; guesses = NULL;
244 
245  printf("Go through all guess source locations...");
246  this->guess_fwd = MALLOC_16(this->nguess,DipoleForward*);
247  for (k = 0; k < this->nguess; k++)
248  this->guess_fwd[k] = NULL;
249  /*
250  * Compute the guesses using the sphere model for speed
251  */
252  orig = f->funcs;
253  if (f->fit_mag_dipoles)
254  f->funcs = f->mag_dipole_funcs;
255  else
256  f->funcs = f->sphere_funcs;
257 
258  for (k = 0; k < this->nguess; k++) {
259  if ((this->guess_fwd[k] = DipoleFitData::dipole_forward_one(f,this->rr[k],NULL)) == NULL)
260  goto bad;
261 #ifdef DEBUG
262  sing = this->guess_fwd[k]->sing;
263  printf("%f %f %f\n",sing[0],sing[1],sing[2]);
264 #endif
265  }
266  f->funcs = orig;
267 
268  printf("[done %d sources]\n",p);
269 
270  return;
271 // return res;
272 
273 bad : {
274  if(guesses)
275  delete guesses;
276 
277  return;
278 // return NULL;
279  }
280 }
281 
282 //=============================================================================================================
283 
284 GuessData::GuessData(const QString &guessname, const QString &guess_surfname, float mindist, float exclude, float grid, DipoleFitData *f, char *guess_save_name)
285 {
286  MneSourceSpaceOld* *sp = NULL;
287  int nsp = 0;
288  int k,p;
289  float guessrad = 0.080f;
290  MneSourceSpaceOld* guesses = NULL;
291 
292  if (!guessname.isEmpty()) {
293  /*
294  * Read the guesses and transform to the appropriate coordinate frame
295  */
296  if (MneSurfaceOrVolume::mne_read_source_spaces(guessname,&sp,&nsp) == FIFF_FAIL)
297  goto bad;
298  if (nsp != 1) {
299  qCritical("Incorrect number of source spaces in guess file");
300  for (k = 0; k < nsp; k++)
301  delete sp[k];
302  FREE_16(sp);
303  goto bad;
304  }
305  printf("Read guesses from %s\n",guessname.toUtf8().constData());
306  guesses = sp[0]; FREE_16(sp);
307  }
308  else {
309  MneSurfaceOld* inner_skull = NULL;
310  int free_inner_skull = FALSE;
311  float r0[3];
312 
313  VEC_COPY_16(r0,f->r0);
314  FiffCoordTransOld::fiff_coord_trans_inv(r0,f->mri_head_t,TRUE);
315  if (f->bem_model) {
316  printf("Using inner skull surface from the BEM (%s)...\n",f->bemname.toUtf8().constData());
317  if ((inner_skull = f->bem_model->fwd_bem_find_surface(FIFFV_BEM_SURF_ID_BRAIN)) == NULL)
318  goto bad;
319  }
320  else if (!guess_surfname.isEmpty()) {
321  printf("Reading inner skull surface from %s...\n",guess_surfname.toUtf8().data());
322  if ((inner_skull = MneSurfaceOrVolume::read_bem_surface(guess_surfname,FIFFV_BEM_SURF_ID_BRAIN,TRUE,NULL)) == NULL)
323  goto bad;
324  free_inner_skull = TRUE;
325  }
326  if ((guesses = (MneSourceSpaceOld*)FwdBemModel::make_guesses(inner_skull,guessrad,r0,grid,exclude,mindist)) == NULL)
327  goto bad;
328  if (free_inner_skull)
329  delete inner_skull;
330  }
331  /*
332  * Save the guesses for future use
333  */
334  if (guesses->nuse == 0) {
335  qCritical("No active guess locations remaining.");
336  goto bad;
337  }
338  if (guess_save_name) {
339  printf("###################DEBUG writing source spaces not yet implemented.");
340  // if (mne_write_source_spaces(guess_save_name,&guesses,1,FALSE) != OK)
341  // goto bad;
342  // printf("Wrote guess locations to %s\n",guess_save_name);
343  }
344  /*
345  * Transform the guess locations to the appropriate coordinate frame
346  */
347  if (MneSurfaceOrVolume::mne_transform_source_spaces_to(f->coord_frame,f->mri_head_t,&guesses,1) != OK)
348  goto bad;
349  printf("Guess locations are now in %s coordinates.\n",FiffCoordTransOld::mne_coord_frame_name(f->coord_frame));
350 
351  this->nguess = guesses->nuse;
352  this->rr = ALLOC_CMATRIX_16(guesses->nuse,3);
353  for (k = 0, p = 0; k < guesses->np; k++)
354  if (guesses->inuse[k]) {
355  VEC_COPY_16(this->rr[p],guesses->rr[k]);
356  p++;
357  }
358  if(guesses)
359  delete guesses;
360  guesses = NULL;
361 
362  this->guess_fwd = MALLOC_16(this->nguess,DipoleForward*);
363  for (k = 0; k < this->nguess; k++)
364  this->guess_fwd[k] = NULL;
365  /*
366  * Compute the guesses using the sphere model for speed
367  */
368  if (!this->compute_guess_fields(f))
369  goto bad;
370 
371  return;
372 
373 bad : {
374  if(guesses)
375  delete guesses;
376  return;
377  }
378 }
379 
380 //=============================================================================================================
381 
383 {
384  FREE_CMATRIX_16(rr);
385  if (guess_fwd) {
386  for (int k = 0; k < nguess; k++)
387  delete guess_fwd[k];
388  FREE_16(guess_fwd);
389  }
390  return;
391 }
392 
393 //=============================================================================================================
394 
396 {
397  dipoleFitFuncs orig = NULL;
398 
399  if (!f) {
400  qCritical("Data missing in compute_guess_fields");
401  return false;
402  }
403  if (!f->noise) {
404  qCritical("Noise covariance missing in compute_guess_fields");
405  return false;
406  }
407  printf("Go through all guess source locations...");
408  orig = f->funcs;
409  if (f->fit_mag_dipoles)
410  f->funcs = f->mag_dipole_funcs;
411  else
412  f->funcs = f->sphere_funcs;
413  for (int k = 0; k < this->nguess; k++) {
414  if ((this->guess_fwd[k] = DipoleFitData::dipole_forward_one(f,this->rr[k],this->guess_fwd[k])) == NULL){
415  if (orig)
416  f->funcs = orig;
417  return false;
418  }
419 #ifdef DEBUG
420  sing = this->guess_fwd[k]->sing;
421  printf("%f %f %f\n",sing[0],sing[1],sing[2]);
422 #endif
423  }
424  f->funcs = orig;
425  printf("[done %d sources]\n",this->nguess);
426 
427  return true;
428 }
dipoleFitFuncs mag_dipole_funcs
GuessData class declaration.
FiffStream class declaration.
MNELIB::MneCovMatrix * noise
Dipole Fit Data implementation.
DipoleForward description.
DipoleForward class declaration.
This defines a surface.
This defines a source space.
DipoleForward ** guess_fwd
Definition: guess_data.h:139
FWDLIB::FwdBemModel * bem_model
FIFFLIB::FiffCoordTransOld * mri_head_t
FiffTag class declaration, which provides fiff tag I/O and processing methods.
MneSurfaceOld class declaration.
MneSourceSpaceOld class declaration.
Dipole Fit Data class declaration.
bool compute_guess_fields(DipoleFitData *f)
Definition: guess_data.cpp:395