MNE-CPP  0.1.9
A Framework for Electrophysiology
fwd_bem_model.h
Go to the documentation of this file.
1 //=============================================================================================================
37 #ifndef FWDBEMMODEL_H
38 #define FWDBEMMODEL_H
39 
40 //=============================================================================================================
41 // INCLUDES
42 //=============================================================================================================
43 
44 #include "fwd_global.h"
45 #include "fwd_coil_set.h"
46 
48 #include <fiff/fiff_dir_node.h>
49 #include <fiff/fiff_tag.h>
50 #include <fiff/fiff_named_matrix.h>
51 
52 //=============================================================================================================
53 // EIGEN INCLUDES
54 //=============================================================================================================
55 
56 #include <Eigen/Core>
57 
58 //=============================================================================================================
59 // QT INCLUDES
60 //=============================================================================================================
61 
62 #include <QSharedPointer>
63 #include <QString>
64 
65 #define FWD_BEM_UNKNOWN -1
66 #define FWD_BEM_CONSTANT_COLL 1
67 #define FWD_BEM_LINEAR_COLL 2
68 
69 #define FWD_BEM_IP_APPROACH_LIMIT 0.1
70 
71 #define FWD_BEM_LIN_FIELD_SIMPLE 1
72 #define FWD_BEM_LIN_FIELD_FERGUSON 2
73 #define FWD_BEM_LIN_FIELD_URANKAR 3
74 
75 //=============================================================================================================
76 // FORWARD DECLARATIONS
77 //=============================================================================================================
78 
79 namespace MNELIB
80 {
81  class MneTriangle;
82  class MneSurfaceOld;
83  class MneSourceSpaceOld;
84  class MneCTFCompDataSet;
85  class MneNamedMatrix;
86 }
87 
88 namespace FIFFLIB {
89  class FiffNamedMatrix;
90 }
91 //=============================================================================================================
92 // DEFINE NAMESPACE FWDLIB
93 //=============================================================================================================
94 
95 namespace FWDLIB
96 {
97 
98 //=============================================================================================================
99 // FWDLIB FORWARD DECLARATIONS
100 //=============================================================================================================
101 
102 class FwdEegSphereModel;
103 
104 //=============================================================================================================
111 {
112 public:
113  typedef QSharedPointer<FwdBemModel> SPtr;
114  typedef QSharedPointer<const FwdBemModel> ConstSPtr;
116  //=========================================================================================================
122  explicit FwdBemModel();
123 
124  //=========================================================================================================
129  virtual ~FwdBemModel();
130 
131  //=========================================================================================================
135  void fwd_bem_free_solution();
136 
137 //char *fwd_bem_make_bem_name(char *name)
139 // * Make a standard BEM file name
140 // */
141 //{
142 // char *s1,*s2;
143 
144 // s1 = strip_from(name,(char*)(".fif"));
145 // s2 = strip_from(s1,(char*)("-sol"));
146 // FREE_3(s1);
147 // s1 = strip_from(s2,(char*)("-bem"));
148 // FREE_3(s2);
149 // s2 = MALLOC_3(strlen(s1)+strlen(BEM_SUFFIX)+1,char);
150 // sprintf(s2,"%s%s",s1,BEM_SUFFIX);
151 // FREE_3(s1);
152 // return s2;
153 //}
154 
155  static QString fwd_bem_make_bem_sol_name(const QString& name);
156 
157  //============================= fwd_bem_model.c =============================
158 
159  static const QString& fwd_bem_explain_surface(int kind);
160 
161  static const QString& fwd_bem_explain_method(int method);
162 
163  static int get_int( FIFFLIB::FiffStream::SPtr& stream, const FIFFLIB::FiffDirNode::SPtr& node,int what,int *res);
164 
165  /*
166  * Return a pointer to a specific surface in a BEM
167  */
168  // Refactored: fwd_bem_find_surface (fwd_bem_model.c)
169  MNELIB::MneSurfaceOld* fwd_bem_find_surface(int kind);
170 
171  static FwdBemModel* fwd_bem_load_surfaces(const QString& name,
172  int *kinds,
173  int nkind);
174 
175  static FwdBemModel* fwd_bem_load_homog_surface(const QString& name);
176 
177  static FwdBemModel* fwd_bem_load_three_layer_surfaces(const QString& name);
178 
179  static int fwd_bem_load_solution(const QString& name, int bem_method, FwdBemModel* m);
180 
181  static int fwd_bem_set_head_mri_t(FwdBemModel* m, FIFFLIB::FiffCoordTransOld* t);
182 
183  //============================= dipole_fit_guesses.c =============================
184 
185  static MNELIB::MneSurfaceOld* make_guesses(MNELIB::MneSurfaceOld* guess_surf, /* Predefined boundary for the guesses */
186  float guessrad, /* Radius for the spherical boundary if the above is missing */
187  float *guess_r0, /* Origin for the spherical boundary */
188  float grid, /* Spacing between guess points */
189  float exclude, /* Exclude points closer than this to the CM of the guess boundary surface */
190  float mindist);
191 
192  //============================= fwd_bem_linear_collocation.c =============================
193 
194  /*
195  * The following approach is based on:
196  *
197  * de Munck JC: "A linear discretization of the volume conductor boundary integral equation using analytically integrated elements",
198  * IEEE Trans Biomed Eng. 1992 39(9) : 986 - 990
199  *
200  */
201 
202  static double calc_beta (double *rk,double *rk1);
203 
204  static void lin_pot_coeff (float *from, /* Origin */
205  MNELIB::MneTriangle* to, /* The destination triangle */
206  double omega[3]);
207 
208  static void correct_auto_elements (MNELIB::MneSurfaceOld* surf,
209  float **mat);
210 
211  static float **fwd_bem_lin_pot_coeff (const QList<MNELIB::MneSurfaceOld*>& surfs);
212 
213  static int fwd_bem_linear_collocation_solution(FwdBemModel* m);
214 
215  //============================= fwd_bem_solution.c =============================
216 
217  static float **fwd_bem_multi_solution (float **solids, /* The solid-angle matrix */
218  float **gamma, /* The conductivity multipliers */
219  int nsurf, /* Number of surfaces */
220  int *ntri);
221 
222  static float **fwd_bem_homog_solution (float **solids,int ntri);
223 
224  static void fwd_bem_ip_modify_solution(float **solution, /* The original solution */
225  float **ip_solution, /* The isolated problem solution */
226  float ip_mult, /* Conductivity ratio */
227  int nsurf, /* Number of surfaces */
228  int *ntri);
229 
230  //============================= fwd_bem_constant_collocation.c =============================
231 
232  static int fwd_bem_check_solids (float **angles,int ntri1,int ntri2, float desired);
233 
234  static float **fwd_bem_solid_angles (const QList<MNELIB::MneSurfaceOld*>& surfs);
235 
236  static int fwd_bem_constant_collocation_solution(FwdBemModel* m);
237 
238  //============================= fwd_bem_model.c =============================
239 
240  static int fwd_bem_compute_solution(FwdBemModel* m,
241  int bem_method);
242 
243  static int fwd_bem_load_recompute_solution(const QString& name,
244  int bem_method,
245  int force_recompute,
246  FwdBemModel* m);
247 
248  //============================= fwd_bem_pot.c =============================
249 
250  static float fwd_bem_inf_field(float *rd, /* Dipole position */
251  float *Q, /* Dipole moment */
252  float *rp, /* Field point */
253  float *dir);
254 
255  static float fwd_bem_inf_pot (float *rd, /* Dipole position */
256  float *Q, /* Dipole moment */
257  float *rp);
258 
259  static int fwd_bem_specify_els(FwdBemModel* m,
260  FwdCoilSet* els);
261 
262  static void fwd_bem_pot_grad_calc (float *rd, /* Dipole position */
263  float *Q, /* Dipole orientation */
264  FwdBemModel* m, /* The model */
265  FwdCoilSet* els, /* Use this electrode set if available */
266  int all_surfs, /* Compute solution on all surfaces? */
267  float *xgrad,
268  float *ygrad,
269  float *zgrad);
270 
271  static void fwd_bem_lin_pot_calc (float *rd, /* Dipole position */
272  float *Q, /* Dipole orientation */
273  FwdBemModel* m, /* The model */
274  FwdCoilSet* els, /* Use this electrode set if available */
275  int all_surfs, /* Compute on all surfaces? */
276  float *pot);
277 
278  static void fwd_bem_lin_pot_grad_calc (float *rd, /* Dipole position */
279  float *Q, /* Dipole orientation */
280  FwdBemModel* m, /* The model */
281  FwdCoilSet* els, /* Use this electrode set if available */
282  int all_surfs, /* Compute on all surfaces? */
283  float *xgrad,
284  float *ygrad,
285  float *zgrad);
286 
287  static void fwd_bem_pot_calc (float *rd, /* Dipole position */
288  float *Q, /* Dipole orientation */
289  FwdBemModel* m, /* The model */
290  FwdCoilSet* els, /* Use this electrode set if available */
291  int all_surfs, /* Compute solution on all surfaces? */
292  float *pot);
293 
294  static int fwd_bem_pot_els (float *rd, /* Dipole position */
295  float *Q, /* Dipole orientation */
296  FwdCoilSet* els, /* Electrode descriptors */
297  float *pot, /* Result */
298  void *client);
299 
300  static int fwd_bem_pot_grad_els (float *rd, /* Dipole position */
301  float *Q, /* Dipole orientation */
302  FwdCoilSet* els, /* Electrode descriptors */
303  float *pot, /* Potentials */
304  float *xgrad, /* Derivatives with respect to dipole position */
305  float *ygrad,
306  float *zgrad,
307  void *client);
308 
309  //============================= fwd_bem_field.c =============================
310 
311  /*
312  * These are some of the integration formulas listed in
313  *
314  * L. Urankar, Common compact analytical formulas for computation of
315  * geometry integrals on a basic Cartesian sub-domain in boundary and
316  * volume integral methods, Engineering Analysis with Boundary Elements,
317  * 7 (3), 1990, 124 - 129.
318  *
319  */
320 
321  static void calc_f (double *xx,
322  double *yy, /* Corner coordinates */
323  double *f0,
324  double *fx,
325  double *fy);
326 
327  static void calc_magic (double u,double z,
328  double A,
329  double B,
330  double *beta,
331  double *D);
332 
333  static void field_integrals (float *from,
335  double *I1p,
336  double *T,double *S1,double *S2,
337  double *f0,double *fx,double *fy);
338 
339  static double one_field_coeff (float *dest, /* The destination field point */
340  float *normal, /* The field direction we are interested in */
341  MNELIB::MneTriangle* tri);
342 
343  static float **fwd_bem_field_coeff(FwdBemModel* m, /* The model */
344  FwdCoilSet* coils);
345 
346  /*
347  * These are the formulas from Ferguson et al
348  * A Complete Linear Discretization for Calculating the Magnetic Field
349  * Using the Boundary-Element Method, IEEE Trans. Biomed. Eng., submitted
350  */
351  static double calc_gamma (double *rk,double *rk1);
352 
353  static void fwd_bem_one_lin_field_coeff_ferg (float *dest, /* The field point */
354  float *dir, /* The interesting direction */
355  MNELIB::MneTriangle* tri, /* The destination triangle */
356  double *res);
357 
358  static void fwd_bem_one_lin_field_coeff_uran(float *dest, /* The field point */
359  float *dir, /* The interesting direction */
360  MNELIB::MneTriangle* tri, /* The destination triangle */
361  double *res);
362 
363  static void fwd_bem_one_lin_field_coeff_simple (float *dest, /* The destination field point */
364  float *normal, /* The field direction we are interested in */
365  MNELIB::MneTriangle* source, /* The source triangle */
366  double *res);
367 
368  typedef void (* linFieldIntFunc)(float *dest,float *dir,MNELIB::MneTriangle* tri, double *res);
369 
370  static float **fwd_bem_lin_field_coeff (FwdBemModel* m, /* The model */
371  FwdCoilSet* coils, /* Coil information */
372  int method);
373 
374  static int fwd_bem_specify_coils(FwdBemModel* m,
375  FwdCoilSet* coils);
376 
377  #define MAG_FACTOR 1e-7 /* \mu_0/4\pi */
378 
379  static void fwd_bem_lin_field_calc(float *rd,
380  float *Q,
381  FwdCoilSet* coils,
382  FwdBemModel* m,
383  float *B);
384 
385  static void fwd_bem_field_calc(float *rd,
386  float *Q,
387  FwdCoilSet* coils,
388  FwdBemModel* m,
389  float *B);
390 
391  static void fwd_bem_field_grad_calc(float *rd,
392  float *Q,
393  FwdCoilSet *coils,
394  FwdBemModel *m,
395  float *xgrad,
396  float *ygrad,
397  float *zgrad);
398 
399  static float fwd_bem_inf_field_der(float *rd, /* Dipole position */
400  float *Q, /* Dipole moment */
401  float *rp, /* Field point */
402  float *dir, /* Which field component */
403  float *comp);
404 
405  static float fwd_bem_inf_pot_der (float *rd, /* Dipole position */
406  float *Q, /* Dipole moment */
407  float *rp, /* Potential point */
408  float *comp);
409 
410  static void fwd_bem_lin_field_grad_calc(float *rd,
411  float *Q,
412  FwdCoilSet *coils,
413  FwdBemModel *m,
414  float *xgrad,
415  float *ygrad,
416  float *zgrad);
417 
418  static int fwd_bem_field(float *rd, /* Dipole position */
419  float *Q, /* Dipole orientation */
420  FwdCoilSet* coils, /* Coil descriptors */
421  float *B, /* Result */
422  void *client);
423 
424  static int fwd_bem_field_grad(float *rd, /* The dipole location */
425  float Q[], /* The dipole components (xyz) */
426  FwdCoilSet* coils, /* The coil definitions */
427  float Bval[], /* Results */
428  float xgrad[], /* The derivatives with respect to */
429  float ygrad[], /* the dipole position coordinates */
430  float zgrad[],
431  void *client);
432 
433  //============================= compute_forward.c =============================
434 
435  static void *meg_eeg_fwd_one_source_space(void *arg);
436 
437  // TODO check if this is the correct class or move
438  static int compute_forward_meg( MNELIB::MneSourceSpaceOld* *spaces,
439  int nspace,
440  FwdCoilSet* coils,
441  FwdCoilSet* comp_coils,
442  MNELIB::MneCTFCompDataSet* comp_data,
443  bool fixed_ori,
444  FwdBemModel* bem_model,
445  Eigen::Vector3f* r0,
446  bool use_threads,
448  FIFFLIB::FiffNamedMatrix& resp_grad,
449  bool bDoGRad);
451  static int compute_forward_eeg( MNELIB::MneSourceSpaceOld* *spaces,
452  int nspace,
453  FwdCoilSet* els,
454  bool fixed_ori,
455  FwdBemModel* bem_model,
456  FwdEegSphereModel* m,
457  bool use_threads,
459  FIFFLIB::FiffNamedMatrix& resp_grad,
460  bool bDoGrad);
462  //============================= fwd_spherefield.c =============================
463  // TODO location of these functions need to be checked -> evtl moving to more suitable space
464  static int fwd_sphere_field(float *rd, /* The dipole location */
465  float Q[], /* The dipole components (xyz) */
466  FwdCoilSet* coils, /* The coil definitions */
467  float Bval[], /* Results */
468  void *client);
469 
470  static int fwd_sphere_field_vec(float *rd, /* The dipole location */
471  FwdCoilSet* coils, /* The coil definitions */
472  float **Bval, /* Results: rows are the fields of the x,y, and z direction dipoles */
473  void *client);
474 
475  static int fwd_sphere_field_grad(float *rd, /* The dipole location */
476  float Q[], /* The dipole components (xyz) */
477  FwdCoilSet* coils, /* The coil definitions */
478  float Bval[], /* Results */
479  float xgrad[], /* The derivatives with respect to */
480  float ygrad[], /* the dipole position coordinates */
481  float zgrad[],
482  void *client);
483 
484  //============================= fwd_mag_dipole_field.c =============================
485  // TODO location of these functions need to be checked -> evtl moving to mor suitable space
486  /*
487  * Compute the field of a magnetic dipole
488  */
489  static int fwd_mag_dipole_field( float *rm, /* The dipole location in the same coordinate system as the coils */
490  float M[], /* The dipole components (xyz) */
491  FwdCoilSet* coils, /* The coil definitions */
492  float Bval[], /* Results */
493  void *client);
494 
495  static int fwd_mag_dipole_field_vec( float *rm, /* The dipole location */
496  FwdCoilSet* coils, /* The coil definitions */
497  float **Bval, /* Results: rows are the fields of the x,y, and z direction dipoles */
498  void *client);
499 
500 public:
501  QString surf_name; /* Name of the file where surfaces were loaded from */
502  QList<MNELIB::MneSurfaceOld*> surfs; /* The interface surfaces from outside towards inside */
503  int *ntri; /* Number of triangles on each surface */
504  int *np; /* Number of vertices on each surface */
505  int nsurf; /* How many */
506  float *sigma; /* The conductivities */
507  float **gamma; /* The gamma factors */
508  float *source_mult; /* These multiply the infinite medium potentials */
509  float *field_mult; /* Multipliers for the magnetic field */
510  int bem_method; /* Which approximation method is used */
511  QString sol_name; /* Name of the file where the solution was loaded from */
512 
513  float **solution; /* The potential solution matrix */
514  float *v0; /* Space for the infinite-medium potentials */
515  int nsol; /* Size of the solution matrix */
516 
517  FIFFLIB::FiffCoordTransOld* head_mri_t; /* Coordinate transformation from head to MRI coordinates */
518 
519  float ip_approach_limit; /* Controls whether we need to use the isolated problem approach */
520  bool use_ip_approach; /* Do we need it */
521 
522 // ### OLD STRUCT ###
523 //typedef struct {
524 // char *surf_name; /* Name of the file where surfaces were loaded from */
525 // FWDLIB::MneSurfaceOld* *surfs; /* The interface surfaces from outside towards inside */
526 // int *ntri; /* Number of triangles on each surface */
527 // int *np; /* Number of vertices on each surface */
528 // int nsurf; /* How many */
529 // float *sigma; /* The conductivities */
530 // float **gamma; /* The gamma factors */
531 // float *source_mult; /* These multiply the infinite medium potentials */
532 // float *field_mult; /* Multipliers for the magnetic field */
533 // int bem_method; /* Which approximation method is used */
534 // char *sol_name; /* Name of the file where the solution was loaded from */
535 
536 // float **solution; /* The potential solution matrix */
537 // float *v0; /* Space for the infinite-medium potentials */
538 // int nsol; /* Size of the solution matrix */
539 
540 // FWDLIB::FiffCoordTransOld* head_mri_t; /* Coordinate transformation from head to MRI coordinates */
541 
542 // float ip_approach_limit; /* Controls whether we need to use the isolated problem approach */
543 // int use_ip_approach; /* Do we need it */
544 //} *FwdBemModel*,FwdBemModel*Rec; /* Holds the BEM model definition */
545 };
546 
547 //=============================================================================================================
548 // INLINE DEFINITIONS
549 //=============================================================================================================
550 } // NAMESPACE FWDLIB
551 
552 #endif // FWDBEMMODEL_H
QSharedPointer< FiffDirNode > SPtr
Definition: fiff_dir_node.h:77
QSharedPointer< const FwdBemModel > ConstSPtr
FwdCoilSet description.
Definition: fwd_coil_set.h:75
forward library export/import macros.
Holds the BEM model definition.
Triangle data.
Definition: mne_triangle.h:75
This defines a surface.
This defines a source space.
One MNE CTF Compensation Data Set description.
Coordinate transformation descriptor.
FiffNamedMatrix class declaration.
FiffDirNode class declaration, which provides fiff dir tree processing methods.
QSharedPointer< FwdBemModel > SPtr
FwdCoilSet class declaration.
QSharedPointer< FiffStream > SPtr
Definition: fiff_stream.h:107
FiffTag class declaration, which provides fiff tag I/O and processing methods.
#define FWDSHARED_EXPORT
Definition: fwd_global.h:57
Electric Current Dipole description.
FiffCoordTransOld class declaration.