MNE-CPP  0.1.9
A Framework for Electrophysiology
fwd_bem_model.cpp
Go to the documentation of this file.
1 //=============================================================================================================
37 //=============================================================================================================
38 // INCLUDES
39 //=============================================================================================================
40 
41 #include "fwd_bem_model.h"
42 #include "fwd_bem_solution.h"
43 #include "fwd_eeg_sphere_model.h"
44 #include <mne/c/mne_surface_old.h>
45 #include <mne/c/mne_triangle.h>
47 
48 #include "fwd_comp_data.h"
49 #include "fwd_bem_model.h"
50 
51 #include "fwd_thread_arg.h"
52 
53 #include <fiff/fiff_stream.h>
54 #include <fiff/fiff_named_matrix.h>
55 
56 #include <QFile>
57 #include <QList>
58 #include <QThread>
59 #include <QtConcurrent>
60 
61 #define _USE_MATH_DEFINES
62 #include <math.h>
63 
64 #include <Eigen/Dense>
65 
66 static float Qx[] = {1.0,0.0,0.0};
67 static float Qy[] = {0.0,1.0,0.0};
68 static float Qz[] = {0.0,0.0,1.0};
69 
70 #ifndef TRUE
71 #define TRUE 1
72 #endif
73 
74 #ifndef FALSE
75 #define FALSE 0
76 #endif
77 
78 #ifndef FAIL
79 #define FAIL -1
80 #endif
81 
82 #ifndef OK
83 #define OK 0
84 #endif
85 
86 #define X_40 0
87 #define Y_40 1
88 #define Z_40 2
89 
90 #define VEC_DIFF_40(from,to,diff) {\
91  (diff)[X_40] = (to)[X_40] - (from)[X_40];\
92  (diff)[Y_40] = (to)[Y_40] - (from)[Y_40];\
93  (diff)[Z_40] = (to)[Z_40] - (from)[Z_40];\
94  }
95 
96 #define VEC_COPY_40(to,from) {\
97  (to)[X_40] = (from)[X_40];\
98  (to)[Y_40] = (from)[Y_40];\
99  (to)[Z_40] = (from)[Z_40];\
100  }
101 
102 #define VEC_DOT_40(x,y) ((x)[X_40]*(y)[X_40] + (x)[Y_40]*(y)[Y_40] + (x)[Z_40]*(y)[Z_40])
103 
104 #define VEC_LEN_40(x) sqrt(VEC_DOT_40(x,x))
105 
106 #define CROSS_PRODUCT_40(x,y,xy) {\
107  (xy)[X_40] = (x)[Y_40]*(y)[Z_40]-(y)[Y_40]*(x)[Z_40];\
108  (xy)[Y_40] = -((x)[X_40]*(y)[Z_40]-(y)[X_40]*(x)[Z_40]);\
109  (xy)[Z_40] = (x)[X_40]*(y)[Y_40]-(y)[X_40]*(x)[Y_40];\
110  }
111 
112 #define MALLOC_40(x,t) (t *)malloc((x)*sizeof(t))
113 
114 #define ALLOC_CMATRIX_40(x,y) mne_cmatrix_40((x),(y))
115 
116 #define FREE_40(x) if ((char *)(x) != NULL) free((char *)(x))
117 
118 #define FREE_CMATRIX_40(m) mne_free_cmatrix_40((m))
119 
120 void mne_free_cmatrix_40 (float **m)
121 {
122  if (m) {
123  FREE_40(*m);
124  FREE_40(m);
125  }
126 }
127 
128 static void matrix_error_40(int kind, int nr, int nc)
129 
130 {
131  if (kind == 1)
132  printf("Failed to allocate memory pointers for a %d x %d matrix\n",nr,nc);
133  else if (kind == 2)
134  printf("Failed to allocate memory for a %d x %d matrix\n",nr,nc);
135  else
136  printf("Allocation error for a %d x %d matrix\n",nr,nc);
137  if (sizeof(void *) == 4) {
138  printf("This is probably because you seem to be using a computer with 32-bit architecture.\n");
139  printf("Please consider moving to a 64-bit platform.");
140  }
141  printf("Cannot continue. Sorry.\n");
142  exit(1);
143 }
144 
145 float **mne_cmatrix_40 (int nr,int nc)
146 
147 {
148  int i;
149  float **m;
150  float *whole;
151 
152  m = MALLOC_40(nr,float *);
153  if (!m) matrix_error_40(1,nr,nc);
154  whole = MALLOC_40(nr*nc,float);
155  if (!whole) matrix_error_40(2,nr,nc);
156 
157  for(i=0;i<nr;i++)
158  m[i] = whole + i*nc;
159  return m;
160 }
161 
162 //float
163 Eigen::MatrixXf toFloatEigenMatrix_40(float **mat, const int m, const int n)
164 {
165  Eigen::MatrixXf eigen_mat(m,n);
166 
167  for ( int i = 0; i < m; ++i)
168  for ( int j = 0; j < n; ++j)
169  eigen_mat(i,j) = mat[i][j];
170 
171  return eigen_mat;
172 }
173 
174 void fromFloatEigenMatrix_40(const Eigen::MatrixXf& from_mat, float **& to_mat, const int m, const int n)
175 {
176  for ( int i = 0; i < m; ++i)
177  for ( int j = 0; j < n; ++j)
178  to_mat[i][j] = from_mat(i,j);
179 }
180 
181 void fromFloatEigenMatrix_40(const Eigen::MatrixXf& from_mat, float **& to_mat)
182 {
183  fromFloatEigenMatrix_40(from_mat, to_mat, from_mat.rows(), from_mat.cols());
184 }
185 
186 float **mne_lu_invert_40(float **mat,int dim)
187 /*
188  * Invert a matrix using the LU decomposition from
189  * LAPACK
190  */
191 {
192  Eigen::MatrixXf eigen_mat = toFloatEigenMatrix_40(mat, dim, dim);
193  Eigen::MatrixXf eigen_mat_inv = eigen_mat.inverse();
194  fromFloatEigenMatrix_40(eigen_mat_inv, mat);
195  return mat;
196 }
197 
198 void mne_transpose_square_40(float **mat, int n)
199 /*
200  * In-place transpose of a square matrix
201  */
202 {
203  int j,k;
204  float val;
205 
206  for (j = 1; j < n; j++)
207  for (k = 0; k < j; k++) {
208  val = mat[j][k];
209  mat[j][k] = mat[k][j];
210  mat[k][j] = val;
211  }
212  return;
213 }
214 
215 float mne_dot_vectors_40(float *v1,
216  float *v2,
217  int nn)
218 
219 {
220 #ifdef BLAS
221  int one = 1;
222  float res = sdot(&nn,v1,&one,v2,&one);
223  return res;
224 #else
225  float res = 0.0;
226  int k;
227 
228  for (k = 0; k < nn; k++)
229  res = res + v1[k]*v2[k];
230  return res;
231 #endif
232 }
233 
234 void mne_add_scaled_vector_to_40(float *v1,float scale, float *v2,int nn)
235 
236 {
237 #ifdef BLAS
238  float fscale = scale;
239  int one = 1;
240  saxpy(&nn,&fscale,v1,&one,v2,&one);
241 #else
242  int k;
243  for (k = 0; k < nn; k++)
244  v2[k] = v2[k] + scale*v1[k];
245 #endif
246  return;
247 }
248 
249 void mne_scale_vector_40(double scale,float *v,int nn)
250 
251 {
252 #ifdef BLAS
253  float fscale = scale;
254  int one = 1;
255  sscal(&nn,&fscale,v,&one);
256 #else
257  int k;
258  for (k = 0; k < nn; k++)
259  v[k] = v[k]*scale;
260 #endif
261 }
262 
263 #include <Eigen/Core>
264 using namespace Eigen;
265 
266 float **mne_mat_mat_mult_40 (float **m1,
267  float **m2,
268  int d1,
269  int d2,
270  int d3)
271 /* Matrix multiplication
272  * result(d1 x d3) = m1(d1 x d2) * m2(d2 x d3) */
273 
274 {
275 #ifdef BLAS
276  float **result = ALLOC_CMATRIX_40(d1,d3);
277  char *transa = "N";
278  char *transb = "N";
279  float zero = 0.0;
280  float one = 1.0;
281  sgemm (transa,transb,&d3,&d1,&d2,
282  &one,m2[0],&d3,m1[0],&d2,&zero,result[0],&d3);
283  return (result);
284 #else
285  float **result = ALLOC_CMATRIX_40(d1,d3);
286  int j,k,p;
287  float sum;
288 
289  // TODO: Hack to include faster Eigen processing. This should eventually be replaced dwith pure Eigen logic.
290  MatrixXf a(d1,d2);
291  MatrixXf b(d2,d3);
292 
293  // ToDo use Eigen::Map
294  for (j = 0; j < d1; j++) {
295  for (k = 0; k < d2; k++) {
296  a(j,k) = m1[j][k];
297  }
298  }
299 
300  for (j = 0; j < d2; j++) {
301  for (k = 0; k < d3; k++) {
302  b(j,k) = m2[j][k];
303  }
304  }
305 
306  MatrixXf resultMat = a * b;
307 
308  for (j = 0; j < d1; j++) {
309  for (k = 0; k < d3; k++) {
310  result[j][k] = resultMat(j,k);
311  }
312  }
313 
314 // for (j = 0; j < d1; j++)
315 // for (k = 0; k < d3; k++) {
316 // sum = 0.0;
317 // for (p = 0; p < d2; p++)
318 // sum = sum + m1[j][p]*m2[p][k];
319 // result[j][k] = sum;
320 // }
321  return (result);
322 #endif
323 }
324 
325 namespace FWDLIB
326 {
327 
328 static struct {
329  int kind;
330  const QString name;
331 } surf_expl[] = { { FIFFV_BEM_SURF_ID_BRAIN , "inner skull" },
332 { FIFFV_BEM_SURF_ID_SKULL , "outer skull" },
333 { FIFFV_BEM_SURF_ID_HEAD , "scalp" },
334 { -1 , "unknown" } };
335 
336 static struct {
337  int method;
338  const QString name;
339 } method_expl[] = { { FWD_BEM_CONSTANT_COLL , "constant collocation" },
340 { FWD_BEM_LINEAR_COLL , "linear collocation" },
341 { -1 , "unknown" } };
342 
343 }
344 
345 #define BEM_SUFFIX "-bem.fif"
346 #define BEM_SOL_SUFFIX "-bem-sol.fif"
347 
348 //============================= misc_util.c =============================
349 
350 static QString strip_from(const QString& s, const QString& suffix)
351 {
352  QString res;
353 
354  if (s.endsWith(suffix)) {
355  res = s;
356  res.chop(suffix.size());
357  }
358  else
359  res = s;
360 
361  return res;
362 }
363 
364 //============================= mne_coord_transforms.c =============================
365 
366 namespace FWDLIB
367 {
368 
369 typedef struct {
370  int frame;
371  const QString name;
373 
374 }
375 
376 const QString mne_coord_frame_name_40(int frame)
377 
378 {
379  static FWDLIB::frameNameRec_40 frames[] = {
380  {FIFFV_COORD_UNKNOWN,"unknown"},
381  {FIFFV_COORD_DEVICE,"MEG device"},
382  {FIFFV_COORD_ISOTRAK,"isotrak"},
383  {FIFFV_COORD_HPI,"hpi"},
384  {FIFFV_COORD_HEAD,"head"},
385  {FIFFV_COORD_MRI,"MRI (surface RAS)"},
386  {FIFFV_MNE_COORD_MRI_VOXEL, "MRI voxel"},
387  {FIFFV_COORD_MRI_SLICE,"MRI slice"},
388  {FIFFV_COORD_MRI_DISPLAY,"MRI display"},
389  {FIFFV_MNE_COORD_CTF_DEVICE,"CTF MEG device"},
390  {FIFFV_MNE_COORD_CTF_HEAD,"CTF/4D/KIT head"},
391  {FIFFV_MNE_COORD_RAS,"RAS (non-zero origin)"},
392  {FIFFV_MNE_COORD_MNI_TAL,"MNI Talairach"},
393  {FIFFV_MNE_COORD_FS_TAL_GTZ,"Talairach (MNI z > 0)"},
394  {FIFFV_MNE_COORD_FS_TAL_LTZ,"Talairach (MNI z < 0)"},
395  {-1,"unknown"}
396  };
397  int k;
398  for (k = 0; frames[k].frame != -1; k++) {
399  if (frame == frames[k].frame)
400  return frames[k].name;
401  }
402  return frames[k].name;
403 }
404 
405 //=============================================================================================================
406 // USED NAMESPACES
407 //=============================================================================================================
408 
409 using namespace Eigen;
410 using namespace FIFFLIB;
411 using namespace MNELIB;
412 using namespace FWDLIB;
413 
414 //=============================================================================================================
415 // DEFINE MEMBER METHODS
416 //=============================================================================================================
417 
418 FwdBemModel::FwdBemModel()
419 :ntri (NULL)
420 ,np (NULL)
421 ,nsurf (0)
422 ,sigma (NULL)
423 ,gamma (NULL)
424 ,source_mult(NULL)
425 ,field_mult (NULL)
426 ,bem_method (FWD_BEM_UNKNOWN)
427 ,solution (NULL)
428 ,nsol (0)
429 ,head_mri_t (NULL)
430 ,v0 (NULL)
431 ,use_ip_approach(false)
432 ,ip_approach_limit(FWD_BEM_IP_APPROACH_LIMIT)
433 {
434 }
435 
436 //=============================================================================================================
437 
439 {
440  for (int k = 0; k < this->nsurf; k++)
441  delete this->surfs[k];
442  FREE_40(this->ntri);
443  FREE_40(this->np);
444  FREE_40(this->sigma);
445  FREE_40(this->source_mult);
446  FREE_40(this->field_mult);
447  FREE_CMATRIX_40(this->gamma);
448  if(this->head_mri_t)
449  delete this->head_mri_t;
450  this->fwd_bem_free_solution();
451 }
452 
453 //=============================================================================================================
454 
456 {
457  FREE_CMATRIX_40(this->solution); this->solution = NULL;
458  this->sol_name.clear();
459  FREE_40(this->v0); this->v0 = NULL;
460  this->bem_method = FWD_BEM_UNKNOWN;
461  this->nsol = 0;
462 
463  return;
464 }
465 
466 //=============================================================================================================
467 
468 QString FwdBemModel::fwd_bem_make_bem_sol_name(const QString& name)
469 /*
470  * Make a standard BEM solution file name
471  */
472 {
473  QString s1, s2;
474 
475  s1 = strip_from(name,".fif");
476  s2 = strip_from(s1,"-sol");
477  s1 = strip_from(s2,"-bem");
478  s2 = QString("%1%2").arg(s1).arg(BEM_SOL_SUFFIX);
479  return s2;
480 }
481 
482 //=============================================================================================================
483 
484 const QString& FwdBemModel::fwd_bem_explain_surface(int kind)
485 {
486  int k;
487 
488  for (k = 0; surf_expl[k].kind >= 0; k++)
489  if (surf_expl[k].kind == kind)
490  return surf_expl[k].name;
491 
492  return surf_expl[k].name;
493 }
494 
495 //=============================================================================================================
496 
497 const QString& FwdBemModel::fwd_bem_explain_method(int method)
498 
499 {
500  int k;
501 
502  for (k = 0; method_expl[k].method >= 0; k++)
503  if (method_expl[k].method == method)
504  return method_expl[k].name;
505 
506  return method_expl[k].name;
507 }
508 
509 //=============================================================================================================
510 
511 int FwdBemModel::get_int(FiffStream::SPtr &stream, const FiffDirNode::SPtr &node, int what, int *res)
512 /*
513  * Wrapper to get int's
514  */
515 {
516  FiffTag::SPtr t_pTag;
517  if(node->find_tag(stream, what, t_pTag)) {
518  if (t_pTag->getType() != FIFFT_INT) {
519  printf("Expected an integer tag : %d (found data type %d instead)\n",what,t_pTag->getType() );
520  return FAIL;
521  }
522  *res = *t_pTag->toInt();
523  return OK;
524  }
525  return FAIL;
526 }
527 
528 //=============================================================================================================
529 
530 MneSurfaceOld *FwdBemModel::fwd_bem_find_surface(int kind)
531 {
532  // if (!model) {
533  // printf("No model specified for fwd_bem_find_surface");
534  // return NULL;
535  // }
536  for (int k = 0; k < this->nsurf; k++)
537  if (this->surfs[k]->id == kind)
538  return this->surfs[k];
539  printf("Desired surface (%d = %s) not found.", kind,fwd_bem_explain_surface(kind).toUtf8().constData());
540  return NULL;
541 }
542 
543 //=============================================================================================================
544 
545 FwdBemModel *FwdBemModel::fwd_bem_load_surfaces(const QString &name, int *kinds, int nkind)
546 /*
547  * Load a set of surfaces
548  */
549 {
550  QList<MneSurfaceOld*> surfs;// = NULL;
551  float *sigma = NULL;
552  float *sigma1;
553  FwdBemModel *m = NULL;
554  int j,k;
555 
556  if (nkind <= 0) {
557  qCritical("No surfaces specified to fwd_bem_load_surfaces");
558  return NULL;
559  }
560 
561 // surfs = MALLOC_40(nkind,MneSurfaceOld*);
562  sigma = MALLOC_40(nkind,float);
563 // for (k = 0; k < nkind; k++)
564 // surfs[k] = NULL;
565 
566  for (k = 0; k < nkind; k++) {
567  surfs.append(MneSurfaceOld::read_bem_surface(name,kinds[k],TRUE,sigma+k));
568  if (surfs[k] == NULL)
569  goto bad;
570  if ((surfs[k] = MneSurfaceOld::read_bem_surface(name,kinds[k],TRUE,sigma+k)) == NULL)
571  goto bad;
572  if (sigma[k] < 0.0) {
573  qCritical("No conductivity available for surface %s",fwd_bem_explain_surface(kinds[k]).toUtf8().constData());
574  goto bad;
575  }
576  if (surfs[k]->coord_frame != FIFFV_COORD_MRI) { /* We make our life much easier with this */
577  qCritical("Surface %s not specified in MRI coordinates.",fwd_bem_explain_surface(kinds[k]).toUtf8().constData());
578  goto bad;
579  }
580  }
581  m = new FwdBemModel;
582 
583  m->surf_name = name;
584  m->nsurf = nkind;
585  m->surfs = surfs;
586  m->sigma = sigma;
587  m->ntri = MALLOC_40(nkind,int);
588  m->np = MALLOC_40(nkind,int);
589  m->gamma = ALLOC_CMATRIX_40(nkind,nkind);
590  m->source_mult = MALLOC_40(nkind,float);
591  m->field_mult = MALLOC_40(nkind,float);
592  /*
593  * Dirty trick for the zero conductivity outside
594  */
595  sigma1 = MALLOC_40(nkind+1,float);
596  sigma1[0] = 0.0;
597  sigma = sigma1+1;
598  for (k = 0; k < m->nsurf; k++)
599  sigma[k] = m->sigma[k];
600  /*
601  * Gamma factors and multipliers
602  */
603  for (j = 0; j < m->nsurf; j++) {
604  m->ntri[j] = m->surfs[j]->ntri;
605  m->np[j] = m->surfs[j]->np;
606  m->source_mult[j] = 2.0/(sigma[j]+sigma[j-1]);
607  m->field_mult[j] = sigma[j]-sigma[j-1];
608  for (k = 0; k < m->nsurf; k++)
609  m->gamma[j][k] = (sigma[k]-sigma[k-1])/(sigma[j]+sigma[j-1]);
610  }
611  FREE_40(sigma1);
612 
613  return m;
614 
615 bad : {
616  FREE_40(sigma);
617  for (k = 0; k < surfs.size(); k++)
618  delete surfs[k];
619 // FREE_40(surfs);
620  surfs.clear();
621  return NULL;
622  }
623 }
624 
625 //=============================================================================================================
626 
627 FwdBemModel *FwdBemModel::fwd_bem_load_homog_surface(const QString &name)
628 /*
629  * Load surfaces for the homogeneous model
630  */
631 {
632  int kinds[] = { FIFFV_BEM_SURF_ID_BRAIN };
633  int nkind = 1;
634 
635  return fwd_bem_load_surfaces(name,kinds,nkind);
636 }
637 
638 //=============================================================================================================
639 
640 FwdBemModel *FwdBemModel::fwd_bem_load_three_layer_surfaces(const QString &name)
641 /*
642  * Load surfaces for three-layer model
643  */
644 {
645  int kinds[] = { FIFFV_BEM_SURF_ID_HEAD, FIFFV_BEM_SURF_ID_SKULL, FIFFV_BEM_SURF_ID_BRAIN };
646  int nkind = 3;
647 
648  return fwd_bem_load_surfaces(name,kinds,nkind);
649 }
650 
651 //=============================================================================================================
652 
653 int FwdBemModel::fwd_bem_load_solution(const QString &name, int bem_method, FwdBemModel *m)
654 /*
655  * Load the potential solution matrix and attach it to the model:
656  *
657  * return values:
658  *
659  * TRUE found a suitable solution
660  * FALSE did not find a suitable solution
661  * FAIL error in reading the solution
662  *
663  */
664 {
665  QFile file(name);
666  FiffStream::SPtr stream(new FiffStream(&file));
667 
668  float **sol = NULL;
669  FiffDirNode::SPtr bem_node;
670  int method;
671  FiffTag::SPtr t_pTag;
672  int nsol;
673 
674  if(!stream->open())
675  goto not_found;
676 
677  /*
678  * Find the BEM data
679  */
680  {
681  QList<FiffDirNode::SPtr> nodes = stream->dirtree()->dir_tree_find(FIFFB_BEM);
682 
683  if (nodes.size() == 0) {
684  printf ("No BEM data in %s",name.toUtf8().constData());
685  goto not_found;
686  }
687  bem_node = nodes[0];
688  }
689  /*
690  * Approximation method
691  */
692  if (get_int(stream,bem_node,FIFF_BEM_APPROX,&method) != OK)
693  goto not_found;
694  if (method == FIFFV_BEM_APPROX_CONST)
695  method = FWD_BEM_CONSTANT_COLL;
696  else if (method == FIFFV_BEM_APPROX_LINEAR)
697  method = FWD_BEM_LINEAR_COLL;
698  else {
699  printf ("Cannot handle BEM approximation method : %d",method);
700  goto bad;
701  }
702  if (bem_method != FWD_BEM_UNKNOWN && method != bem_method) {
703  printf("Approximation method in file : %d desired : %d",method,bem_method);
704  goto not_found;
705  }
706  {
707  int dim,k;
708 
709  if (!bem_node->find_tag(stream, FIFF_BEM_POT_SOLUTION, t_pTag))
710  goto bad;
711  qint32 ndim;
712  QVector<qint32> dims;
713  t_pTag->getMatrixDimensions(ndim, dims);
714 
715  if (ndim != 2) {
716  printf("Expected a two-dimensional solution matrix instead of a %d dimensional one",ndim);
717  goto bad;
718  }
719  for (k = 0, dim = 0; k < m->nsurf; k++)
720  dim = dim + ((method == FWD_BEM_LINEAR_COLL) ? m->surfs[k]->np : m->surfs[k]->ntri);
721  if (dims[0] != dim || dims[1] != dim) {
722  printf("Expected a %d x %d solution matrix instead of a %d x %d one",dim,dim,dims[0],dims[1]);
723  goto not_found;
724  }
725 
726  MatrixXf tmp_sol = t_pTag->toFloatMatrix().transpose();
727  sol = ALLOC_CMATRIX_40(tmp_sol.rows(),tmp_sol.cols());
728  fromFloatEigenMatrix_40(tmp_sol, sol);
729  nsol = dims[1];
730  }
731  if(m)
733  m->sol_name = name;
734  m->solution = sol;
735  m->nsol = nsol;
736  m->bem_method = method;
737  stream->close();
738 
739  return TRUE;
740 
741 bad : {
742  stream->close();
743  FREE_CMATRIX_40(sol);
744  return FAIL;
745  }
746 
747 not_found : {
748  stream->close();
749  FREE_CMATRIX_40(sol);
750  return FALSE;
751  }
752 }
753 
754 //=============================================================================================================
755 
756 int FwdBemModel::fwd_bem_set_head_mri_t(FwdBemModel *m, FiffCoordTransOld *t)
757 /*
758  * Set the coordinate transformation
759  */
760 {
761  if (t->from == FIFFV_COORD_HEAD && t->to == FIFFV_COORD_MRI) {
762  if(m->head_mri_t)
763  delete m->head_mri_t;
764  m->head_mri_t = new FiffCoordTransOld( *t );
765  return OK;
766  }
767  else if (t->from == FIFFV_COORD_MRI && t->to == FIFFV_COORD_HEAD) {
768  if(m->head_mri_t)
769  delete m->head_mri_t;
770  m->head_mri_t = t->fiff_invert_transform();
771  return OK;
772  }
773  else {
774  printf ("Improper coordinate transform delivered to fwd_bem_set_head_mri_t");
775  return FAIL;
776  }
777 }
778 
779 //=============================================================================================================
780 
781 MneSurfaceOld* FwdBemModel::make_guesses(MneSurfaceOld* guess_surf, float guessrad, float *guess_r0, float grid, float exclude, float mindist) /* Exclude points closer than this to
782  * the guess boundary surface */
783 /*
784  * Make a guess space inside a sphere
785  */
786 {
787  char *bemname = NULL;
788  MneSurfaceOld* sphere = NULL;
789  MneSurfaceOld* res = NULL;
790  int k;
791  float dist;
792  float r0[] = { 0.0, 0.0, 0.0 };
793 
794  if (!guess_r0)
795  guess_r0 = r0;
796 
797  if (!guess_surf) {
798  printf("Making a spherical guess space with radius %7.1f mm...\n",1000*guessrad);
799  //#ifdef USE_SHARE_PATH
800  // if ((bemname = mne_compose_mne_name("share/mne","icos.fif")) == NULL)
801  //#else
802  // if ((bemname = mne_compose_mne_name("setup/mne","icos.fif")) == NULL)
803  //#endif
804  // goto out;
805 
806  // QFile bemFile("/usr/pubsw/packages/mne/stable/share/mne/icos.fif");
807 
808  QFile bemFile(QString(QCoreApplication::applicationDirPath() + "/resources/general/surf2bem/icos.fif"));
809  if ( !QCoreApplication::startingUp() )
810  bemFile.setFileName(QCoreApplication::applicationDirPath() + QString("/resources/general/surf2bem/icos.fif"));
811  else if (!bemFile.exists())
812  bemFile.setFileName("./resources/general/surf2bem/icos.fif");
813 
814  if( !bemFile.exists () ){
815  qDebug() << bemFile.fileName() << "does not exists.";
816  goto out;
817  }
818 
819  bemname = MALLOC_40(strlen(bemFile.fileName().toUtf8().data())+1,char);
820  strcpy(bemname,bemFile.fileName().toUtf8().data());
821 
822  if ((sphere = MneSourceSpaceOld::read_bem_surface(bemname,9003,FALSE,NULL)) == NULL)
823  goto out;
824 
825  for (k = 0; k < sphere->np; k++) {
826  dist = VEC_LEN_40(sphere->rr[k]);
827  sphere->rr[k][X_40] = guessrad*sphere->rr[k][X_40]/dist + guess_r0[X_40];
828  sphere->rr[k][Y_40] = guessrad*sphere->rr[k][Y_40]/dist + guess_r0[Y_40];
829  sphere->rr[k][Z_40] = guessrad*sphere->rr[k][Z_40]/dist + guess_r0[Z_40];
830  }
831  if (MneSurfaceOrVolume::mne_source_space_add_geometry_info((MneSourceSpaceOld*)sphere,TRUE) == FAIL)
832  goto out;
833  guess_surf = sphere;
834  }
835  else {
836  printf("Guess surface (%d = %s) is in %s coordinates\n",
837  guess_surf->id,FwdBemModel::fwd_bem_explain_surface(guess_surf->id).toUtf8().constData(),
838  mne_coord_frame_name_40(guess_surf->coord_frame).toUtf8().constData());
839  }
840  printf("Filtering (grid = %6.f mm)...\n",1000*grid);
841  res = (MneSurfaceOld*)MneSurfaceOrVolume::make_volume_source_space(guess_surf,grid,exclude,mindist);
842 
843 out : {
844  FREE_40(bemname);
845  if(sphere)
846  delete sphere;
847  return res;
848  }
849 }
850 
851 //=============================================================================================================
852 
853 double FwdBemModel::calc_beta(double *rk, double *rk1)
854 
855 {
856  double rkk1[3];
857  double size;
858  double res;
859 
860  VEC_DIFF_40(rk,rk1,rkk1);
861  size = VEC_LEN_40(rkk1);
862 
863  res = log((VEC_LEN_40(rk)*size + VEC_DOT_40(rk,rkk1))/
864  (VEC_LEN_40(rk1)*size + VEC_DOT_40(rk1,rkk1)))/size;
865  return (res);
866 }
867 
868 //=============================================================================================================
869 
870 void FwdBemModel::lin_pot_coeff(float *from, MneTriangle* to, double omega[]) /* The final result */
871 /*
872  * The linear potential matrix element computations
873  */
874 {
875  double y1[3],y2[3],y3[3]; /* Corners with origin at from */
876  double *y[5];
877  double **yy;
878  double l1,l2,l3; /* Lengths of y1, y2, and y3 */
879  double solid; /* The standard solid angle */
880  double vec_omega[3]; /* The cross-product integral */
881  double cross[3]; /* y1 x y2 */
882  double triple; /* VEC_DOT_40(y1 x y2,y3) */
883  double ss;
884  double beta[3],bbeta[3];
885  int j,k;
886  double z[3];
887  double n2,area2;
888  double diff[3];
889  static const double solid_eps = 4.0*M_PI/1.0E6;
890  /*
891  * This circularity makes things easy for us...
892  */
893  y[0] = y3;
894  y[1] = y1;
895  y[2] = y2;
896  y[3] = y3;
897  y[4] = y1;
898  yy = y + 1; /* yy can have index -1! */
899  /*
900  * The standard solid angle computation
901  */
902  VEC_DIFF_40(from,to->r1,y1);
903  VEC_DIFF_40(from,to->r2,y2);
904  VEC_DIFF_40(from,to->r3,y3);
905 
906  CROSS_PRODUCT_40(y1,y2,cross);
907  triple = VEC_DOT_40(cross,y3);
908 
909  l1 = VEC_LEN_40(y1);
910  l2 = VEC_LEN_40(y2);
911  l3 = VEC_LEN_40(y3);
912  ss = (l1*l2*l3+VEC_DOT_40(y1,y2)*l3+VEC_DOT_40(y1,y3)*l2+VEC_DOT_40(y2,y3)*l1);
913  solid = 2.0*atan2(triple,ss);
914  if (std::fabs(solid) < solid_eps) {
915  for (k = 0; k < 3; k++)
916  omega[k] = 0.0;
917  }
918  else {
919  /*
920  * Calculate the magic vector vec_omega
921  */
922  for (j = 0; j < 3; j++)
923  beta[j] = calc_beta(yy[j],yy[j+1]);
924  bbeta[0] = beta[2] - beta[0];
925  bbeta[1] = beta[0] - beta[1];
926  bbeta[2] = beta[1] - beta[2];
927 
928  for (j = 0; j < 3; j++)
929  vec_omega[j] = 0.0;
930  for (j = 0; j < 3; j++)
931  for (k = 0; k < 3; k++)
932  vec_omega[k] = vec_omega[k] + bbeta[j]*yy[j][k];
933  /*
934  * Put it all together...
935  */
936  area2 = 2.0*to->area;
937  n2 = 1.0/(area2*area2);
938  for (k = 0; k < 3; k++) {
939  CROSS_PRODUCT_40(yy[k+1],yy[k-1],z);
940  VEC_DIFF_40(yy[k+1],yy[k-1],diff);
941  omega[k] = n2*(-area2*VEC_DOT_40(z,to->nn)*solid +
942  triple*VEC_DOT_40(diff,vec_omega));
943  }
944  }
945 #ifdef CHECK
946  /*
947  * Check it out!
948  *
949  * omega1 + omega2 + omega3 = solid
950  */
951  rel1 = (solid + omega[X_40]+omega[Y_40]+omega[Z_40])/solid;
952  /*
953  * The other way of evaluating...
954  */
955  for (j = 0; j < 3; j++)
956  check[j] = 0;
957  for (k = 0; k < 3; k++) {
958  CROSS_PRODUCT_40(to->nn[to],yy[k],z);
959  for (j = 0; j < 3; j++)
960  check[j] = check[j] + omega[k]*z[j];
961  }
962  for (j = 0; j < 3; j++)
963  check[j] = -area2*check[j]/triple;
964  fprintf (stderr,"(%g,%g,%g) =? (%g,%g,%g)\n",
965  check[X_40],check[Y_40],check[Z_40],
966  vec_omega[X_40],vec_omega[Y_40],vec_omega[Z_40]);
967  for (j = 0; j < 3; j++)
968  check[j] = check[j] - vec_omega[j];
969  rel2 = sqrt(VEC_DOT_40(check,check)/VEC_DOT_40(vec_omega,vec_omega));
970  fprintf (stderr,"err1 = %g, err2 = %g\n",100*rel1,100*rel2);
971 #endif
972  return;
973 }
974 
975 //=============================================================================================================
976 
977 void FwdBemModel::correct_auto_elements(MneSurfaceOld *surf, float **mat)
978 /*
979  * Improve auto-element approximation...
980  */
981 {
982  float *row;
983  float sum,miss;
984  int nnode = surf->np;
985  int ntri = surf->ntri;
986  int nmemb;
987  int j,k;
988  float pi2 = 2.0*M_PI;
989  MneTriangle* tri;
990 
991 #ifdef SIMPLE
992  for (j = 0; j < nnode; j++) {
993  row = mat[j];
994  sum = 0.0;
995  for (k = 0; k < nnode; k++)
996  sum = sum + row[k];
997  fprintf (stderr,"row %d sum = %g missing = %g\n",j+1,sum/pi2,
998  1.0-sum/pi2);
999  row[j] = pi2 - sum;
1000  }
1001 #else
1002  for (j = 0; j < nnode; j++) {
1003  /*
1004  * How much is missing?
1005  */
1006  row = mat[j];
1007  sum = 0.0;
1008  for (k = 0; k < nnode; k++)
1009  sum = sum + row[k];
1010  miss = pi2-sum;
1011  nmemb = surf->nneighbor_tri[j];
1012  /*
1013  * The node itself receives one half
1014  */
1015  row[j] = miss/2.0;
1016  /*
1017  * The rest is divided evenly among the member nodes...
1018  */
1019  miss = miss/(4.0*nmemb);
1020  for (k = 0,tri = surf->tris; k < ntri; k++,tri++) {
1021  if (tri->vert[0] == j) {
1022  row[tri->vert[1]] = row[tri->vert[1]] + miss;
1023  row[tri->vert[2]] = row[tri->vert[2]] + miss;
1024  }
1025  else if (tri->vert[1] == j) {
1026  row[tri->vert[0]] = row[tri->vert[0]] + miss;
1027  row[tri->vert[2]] = row[tri->vert[2]] + miss;
1028  }
1029  else if (tri->vert[2] == j) {
1030  row[tri->vert[0]] = row[tri->vert[0]] + miss;
1031  row[tri->vert[1]] = row[tri->vert[1]] + miss;
1032  }
1033  }
1034  /*
1035  * Just check it it out...
1036  *
1037  for (k = 0, sum = 0; k < nnode; k++)
1038  sum = sum + row[k];
1039  fprintf (stderr,"row %d sum = %g\n",j+1,sum/pi2);
1040  */
1041  }
1042 #endif
1043  return;
1044 }
1045 
1046 //=============================================================================================================
1047 
1048 float **FwdBemModel::fwd_bem_lin_pot_coeff(const QList<MneSurfaceOld*>& surfs)
1049 /*
1050  * Calculate the coefficients for linear collocation approach
1051  */
1052 {
1053  float **mat = NULL;
1054  float **sub_mat = NULL;
1055  int np1,np2,ntri,np_tot,np_max;
1056  float **nodes;
1057  MneTriangle* tri;
1058  double omega[3];
1059  double *row = NULL;
1060  int j,k,p,q,c;
1061  int joff,koff;
1062  MneSurfaceOld* surf1;
1063  MneSurfaceOld* surf2;
1064 
1065  for (p = 0, np_tot = np_max = 0; p < surfs.size(); p++) {
1066  np_tot += surfs[p]->np;
1067  if (surfs[p]->np > np_max)
1068  np_max = surfs[p]->np;
1069  }
1070 
1071  mat = ALLOC_CMATRIX_40(np_tot,np_tot);
1072  for (j = 0; j < np_tot; j++)
1073  for (k = 0; k < np_tot; k++)
1074  mat[j][k] = 0.0;
1075  row = MALLOC_40(np_max,double);
1076  sub_mat = MALLOC_40(np_max,float *);
1077  for (p = 0, joff = 0; p < surfs.size(); p++, joff = joff + np1) {
1078  surf1 = surfs[p];
1079  np1 = surf1->np;
1080  nodes = surf1->rr;
1081  for (q = 0, koff = 0; q < surfs.size(); q++, koff = koff + np2) {
1082  surf2 = surfs[q];
1083  np2 = surf2->np;
1084  ntri = surf2->ntri;
1085 
1086  printf("\t\t%s (%d) -> %s (%d) ... ",
1087  fwd_bem_explain_surface(surf1->id).toUtf8().constData(),np1,
1088  fwd_bem_explain_surface(surf2->id).toUtf8().constData(),np2);
1089 
1090  for (j = 0; j < np1; j++) {
1091  for (k = 0; k < np2; k++)
1092  row[k] = 0.0;
1093  for (k = 0, tri = surf2->tris; k < ntri; k++,tri++) {
1094  /*
1095  * No contribution from a triangle that
1096  * this vertex belongs to
1097  */
1098  if (p == q && (tri->vert[0] == j || tri->vert[1] == j || tri->vert[2] == j))
1099  continue;
1100  /*
1101  * Otherwise do the hard job
1102  */
1103  lin_pot_coeff (nodes[j],tri,omega);
1104  for (c = 0; c < 3; c++)
1105  row[tri->vert[c]] = row[tri->vert[c]] - omega[c];
1106  }
1107  for (k = 0; k < np2; k++)
1108  mat[j+joff][k+koff] = row[k];
1109  }
1110  if (p == q) {
1111  for (j = 0; j < np1; j++)
1112  sub_mat[j] = mat[j+joff]+koff;
1113  correct_auto_elements (surf1,sub_mat);
1114  }
1115  printf("[done]\n");
1116  }
1117  }
1118  FREE_40(row);
1119  FREE_40(sub_mat);
1120  return(mat);
1121 }
1122 
1123 //=============================================================================================================
1124 
1125 int FwdBemModel::fwd_bem_linear_collocation_solution(FwdBemModel *m)
1126 /*
1127  * Compute the linear collocation potential solution
1128  */
1129 {
1130  float **coeff = NULL;
1131  float ip_mult;
1132  int k;
1133 
1134  if(m)
1135  m->fwd_bem_free_solution();
1136 
1137  printf("\nComputing the linear collocation solution...\n");
1138  fprintf (stderr,"\tMatrix coefficients...\n");
1139  if ((coeff = fwd_bem_lin_pot_coeff (m->surfs)) == NULL)
1140  goto bad;
1141 
1142  for (k = 0, m->nsol = 0; k < m->nsurf; k++)
1143  m->nsol += m->surfs[k]->np;
1144 
1145  fprintf (stderr,"\tInverting the coefficient matrix...\n");
1146  if ((m->solution = fwd_bem_multi_solution (coeff,m->gamma,m->nsurf,m->np)) == NULL)
1147  goto bad;
1148 
1149  /*
1150  * IP approach?
1151  */
1152  if ((m->nsurf == 3) &&
1153  (ip_mult = m->sigma[m->nsurf-2]/m->sigma[m->nsurf-1]) <= m->ip_approach_limit) {
1154  float **ip_solution = NULL;
1155 
1156  fprintf (stderr,"IP approach required...\n");
1157 
1158  fprintf (stderr,"\tMatrix coefficients (homog)...\n");
1159  QList<MneSurfaceOld*> last_surfs;
1160  last_surfs << m->surfs.last();
1161  if ((coeff = fwd_bem_lin_pot_coeff(last_surfs))== NULL)//m->surfs+m->nsurf-1,1)) == NULL)
1162  goto bad;
1163 
1164  fprintf (stderr,"\tInverting the coefficient matrix (homog)...\n");
1165  if ((ip_solution = fwd_bem_homog_solution (coeff,m->surfs[m->nsurf-1]->np)) == NULL)
1166  goto bad;
1167 
1168  fprintf (stderr,"\tModify the original solution to incorporate IP approach...\n");
1169 
1170  fwd_bem_ip_modify_solution(m->solution,ip_solution,ip_mult,m->nsurf,m->np);
1171  FREE_CMATRIX_40(ip_solution);
1172 
1173  }
1174  m->bem_method = FWD_BEM_LINEAR_COLL;
1175  printf("Solution ready.\n");
1176  return OK;
1177 
1178 bad : {
1179  if(m)
1180  m->fwd_bem_free_solution();
1181  FREE_CMATRIX_40(coeff);
1182  return FAIL;
1183  }
1184 }
1185 
1186 //=============================================================================================================
1187 
1188 float **FwdBemModel::fwd_bem_multi_solution(float **solids, float **gamma, int nsurf, int *ntri) /* Number of triangles or nodes on each surface */
1189 /*
1190  * Invert I - solids/(2*M_PI)
1191  * Take deflation into account
1192  * The matrix is destroyed after inversion
1193  * This is the general multilayer case
1194  */
1195 {
1196  int j,k,p,q;
1197  float defl;
1198  float pi2 = 1.0/(2*M_PI);
1199  float mult;
1200  int joff,koff,jup,kup,ntot;
1201 
1202  for (j = 0,ntot = 0; j < nsurf; j++)
1203  ntot += ntri[j];
1204  defl = 1.0/ntot;
1205  /*
1206  * Modify the matrix
1207  */
1208  for (p = 0, joff = 0; p < nsurf; p++) {
1209  jup = ntri[p] + joff;
1210  for (q = 0, koff = 0; q < nsurf; q++) {
1211  kup = ntri[q] + koff;
1212  mult = (gamma == NULL) ? pi2 : pi2*gamma[p][q];
1213  for (j = joff; j < jup; j++)
1214  for (k = koff; k < kup; k++)
1215  solids[j][k] = defl - solids[j][k]*mult;
1216  koff = kup;
1217  }
1218  joff = jup;
1219  }
1220  for (k = 0; k < ntot; k++)
1221  solids[k][k] = solids[k][k] + 1.0;
1222 
1223  return (mne_lu_invert_40(solids,ntot));
1224 }
1225 
1226 //=============================================================================================================
1227 
1228 float **FwdBemModel::fwd_bem_homog_solution(float **solids, int ntri)
1229 /*
1230  * Invert I - solids/(2*M_PI)
1231  * Take deflation into account
1232  * The matrix is destroyed after inversion
1233  * This is the homogeneous model case
1234  */
1235 {
1236  return fwd_bem_multi_solution (solids,NULL,1,&ntri);
1237 }
1238 
1239 //=============================================================================================================
1240 
1241 void FwdBemModel::fwd_bem_ip_modify_solution(float **solution, float **ip_solution, float ip_mult, int nsurf, int *ntri) /* Number of triangles (nodes) on each surface */
1242 /*
1243  * Modify the solution according to the IP approach
1244  */
1245 {
1246  int s;
1247  int j,k,joff,koff,ntot,nlast;
1248  float mult;
1249  float *row = NULL;
1250  float **sub = NULL;
1251 
1252  for (s = 0, koff = 0; s < nsurf-1; s++)
1253  koff = koff + ntri[s];
1254  nlast = ntri[nsurf-1];
1255  ntot = koff + nlast;
1256 
1257  row = MALLOC_40(nlast,float);
1258  sub = MALLOC_40(ntot,float *);
1259  mult = (1.0 + ip_mult)/ip_mult;
1260 
1261  printf("\t\tCombining...");
1262 #ifndef OLD
1263  printf("t ");
1264  mne_transpose_square_40(ip_solution,nlast);
1265 #endif
1266  for (s = 0, joff = 0; s < nsurf; s++) {
1267  printf("%d3 ",s+1);
1268  /*
1269  * Pick the correct submatrix
1270  */
1271  for (j = 0; j < ntri[s]; j++)
1272  sub[j] = solution[j+joff]+koff;
1273  /*
1274  * Multiply
1275  */
1276 #ifdef OLD
1277  for (j = 0; j < ntri[s]; j++) {
1278  for (k = 0; k < nlast; k++) {
1279  res = mne_dot_vectors_skip_skip(sub[j],1,ip_solution[0]+k,nlast,nlast);
1280  row[k] = sub[j][k] - 2.0*res;
1281  }
1282  for (k = 0; k < nlast; k++)
1283  sub[j][k] = row[k];
1284  }
1285 #else
1286  for (j = 0; j < ntri[s]; j++) {
1287  for (k = 0; k < nlast; k++)
1288  row[k] = mne_dot_vectors_40(sub[j],ip_solution[k],nlast);
1289  mne_add_scaled_vector_to_40(row,-2.0,sub[j],nlast);
1290  }
1291 #endif
1292  joff = joff+ntri[s];
1293  }
1294 #ifndef OLD
1295  printf("t ");
1296  mne_transpose_square_40(ip_solution,nlast);
1297 #endif
1298  printf("33 ");
1299  /*
1300  * The lower right corner is a special case
1301  */
1302  for (j = 0; j < nlast; j++)
1303  for (k = 0; k < nlast; k++)
1304  sub[j][k] = sub[j][k] + mult*ip_solution[j][k];
1305  /*
1306  * Final scaling
1307  */
1308  printf("done.\n\t\tScaling...");
1309  mne_scale_vector_40(ip_mult,solution[0],ntot*ntot);
1310  printf("done.\n");
1311  FREE_40(row); FREE_40(sub);
1312  return;
1313 }
1314 
1315 //=============================================================================================================
1316 
1317 int FwdBemModel::fwd_bem_check_solids(float **angles, int ntri1, int ntri2, float desired)
1318 /*
1319  * Check the angle computations
1320  */
1321 {
1322  float *sums = MALLOC_40(ntri1,float);
1323  float sum;
1324  int j,k;
1325  int res = 0;
1326 
1327  for (j = 0; j < ntri1; j++) {
1328  sum = 0;
1329  for (k = 0; k < ntri2; k++)
1330  sum = sum + angles[j][k];
1331  sums[j] = sum/(2*M_PI);
1332  }
1333  for (j = 0; j < ntri1; j++)
1334  /*
1335  * Three cases:
1336  * same surface: sum = 2*pi
1337  * to outer: sum = 4*pi
1338  * to inner: sum = 0*pi;
1339  */
1340  if (std::fabs(sums[j]-desired) > 1e-4) {
1341  printf("solid angle matrix: rowsum[%d] = 2PI*%g",
1342  j+1,sums[j]);
1343  res = -1;
1344  break;
1345  }
1346  FREE_40(sums);
1347  return res;
1348 }
1349 
1350 //=============================================================================================================
1351 
1352 float **FwdBemModel::fwd_bem_solid_angles(const QList<MneSurfaceOld*>& surfs)
1353 /*
1354  * Compute the solid angle matrix
1355  */
1356 {
1357  MneSurfaceOld* surf1;
1358  MneSurfaceOld* surf2;
1359  MneTriangle* tri;
1360  int ntri1,ntri2,ntri_tot;
1361  int j,k,p,q;
1362  int joff,koff;
1363  float **solids;
1364  float result;
1365  float **sub_solids = NULL;
1366  float desired;
1367 
1368  for (p = 0,ntri_tot = 0; p < surfs.size(); p++)
1369  ntri_tot += surfs[p]->ntri;
1370 
1371  sub_solids = MALLOC_40(ntri_tot,float *);
1372  solids = ALLOC_CMATRIX_40(ntri_tot,ntri_tot);
1373  for (p = 0, joff = 0; p < surfs.size(); p++, joff = joff + ntri1) {
1374  surf1 = surfs[p];
1375  ntri1 = surf1->ntri;
1376  for (q = 0, koff = 0; q < surfs.size(); q++, koff = koff + ntri2) {
1377  surf2 = surfs[q];
1378  ntri2 = surf2->ntri;
1379  printf("\t\t%s (%d) -> %s (%d) ... ",fwd_bem_explain_surface(surf1->id).toUtf8().constData(),ntri1,fwd_bem_explain_surface(surf2->id).toUtf8().constData(),ntri2);
1380  for (j = 0; j < ntri1; j++)
1381  for (k = 0, tri = surf2->tris; k < ntri2; k++, tri++) {
1382  if (p == q && j == k)
1383  result = 0.0;
1384  else
1385  result = MneSurfaceOrVolume::solid_angle (surf1->tris[j].cent,tri);
1386  solids[j+joff][k+koff] = result;
1387  }
1388  for (j = 0; j < ntri1; j++)
1389  sub_solids[j] = solids[j+joff]+koff;
1390  printf("[done]\n");
1391  if (p == q)
1392  desired = 1;
1393  else if (p < q)
1394  desired = 0;
1395  else
1396  desired = 2;
1397  if (fwd_bem_check_solids(sub_solids,ntri1,ntri2,desired) == FAIL) {
1398  FREE_CMATRIX_40(solids);
1399  FREE_40(sub_solids);
1400  return NULL;
1401  }
1402  }
1403  }
1404  FREE_40(sub_solids);
1405  return (solids);
1406 }
1407 
1408 //=============================================================================================================
1409 
1410 int FwdBemModel::fwd_bem_constant_collocation_solution(FwdBemModel *m)
1411 /*
1412  * Compute the solution for the constant collocation approach
1413  */
1414 {
1415  float **solids = NULL;
1416  int k;
1417  float ip_mult;
1418 
1419  if(m)
1420  m->fwd_bem_free_solution();
1421 
1422  printf("\nComputing the constant collocation solution...\n");
1423  printf("\tSolid angles...\n");
1424  if ((solids = fwd_bem_solid_angles(m->surfs)) == NULL)
1425  goto bad;
1426 
1427  for (k = 0, m->nsol = 0; k < m->nsurf; k++)
1428  m->nsol += m->surfs[k]->ntri;
1429 
1430  fprintf (stderr,"\tInverting the coefficient matrix...\n");
1431  if ((m->solution = fwd_bem_multi_solution (solids,m->gamma,m->nsurf,m->ntri)) == NULL)
1432  goto bad;
1433  /*
1434  * IP approach?
1435  */
1436  if ((m->nsurf == 3) &&
1437  (ip_mult = m->sigma[m->nsurf-2]/m->sigma[m->nsurf-1]) <= m->ip_approach_limit) {
1438  float **ip_solution = NULL;
1439 
1440  fprintf (stderr,"IP approach required...\n");
1441 
1442  fprintf (stderr,"\tSolid angles (homog)...\n");
1443  QList<MneSurfaceOld*> last_surfs;
1444  last_surfs << m->surfs.last();
1445  if ((solids = fwd_bem_solid_angles (last_surfs)) == NULL)//m->surfs+m->nsurf-1,1)) == NULL)
1446  goto bad;
1447 
1448  fprintf (stderr,"\tInverting the coefficient matrix (homog)...\n");
1449  if ((ip_solution = fwd_bem_homog_solution (solids,m->surfs[m->nsurf-1]->ntri)) == NULL)
1450  goto bad;
1451 
1452  fprintf (stderr,"\tModify the original solution to incorporate IP approach...\n");
1453  fwd_bem_ip_modify_solution(m->solution,ip_solution,ip_mult,m->nsurf,m->ntri);
1454  FREE_CMATRIX_40(ip_solution);
1455  }
1456  m->bem_method = FWD_BEM_CONSTANT_COLL;
1457  fprintf (stderr,"Solution ready.\n");
1458 
1459  return OK;
1460 
1461 bad : {
1462  if(m)
1463  m->fwd_bem_free_solution();
1464  FREE_CMATRIX_40(solids);
1465  return FAIL;
1466  }
1467 }
1468 
1469 //=============================================================================================================
1470 
1471 int FwdBemModel::fwd_bem_compute_solution(FwdBemModel *m, int bem_method)
1472 /*
1473  * Compute the solution
1474  */
1475 {
1476  /*
1477  * Compute the solution
1478  */
1479  if (bem_method == FWD_BEM_LINEAR_COLL)
1480  return fwd_bem_linear_collocation_solution(m);
1481  else if (bem_method == FWD_BEM_CONSTANT_COLL)
1482  return fwd_bem_constant_collocation_solution(m);
1483 
1484  if(m)
1485  m->fwd_bem_free_solution();
1486  printf ("Unknown BEM method: %d\n",bem_method);
1487  return FAIL;
1488 }
1489 
1490 //=============================================================================================================
1491 
1492 int FwdBemModel::fwd_bem_load_recompute_solution(const QString& name, int bem_method, int force_recompute, FwdBemModel *m)
1493 /*
1494  * Load or recompute the potential solution matrix
1495  */
1496 {
1497  int solres;
1498 
1499  if (!m) {
1500  printf ("No model specified for fwd_bem_load_recompute_solution");
1501  return FAIL;
1502  }
1503 
1504  if (!force_recompute) {
1505  if(m)
1506  m->fwd_bem_free_solution();
1507  solres = fwd_bem_load_solution(name,bem_method,m);
1508  if (solres == TRUE) {
1509  printf("\nLoaded %s BEM solution from %s\n",fwd_bem_explain_method(m->bem_method).toUtf8().constData(),name.toUtf8().constData());
1510  return OK;
1511  }
1512  else if (solres == FAIL)
1513  return FAIL;
1514 #ifdef DEBUG
1515  else
1516  printf("Desired BEM solution not available in %s (%s)\n",name,err_get_error());
1517 #endif
1518  }
1519  if (bem_method == FWD_BEM_UNKNOWN)
1520  bem_method = FWD_BEM_LINEAR_COLL;
1521  return fwd_bem_compute_solution(m,bem_method);
1522 }
1523 
1524 //=============================================================================================================
1525 
1526 float FwdBemModel::fwd_bem_inf_field(float *rd, float *Q, float *rp, float *dir) /* Which field component */
1527 /*
1528  * Infinite-medium magnetic field
1529  * (without \mu_0/4\pi)
1530  */
1531 {
1532  float diff[3],diff2,cross[3];
1533 
1534  VEC_DIFF_40(rd,rp,diff);
1535  diff2 = VEC_DOT_40(diff,diff);
1536  CROSS_PRODUCT_40(Q,diff,cross);
1537 
1538  return (VEC_DOT_40(cross,dir)/(diff2*sqrt(diff2)));
1539 }
1540 
1541 //=============================================================================================================
1542 
1543 float FwdBemModel::fwd_bem_inf_pot(float *rd, float *Q, float *rp) /* Potential point */
1544 /*
1545  * The infinite medium potential
1546  */
1547 {
1548  float diff[3];
1549  float diff2;
1550  VEC_DIFF_40(rd,rp,diff);
1551  diff2 = VEC_DOT_40(diff,diff);
1552  return (VEC_DOT_40(Q,diff)/(4.0*M_PI*diff2*sqrt(diff2)));
1553 }
1554 
1555 //=============================================================================================================
1556 
1557 int FwdBemModel::fwd_bem_specify_els(FwdBemModel* m, FwdCoilSet *els)
1558 /*
1559  * Set up for computing the solution at a set of electrodes
1560  */
1561 {
1562  FwdCoil* el;
1563  MneSurfaceOld* scalp;
1564  int k,p,q,v;
1565  float *one_sol,*pick_sol;
1566  float r[3],w[3],dist;
1567  int best;
1568  MneTriangle* tri;
1569  float x,y,z;
1570  FwdBemSolution* sol;
1571 
1572  if (!m) {
1573  printf("Model missing in fwd_bem_specify_els");
1574  goto bad;
1575  }
1576  if (!m->solution) {
1577  printf("Solution not computed in fwd_bem_specify_els");
1578  goto bad;
1579  }
1580  if (!els || els->ncoil == 0)
1581  return OK;
1582  els->fwd_free_coil_set_user_data();
1583  /*
1584  * Hard work follows
1585  */
1586  els->user_data = sol = new FwdBemSolution();
1587  els->user_data_free = FwdBemSolution::fwd_bem_free_coil_solution;
1588 
1589  sol->ncoil = els->ncoil;
1590  sol->np = m->nsol;
1591  sol->solution = ALLOC_CMATRIX_40(sol->ncoil,sol->np);
1592  /*
1593  * Go through all coils
1594  */
1595  for (k = 0; k < els->ncoil; k++) {
1596  el = els->coils[k];
1597  one_sol = sol->solution[k];
1598  for (q = 0; q < m->nsol; q++)
1599  one_sol[q] = 0.0;
1600  scalp = m->surfs[0];
1601  /*
1602  * Go through all 'integration points'
1603  */
1604  for (p = 0; p < el->np; p++) {
1605  VEC_COPY_40(r,el->rmag[p]);
1606  if (m->head_mri_t != NULL)
1607  FiffCoordTransOld::fiff_coord_trans(r,m->head_mri_t,FIFFV_MOVE);
1608  best = MneSurfaceOrVolume::mne_project_to_surface(scalp,NULL,r,FALSE,&dist);
1609  if (best < 0) {
1610  printf("One of the electrodes could not be projected onto the scalp surface. How come?");
1611  goto bad;
1612  }
1613  if (m->bem_method == FWD_BEM_CONSTANT_COLL) {
1614  /*
1615  * Simply pick the value at the triangle
1616  */
1617  pick_sol = m->solution[best];
1618  for (q = 0; q < m->nsol; q++)
1619  one_sol[q] += el->w[p]*pick_sol[q];
1620  }
1621  else if (m->bem_method == FWD_BEM_LINEAR_COLL) {
1622  /*
1623  * Calculate a linear interpolation between the vertex values
1624  */
1625  tri = scalp->tris+best;
1626  MneSurfaceOrVolume::mne_triangle_coords(r,scalp,best,&x,&y,&z);
1627 
1628  w[X_40] = el->w[p]*(1.0 - x - y);
1629  w[Y_40] = el->w[p]*x;
1630  w[Z_40] = el->w[p]*y;
1631  for (v = 0; v < 3; v++) {
1632  pick_sol = m->solution[tri->vert[v]];
1633  for (q = 0; q < m->nsol; q++)
1634  one_sol[q] += w[v]*pick_sol[q];
1635  }
1636  }
1637  else {
1638  printf("Unknown BEM approximation method : %d\n",m->bem_method);
1639  goto bad;
1640  }
1641  }
1642  }
1643  return OK;
1644 
1645 bad : {
1646  els->fwd_free_coil_set_user_data();
1647  return FAIL;
1648  }
1649 }
1650 
1651 //=============================================================================================================
1652 
1653 void FwdBemModel::fwd_bem_pot_grad_calc(float *rd, float *Q, FwdBemModel* m, FwdCoilSet* els, int all_surfs, float *xgrad, float *ygrad, float *zgrad)
1654 /*
1655  * Compute the potentials due to a current dipole
1656  */
1657 {
1658  MneTriangle* tri;
1659  int ntri;
1660  int s,k,p,nsol,pp;
1661  float mult;
1662  float *v0,ee[3];
1663  float **solution;
1664  float mri_rd[3],mri_Q[3];
1665 
1666  float *grads[3];
1667  float *grad;
1668 
1669  grads[0] = xgrad;
1670  grads[1] = ygrad;
1671  grads[2] = zgrad;
1672 
1673  if (!m->v0)
1674  m->v0 = MALLOC_40(m->nsol,float);
1675  v0 = m->v0;
1676 
1677  VEC_COPY_40(mri_rd,rd);
1678  VEC_COPY_40(mri_Q,Q);
1679  if (m->head_mri_t) {
1680  FiffCoordTransOld::fiff_coord_trans(mri_rd,m->head_mri_t,FIFFV_MOVE);
1681  FiffCoordTransOld::fiff_coord_trans(mri_Q,m->head_mri_t,FIFFV_NO_MOVE);
1682  }
1683  for (pp = 0; pp < 3; pp++) {
1684  grad = grads[pp];
1685 
1686  for (p = 0; p < 3; p++)
1687  ee[p] = p == pp ? 1.0 : 0.0;
1688  if (m->head_mri_t)
1689  FiffCoordTransOld::fiff_coord_trans(ee,m->head_mri_t,FIFFV_NO_MOVE);
1690 
1691  for (s = 0, p = 0; s < m->nsurf; s++) {
1692  ntri = m->surfs[s]->ntri;
1693  tri = m->surfs[s]->tris;
1694  mult = m->source_mult[s];
1695  for (k = 0; k < ntri; k++, tri++)
1696  v0[p++] = mult*fwd_bem_inf_pot_der(mri_rd,mri_Q,tri->cent,ee);
1697  }
1698  if (els) {
1699  FwdBemSolution* sol = (FwdBemSolution*)els->user_data;
1700  solution = sol->solution;
1701  nsol = sol->ncoil;
1702  }
1703  else {
1704  solution = m->solution;
1705  nsol = all_surfs ? m->nsol : m->surfs[0]->ntri;
1706  }
1707  for (k = 0; k < nsol; k++)
1708  grad[k] = mne_dot_vectors_40(solution[k],v0,m->nsol);
1709  }
1710  return;
1711 }
1712 
1713 //=============================================================================================================
1714 
1715 void FwdBemModel::fwd_bem_lin_pot_calc(float *rd, float *Q, FwdBemModel *m, FwdCoilSet *els, int all_surfs, float *pot) /* Put the result here */
1716 /*
1717  * Compute the potentials due to a current dipole
1718  * using the linear potential approximation
1719  */
1720 {
1721  float **rr;
1722  int np;
1723  int s,k,p,nsol;
1724  float mult,mri_rd[3],mri_Q[3];
1725 
1726  float *v0;
1727  float **solution;
1728 
1729  if (!m->v0)
1730  m->v0 = MALLOC_40(m->nsol,float);
1731  v0 = m->v0;
1732 
1733  VEC_COPY_40(mri_rd,rd);
1734  VEC_COPY_40(mri_Q,Q);
1735  if (m->head_mri_t) {
1736  FiffCoordTransOld::fiff_coord_trans(mri_rd,m->head_mri_t,FIFFV_MOVE);
1737  FiffCoordTransOld::fiff_coord_trans(mri_Q,m->head_mri_t,FIFFV_NO_MOVE);
1738  }
1739  for (s = 0, p = 0; s < m->nsurf; s++) {
1740  np = m->surfs[s]->np;
1741  rr = m->surfs[s]->rr;
1742  mult = m->source_mult[s];
1743  for (k = 0; k < np; k++)
1744  v0[p++] = mult*fwd_bem_inf_pot(mri_rd,mri_Q,rr[k]);
1745  }
1746  if (els) {
1747  FwdBemSolution* sol = (FwdBemSolution*)els->user_data;
1748  solution = sol->solution;
1749  nsol = sol->ncoil;
1750  }
1751  else {
1752  solution = m->solution;
1753  nsol = all_surfs ? m->nsol : m->surfs[0]->np;
1754  }
1755  for (k = 0; k < nsol; k++)
1756  pot[k] = mne_dot_vectors_40(solution[k],v0,m->nsol);
1757  return;
1758 }
1759 
1760 //=============================================================================================================
1761 
1762 void FwdBemModel::fwd_bem_lin_pot_grad_calc(float *rd, float *Q, FwdBemModel *m, FwdCoilSet *els, int all_surfs, float *xgrad, float *ygrad, float *zgrad)
1763 /*
1764  * Compute the derivaties of potentials due to a current dipole with respect to the dipole position
1765  * using the linear potential approximation
1766  */
1767 {
1768  float **rr;
1769  int np;
1770  int s,k,p,nsol,pp;
1771  float mult,mri_rd[3],mri_Q[3];
1772 
1773  float *v0,ee[3];
1774  float **solution;
1775 
1776  float *grads[3];
1777  float *grad;
1778 
1779  grads[0] = xgrad;
1780  grads[1] = ygrad;
1781  grads[2] = zgrad;
1782 
1783  if (!m->v0)
1784  m->v0 = MALLOC_40(m->nsol,float);
1785  v0 = m->v0;
1786 
1787  VEC_COPY_40(mri_rd,rd);
1788  VEC_COPY_40(mri_Q,Q);
1789  if (m->head_mri_t) {
1790  FiffCoordTransOld::fiff_coord_trans(mri_rd,m->head_mri_t,FIFFV_MOVE);
1791  FiffCoordTransOld::fiff_coord_trans(mri_Q,m->head_mri_t,FIFFV_NO_MOVE);
1792  }
1793  for (pp = 0; pp < 3; pp++) {
1794  grad = grads[pp];
1795 
1796  for (p = 0; p < 3; p++)
1797  ee[p] = p == pp ? 1.0 : 0.0;
1798  if (m->head_mri_t)
1799  FiffCoordTransOld::fiff_coord_trans(ee,m->head_mri_t,FIFFV_NO_MOVE);
1800 
1801  for (s = 0, p = 0; s < m->nsurf; s++) {
1802  np = m->surfs[s]->np;
1803  rr = m->surfs[s]->rr;
1804  mult = m->source_mult[s];
1805  for (k = 0; k < np; k++)
1806  v0[p++] = mult*fwd_bem_inf_pot_der(mri_rd,mri_Q,rr[k],ee);
1807  }
1808  if (els) {
1809  FwdBemSolution* sol = (FwdBemSolution*)els->user_data;
1810  solution = sol->solution;
1811  nsol = sol->ncoil;
1812  }
1813  else {
1814  solution = m->solution;
1815  nsol = all_surfs ? m->nsol : m->surfs[0]->np;
1816  }
1817  for (k = 0; k < nsol; k++)
1818  grad[k] = mne_dot_vectors_40(solution[k],v0,m->nsol);
1819  }
1820  return;
1821 }
1822 
1823 //=============================================================================================================
1824 
1825 void FwdBemModel::fwd_bem_pot_calc(float *rd, float *Q, FwdBemModel *m, FwdCoilSet *els, int all_surfs, float *pot)
1826 /*
1827  * Compute the potentials due to a current dipole
1828  */
1829 {
1830  MneTriangle* tri;
1831  int ntri;
1832  int s,k,p,nsol;
1833  float mult;
1834  float *v0;
1835  float **solution;
1836  float mri_rd[3],mri_Q[3];
1837 
1838  if (!m->v0)
1839  m->v0 = MALLOC_40(m->nsol,float);
1840  v0 = m->v0;
1841 
1842  VEC_COPY_40(mri_rd,rd);
1843  VEC_COPY_40(mri_Q,Q);
1844  if (m->head_mri_t) {
1845  FiffCoordTransOld::fiff_coord_trans(mri_rd,m->head_mri_t,FIFFV_MOVE);
1846  FiffCoordTransOld::fiff_coord_trans(mri_Q,m->head_mri_t,FIFFV_NO_MOVE);
1847  }
1848  for (s = 0, p = 0; s < m->nsurf; s++) {
1849  ntri = m->surfs[s]->ntri;
1850  tri = m->surfs[s]->tris;
1851  mult = m->source_mult[s];
1852  for (k = 0; k < ntri; k++, tri++)
1853  v0[p++] = mult*fwd_bem_inf_pot(mri_rd,mri_Q,tri->cent);
1854  }
1855  if (els) {
1856  FwdBemSolution* sol = (FwdBemSolution*)els->user_data;
1857  solution = sol->solution;
1858  nsol = sol->ncoil;
1859  }
1860  else {
1861  solution = m->solution;
1862  nsol = all_surfs ? m->nsol : m->surfs[0]->ntri;
1863  }
1864  for (k = 0; k < nsol; k++)
1865  pot[k] = mne_dot_vectors_40(solution[k],v0,m->nsol);
1866  return;
1867 }
1868 
1869 //=============================================================================================================
1870 
1871 int FwdBemModel::fwd_bem_pot_els(float *rd, float *Q, FwdCoilSet *els, float *pot, void *client) /* The model */
1872 /*
1873  * This version calculates the potential on all surfaces
1874  */
1875 {
1876  FwdBemModel* m = (FwdBemModel*)client;
1877  FwdBemSolution* sol = (FwdBemSolution*)els->user_data;
1878 
1879  if (!m) {
1880  printf("No BEM model specified to fwd_bem_pot_els");
1881  return FAIL;
1882  }
1883  if (!m->solution) {
1884  printf("No solution available for fwd_bem_pot_els");
1885  return FAIL;
1886  }
1887  if (!sol || sol->ncoil != els->ncoil) {
1888  printf("No appropriate electrode-specific data available in fwd_bem_pot_coils");
1889  return FAIL;
1890  }
1891  if (m->bem_method == FWD_BEM_CONSTANT_COLL)
1892  fwd_bem_pot_calc(rd,Q,m,els,FALSE,pot);
1893  else if (m->bem_method == FWD_BEM_LINEAR_COLL)
1894  fwd_bem_lin_pot_calc(rd,Q,m,els,FALSE,pot);
1895  else {
1896  printf("Unknown BEM method : %d",m->bem_method);
1897  return FAIL;
1898  }
1899  return OK;
1900 }
1901 
1902 //=============================================================================================================
1903 
1904 int FwdBemModel::fwd_bem_pot_grad_els(float *rd, float *Q, FwdCoilSet *els, float *pot, float *xgrad, float *ygrad, float *zgrad, void *client) /* The model */
1905 /*
1906  * This version calculates the potential on all surfaces
1907  */
1908 {
1909  FwdBemModel* m = (FwdBemModel*)client;
1910  FwdBemSolution* sol = (FwdBemSolution*)els->user_data;
1911 
1912  if (!m) {
1913  qCritical("No BEM model specified to fwd_bem_pot_els");
1914  return FAIL;
1915  }
1916  if (!m->solution) {
1917  qCritical("No solution available for fwd_bem_pot_els");
1918  return FAIL;
1919  }
1920  if (!sol || sol->ncoil != els->ncoil) {
1921  qCritical("No appropriate electrode-specific data available in fwd_bem_pot_coils");
1922  return FAIL;
1923  }
1924  if (m->bem_method == FWD_BEM_CONSTANT_COLL) {
1925  if (pot)
1926  fwd_bem_pot_calc(rd,Q,m,els,FALSE,pot);
1927  fwd_bem_pot_grad_calc(rd,Q,m,els,FALSE,xgrad,ygrad,zgrad);
1928  }
1929  else if (m->bem_method == FWD_BEM_LINEAR_COLL) {
1930  if (pot)
1931  fwd_bem_lin_pot_calc(rd,Q,m,els,FALSE,pot);
1932  fwd_bem_lin_pot_grad_calc(rd,Q,m,els,FALSE,xgrad,ygrad,zgrad);
1933  }
1934  else {
1935  qCritical("Unknown BEM method : %d",m->bem_method);
1936  return FAIL;
1937  }
1938  return OK;
1939 }
1940 
1941 //=============================================================================================================
1942 
1943 #define ARSINH(x) log((x) + sqrt(1.0+(x)*(x)))
1944 
1945 void FwdBemModel::calc_f(double *xx, double *yy, double *f0, double *fx, double *fy) /* The weights in the linear approximation */
1946 {
1947  double det = -xx[Y_40]*yy[X_40] + xx[Z_40]*yy[X_40] +
1948  xx[X_40]*yy[Y_40] - xx[Z_40]*yy[Y_40] - xx[X_40]*yy[Z_40] + xx[Y_40]*yy[Z_40];
1949  int k;
1950 
1951  f0[X_40] = -xx[Z_40]*yy[Y_40] + xx[Y_40]*yy[Z_40];
1952  f0[Y_40] = xx[Z_40]*yy[X_40] - xx[X_40]*yy[Z_40];
1953  f0[Z_40] = -xx[Y_40]*yy[X_40] + xx[X_40]*yy[Y_40];
1954 
1955  fx[X_40] = yy[Y_40] - yy[Z_40];
1956  fx[Y_40] = -yy[X_40] + yy[Z_40];
1957  fx[Z_40] = yy[X_40] - yy[Y_40];
1958 
1959  fy[X_40] = -xx[Y_40] + xx[Z_40];
1960  fy[Y_40] = xx[X_40] - xx[Z_40];
1961  fy[Z_40] = -xx[X_40] + xx[Y_40];
1962 
1963  for (k = 0; k < 3; k++) {
1964  f0[k] = f0[k]/det;
1965  fx[k] = fx[k]/det;
1966  fy[k] = fy[k]/det;
1967  }
1968 }
1969 
1970 //=============================================================================================================
1971 
1972 void FwdBemModel::calc_magic(double u, double z, double A, double B, double *beta, double *D)
1973 /*
1974  * Calculate Urankar's magic numbers
1975  */
1976 {
1977  double B2 = 1.0 + B*B;
1978  double ABu = A + B*u;
1979  *D = sqrt(u*u + z*z + ABu*ABu);
1980  beta[0] = ABu/sqrt(u*u + z*z);
1981  beta[1] = (A*B + B2*u)/sqrt(A*A + B2*z*z);
1982  beta[2] = (B*z*z - A*u)/(z*(*D));
1983 }
1984 
1985 //=============================================================================================================
1986 
1987 void FwdBemModel::field_integrals(float *from, MneTriangle* to, double *I1p, double *T, double *S1, double *S2, double *f0, double *fx, double *fy)
1988 {
1989  double y1[3],y2[3],y3[3];
1990  double xx[4],yy[4];
1991  double A,B,z,dx;
1992  double beta[3],I1,Tx,Ty,Txx,Tyy,Sxx,mult;
1993  double S1x,S1y,S2x,S2y;
1994  double D1,B2;
1995  int k;
1996  /*
1997  * Preliminaries...
1998  *
1999  * 1. Move origin to viewpoint...
2000  *
2001  */
2002  VEC_DIFF_40(from,to->r1,y1);
2003  VEC_DIFF_40(from,to->r2,y2);
2004  VEC_DIFF_40(from,to->r3,y3);
2005  /*
2006  * 2. Calculate local xy coordinates...
2007  */
2008  xx[0] = VEC_DOT_40(y1,to->ex);
2009  xx[1] = VEC_DOT_40(y2,to->ex);
2010  xx[2] = VEC_DOT_40(y3,to->ex);
2011  xx[3] = xx[0];
2012 
2013  yy[0] = VEC_DOT_40(y1,to->ey);
2014  yy[1] = VEC_DOT_40(y2,to->ey);
2015  yy[2] = VEC_DOT_40(y3,to->ey);
2016  yy[3] = yy[0];
2017 
2018  calc_f (xx,yy,f0,fx,fy);
2019  /*
2020  * 3. Distance of the plane from origin...
2021  */
2022  z = VEC_DOT_40(y1,to->nn);
2023  /*
2024  * Put together the line integral...
2025  * We use the convention where the local y-axis
2026  * is parallel to the last side and, therefore, dx = 0
2027  * on that side. We can thus omit the last side from this
2028  * computation in some cases.
2029  */
2030  I1 = 0.0;
2031  Tx = 0.0;
2032  Ty = 0.0;
2033  S1x = 0.0;
2034  S1y = 0.0;
2035  S2x = 0.0;
2036  S2y = 0.0;
2037  for (k = 0; k < 2; k++) {
2038  dx = xx[k+1] - xx[k];
2039  A = (yy[k]*xx[k+1] - yy[k+1]*xx[k])/dx;
2040  B = (yy[k+1]-yy[k])/dx;
2041  B2 = (1.0 + B*B);
2042  /*
2043  * Upper limit
2044  */
2045  calc_magic (xx[k+1],z,A,B,beta,&D1);
2046  I1 = I1 - xx[k+1]*ARSINH(beta[0]) - (A/sqrt(1.0+B*B))*ARSINH(beta[1])
2047  - z*atan(beta[2]);
2048  Txx = ARSINH(beta[1])/sqrt(B2);
2049  Tx = Tx + Txx;
2050  Ty = Ty + B*Txx;
2051  Sxx = (D1 - A*B*Txx)/B2;
2052  S1x = S1x + Sxx;
2053  S1y = S1y + B*Sxx;
2054  Sxx = (B*D1 + A*Txx)/B2;
2055  S2x = S2x + Sxx;
2056  /*
2057  * Lower limit
2058  */
2059  calc_magic (xx[k],z,A,B,beta,&D1);
2060  I1 = I1 + xx[k]*ARSINH(beta[0]) + (A/sqrt(1.0+B*B))*ARSINH(beta[1])
2061  + z*atan(beta[2]);
2062  Txx = ARSINH(beta[1])/sqrt(B2);
2063  Tx = Tx - Txx;
2064  Ty = Ty - B*Txx;
2065  Sxx = (D1 - A*B*Txx)/B2;
2066  S1x = S1x - Sxx;
2067  S1y = S1y - B*Sxx;
2068  Sxx = (B*D1 + A*Txx)/B2;
2069  S2x = S2x - Sxx;
2070  }
2071  /*
2072  * Handle last side (dx = 0) in a special way;
2073  */
2074  mult = 1.0/sqrt(xx[k]*xx[k]+z*z);
2075  /*
2076  * Upper...
2077  */
2078  Tyy = ARSINH(mult*yy[k+1]);
2079  Ty = Ty + Tyy;
2080  S1y = S1y + xx[k]*Tyy;
2081  /*
2082  * Lower...
2083  */
2084  Tyy = ARSINH(mult*yy[k]);
2085  Ty = Ty - Tyy;
2086  S1y = S1y - xx[k]*Tyy;
2087  /*
2088  * Set return values
2089  */
2090  *I1p = I1;
2091  T[X_40] = Tx;
2092  T[Y_40] = Ty;
2093  S1[X_40] = S1x;
2094  S1[Y_40] = S1y;
2095  S2[X_40] = S2x;
2096  S2[Y_40] = -S1x;
2097  return;
2098 }
2099 
2100 //=============================================================================================================
2101 
2102 double FwdBemModel::one_field_coeff(float *dest, float *normal, MneTriangle* tri)
2103 /*
2104  * Compute the integral over one triangle.
2105  * This looks magical but it is not.
2106  */
2107 {
2108  double *yy[4];
2109  double y1[3],y2[3],y3[3];
2110  double beta[3];
2111  double bbeta[3];
2112  double coeff[3];
2113  int j,k;
2114 
2115  yy[0] = y1;
2116  yy[1] = y2;
2117  yy[2] = y3;
2118  yy[3] = y1;
2119  VEC_DIFF_40(dest,tri->r1,y1);
2120  VEC_DIFF_40(dest,tri->r2,y2);
2121  VEC_DIFF_40(dest,tri->r3,y3);
2122  for (j = 0; j < 3; j++)
2123  beta[j] = calc_beta(yy[j],yy[j+1]);
2124  bbeta[0] = beta[2] - beta[0];
2125  bbeta[1] = beta[0] - beta[1];
2126  bbeta[2] = beta[1] - beta[2];
2127 
2128  for (j = 0; j < 3; j++)
2129  coeff[j] = 0.0;
2130  for (j = 0; j < 3; j++)
2131  for (k = 0; k < 3; k++)
2132  coeff[k] = coeff[k] + yy[j][k]*bbeta[j];
2133  return (VEC_DOT_40(coeff,normal));
2134 }
2135 
2136 //=============================================================================================================
2137 
2138 float **FwdBemModel::fwd_bem_field_coeff(FwdBemModel *m, FwdCoilSet *coils) /* Gradiometer coil positions */
2139 /*
2140  * Compute the weighting factors to obtain the magnetic field
2141  */
2142 {
2143  MneSurfaceOld* surf;
2144  MneTriangle* tri;
2145  FwdCoil* coil;
2146  FwdCoilSet* tcoils = NULL;
2147  int ntri;
2148  float **coeff = NULL;
2149  int j,k,p,s,off;
2150  double res;
2151  double mult;
2152 
2153  if (m->solution == NULL) {
2154  printf("Solution matrix missing in fwd_bem_field_coeff");
2155  return NULL;
2156  }
2157  if (m->bem_method != FWD_BEM_CONSTANT_COLL) {
2158  printf("BEM method should be constant collocation for fwd_bem_field_coeff");
2159  return NULL;
2160  }
2161  if (coils->coord_frame != FIFFV_COORD_MRI) {
2162  if (coils->coord_frame == FIFFV_COORD_HEAD) {
2163  if (!m->head_mri_t) {
2164  printf("head -> mri coordinate transform missing in fwd_bem_field_coeff");
2165  return NULL;
2166  }
2167  else {
2168  if (!coils) {
2169  qWarning("No coils to duplicate");
2170  return NULL;
2171  }
2172  /*
2173  * Make a transformed duplicate
2174  */
2175  if ((tcoils = coils->dup_coil_set(m->head_mri_t)) == NULL)
2176  return NULL;
2177  coils = tcoils;
2178  }
2179  }
2180  else {
2181  printf("Incompatible coil coordinate frame %d for fwd_bem_field_coeff",coils->coord_frame);
2182  return NULL;
2183  }
2184  }
2185  ntri = m->nsol;
2186  coeff = ALLOC_CMATRIX_40(coils->ncoil,ntri);
2187 
2188  for (s = 0, off = 0; s < m->nsurf; s++) {
2189  surf = m->surfs[s];
2190  ntri = surf->ntri;
2191  tri = surf->tris;
2192  mult = m->field_mult[s];
2193 
2194  for (k = 0; k < ntri; k++,tri++) {
2195  for (j = 0; j < coils->ncoil; j++) {
2196  coil = coils->coils[j];
2197  res = 0.0;
2198  for (p = 0; p < coil->np; p++)
2199  res = res + coil->w[p]*one_field_coeff(coil->rmag[p],coil->cosmag[p],tri);
2200  coeff[j][k+off] = mult*res;
2201  }
2202  }
2203  off = off + ntri;
2204  }
2205  delete tcoils;
2206  return coeff;
2207 }
2208 
2209 //=============================================================================================================
2210 
2211 double FwdBemModel::calc_gamma(double *rk, double *rk1)
2212 {
2213  double rkk1[3];
2214  double size;
2215  double res;
2216 
2217  VEC_DIFF_40(rk,rk1,rkk1);
2218  size = VEC_LEN_40(rkk1);
2219 
2220  res = log((VEC_LEN_40(rk1)*size + VEC_DOT_40(rk1,rkk1))/
2221  (VEC_LEN_40(rk)*size + VEC_DOT_40(rk,rkk1)))/size;
2222  return (res);
2223 }
2224 
2225 //=============================================================================================================
2226 
2227 void FwdBemModel::fwd_bem_one_lin_field_coeff_ferg(float *dest, float *dir, MneTriangle* tri, double *res) /* The results */
2228 {
2229  double c[3]; /* Component of dest vector normal to
2230  * the triangle plane */
2231  double A[3]; /* Projection of dest onto the triangle */
2232  double c1[3],c2[3],c3[3];
2233  double y1[3],y2[3],y3[3];
2234  double *yy[4],*cc[4];
2235  double rjk[3][3];
2236  double cross[3],triple,l1,l2,l3,solid,clen;
2237  double common,sum,beta,gamma;
2238  int k;
2239 
2240  yy[0] = y1; cc[0] = c1;
2241  yy[1] = y2; cc[1] = c2;
2242  yy[2] = y3; cc[2] = c3;
2243  yy[3] = y1; cc[3] = c1;
2244 
2245  VEC_DIFF_40(tri->r2,tri->r3,rjk[0]);
2246  VEC_DIFF_40(tri->r3,tri->r1,rjk[1]);
2247  VEC_DIFF_40(tri->r1,tri->r2,rjk[2]);
2248 
2249  for (k = 0; k < 3; k++) {
2250  y1[k] = tri->r1[k] - dest[k];
2251  y2[k] = tri->r2[k] - dest[k];
2252  y3[k] = tri->r3[k] - dest[k];
2253  }
2254  clen = VEC_DOT_40(y1,tri->nn);
2255  for (k = 0; k < 3; k++) {
2256  c[k] = clen*tri->nn[k];
2257  A[k] = dest[k] + c[k];
2258  c1[k] = tri->r1[k] - A[k];
2259  c2[k] = tri->r2[k] - A[k];
2260  c3[k] = tri->r3[k] - A[k];
2261  }
2262  /*
2263  * beta and gamma...
2264  */
2265  for (sum = 0.0, k = 0; k < 3; k++) {
2266  CROSS_PRODUCT_40(cc[k],cc[k+1],cross);
2267  beta = VEC_DOT_40(cross,tri->nn);
2268  gamma = calc_gamma (yy[k],yy[k+1]);
2269  sum = sum + beta*gamma;
2270  }
2271  /*
2272  * Solid angle...
2273  */
2274  CROSS_PRODUCT_40(y1,y2,cross);
2275  triple = VEC_DOT_40(cross,y3);
2276 
2277  l1 = VEC_LEN_40(y1);
2278  l2 = VEC_LEN_40(y2);
2279  l3 = VEC_LEN_40(y3);
2280  solid = 2.0*atan2(triple,
2281  (l1*l2*l3+
2282  VEC_DOT_40(y1,y2)*l3+
2283  VEC_DOT_40(y1,y3)*l2+
2284  VEC_DOT_40(y2,y3)*l1));
2285  /*
2286  * Now we are ready to assemble it all together
2287  */
2288  common = (sum-clen*solid)/(2.0*tri->area);
2289  for (k = 0; k < 3; k++)
2290  res[k] = -VEC_DOT_40(rjk[k],dir)*common;
2291  return;
2292 }
2293 
2294 //=============================================================================================================
2295 
2296 void FwdBemModel::fwd_bem_one_lin_field_coeff_uran(float *dest, float *dir, MneTriangle* tri, double *res) /* The results */
2297 {
2298  double I1,T[2],S1[2],S2[2];
2299  double f0[3],fx[3],fy[3];
2300  double res_x,res_y;
2301  double x_fac,y_fac;
2302  int k;
2303  double len;
2304  /*
2305  * Compute the component integrals
2306  */
2307  field_integrals (dest,tri,&I1,T,S1,S2,f0,fx,fy);
2308  /*
2309  * Compute the coefficient for each node...
2310  */
2311  len = VEC_LEN_40(dir);
2312  dir[X_40] = dir[X_40]/len;
2313  dir[Y_40] = dir[Y_40]/len;
2314  dir[Z_40] = dir[Z_40]/len;
2315 
2316  x_fac = -VEC_DOT_40(dir,tri->ex);
2317  y_fac = -VEC_DOT_40(dir,tri->ey);
2318  for (k = 0; k < 3; k++) {
2319  res_x = f0[k]*T[X_40] + fx[k]*S1[X_40] + fy[k]*S2[X_40] + fy[k]*I1;
2320  res_y = f0[k]*T[Y_40] + fx[k]*S1[Y_40] + fy[k]*S2[Y_40] - fx[k]*I1;
2321  res[k] = x_fac*res_x + y_fac*res_y;
2322  }
2323  return;
2324 }
2325 
2326 //=============================================================================================================
2327 
2328 void FwdBemModel::fwd_bem_one_lin_field_coeff_simple(float *dest, float *normal, MneTriangle* source, double *res) /* The result for each triangle node */
2329 /*
2330  * Simple version...
2331  */
2332 {
2333  float diff[3];
2334  float vec_result[3];
2335  float dl;
2336  int k;
2337  float *rr[3];
2338 
2339  rr[0] = source->r1;
2340  rr[1] = source->r2;
2341  rr[2] = source->r3;
2342 
2343  for (k = 0; k < 3; k++) {
2344  VEC_DIFF_40(rr[k],dest,diff);
2345  dl = VEC_DOT_40(diff,diff);
2346  CROSS_PRODUCT_40(diff,source->nn,vec_result);
2347  res[k] = source->area*VEC_DOT_40(vec_result,normal)/(3.0*dl*sqrt(dl));
2348  }
2349  return;
2350 }
2351 
2352 //=============================================================================================================
2353 
2354 float **FwdBemModel::fwd_bem_lin_field_coeff(FwdBemModel *m, FwdCoilSet *coils, int method) /* Which integration formula to use */
2355 /*
2356  * Compute the weighting factors to obtain the magnetic field
2357  * in the linear potential approximation
2358  */
2359 {
2360  MneSurfaceOld* surf;
2361  MneTriangle* tri;
2362  FwdCoil* coil;
2363  FwdCoilSet* tcoils = NULL;
2364  int ntri;
2365  float **coeff = NULL;
2366  int j,k,p,pp,off,s;
2367  double res[3],one[3];
2368  float mult;
2369  linFieldIntFunc func;
2370 
2371  if (m->solution == NULL) {
2372  printf("Solution matrix missing in fwd_bem_lin_field_coeff");
2373  return NULL;
2374  }
2375  if (m->bem_method != FWD_BEM_LINEAR_COLL) {
2376  printf("BEM method should be linear collocation for fwd_bem_lin_field_coeff");
2377  return NULL;
2378  }
2379  if (coils->coord_frame != FIFFV_COORD_MRI) {
2380  if (coils->coord_frame == FIFFV_COORD_HEAD) {
2381  if (!m->head_mri_t) {
2382  printf("head -> mri coordinate transform missing in fwd_bem_lin_field_coeff");
2383  return NULL;
2384  }
2385  else {
2386  if (!coils) {
2387  qWarning("No coils to duplicate");
2388  return NULL;
2389  }
2390  /*
2391  * Make a transformed duplicate
2392  */
2393  if ((tcoils = coils->dup_coil_set(m->head_mri_t)) == NULL)
2394  return NULL;
2395  coils = tcoils;
2396  }
2397  }
2398  else {
2399  printf("Incompatible coil coordinate frame %d for fwd_bem_field_coeff",coils->coord_frame);
2400  return NULL;
2401  }
2402  }
2403  if (method == FWD_BEM_LIN_FIELD_FERGUSON)
2404  func = fwd_bem_one_lin_field_coeff_ferg;
2405  else if (method == FWD_BEM_LIN_FIELD_URANKAR)
2406  func = fwd_bem_one_lin_field_coeff_uran;
2407  else
2408  func = fwd_bem_one_lin_field_coeff_simple;
2409 
2410  coeff = ALLOC_CMATRIX_40(coils->ncoil,m->nsol);
2411  for (k = 0; k < m->nsol; k++)
2412  for (j = 0; j < coils->ncoil; j++)
2413  coeff[j][k] = 0.0;
2414  /*
2415  * Process each of the surfaces
2416  */
2417  for (s = 0, off = 0; s < m->nsurf; s++) {
2418  surf = m->surfs[s];
2419  ntri = surf->ntri;
2420  tri = surf->tris;
2421  mult = m->field_mult[s];
2422 
2423  for (k = 0; k < ntri; k++,tri++) {
2424  for (j = 0; j < coils->ncoil; j++) {
2425  coil = coils->coils[j];
2426  for (pp = 0; pp < 3; pp++)
2427  res[pp] = 0;
2428  /*
2429  * Accumulate the coefficients for each triangle node...
2430  */
2431  for (p = 0; p < coil->np; p++) {
2432  func(coil->rmag[p],coil->cosmag[p],tri,one);
2433  for (pp = 0; pp < 3; pp++)
2434  res[pp] = res[pp] + coil->w[p]*one[pp];
2435  }
2436  /*
2437  * Add these to the corresponding coefficient matrix
2438  * elements...
2439  */
2440  for (pp = 0; pp < 3; pp++)
2441  coeff[j][tri->vert[pp]+off] = coeff[j][tri->vert[pp]+off] + mult*res[pp];
2442  }
2443  }
2444  off = off + surf->np;
2445  }
2446  /*
2447  * Discard the duplicate
2448  */
2449  delete tcoils;
2450  return (coeff);
2451 }
2452 
2453 //=============================================================================================================
2454 
2455 int FwdBemModel::fwd_bem_specify_coils(FwdBemModel *m, FwdCoilSet *coils)
2456 /*
2457  * Set up for computing the solution at a set of coils
2458  */
2459 {
2460  float **sol = NULL;
2461  FwdBemSolution* csol;
2462 
2463  if (!m) {
2464  printf("Model missing in fwd_bem_specify_coils");
2465  goto bad;
2466  }
2467  if (!m->solution) {
2468  printf("Solution not computed in fwd_bem_specify_coils");
2469  goto bad;
2470  }
2471  if(coils)
2472  coils->fwd_free_coil_set_user_data();
2473  if (!coils || coils->ncoil == 0)
2474  return OK;
2475  if (m->bem_method == FWD_BEM_CONSTANT_COLL)
2476  sol = fwd_bem_field_coeff(m,coils);
2477  else if (m->bem_method == FWD_BEM_LINEAR_COLL)
2478  sol = fwd_bem_lin_field_coeff(m,coils,FWD_BEM_LIN_FIELD_SIMPLE);
2479  else {
2480  printf("Unknown BEM method in fwd_bem_specify_coils : %d",m->bem_method);
2481  goto bad;
2482  }
2483  coils->user_data = csol = new FwdBemSolution();
2484  coils->user_data_free = FwdBemSolution::fwd_bem_free_coil_solution;
2485 
2486  csol->ncoil = coils->ncoil;
2487  csol->np = m->nsol;
2488  csol->solution = mne_mat_mat_mult_40(sol,
2489  m->solution,
2490  coils->ncoil,
2491  m->nsol,
2492  m->nsol);//TODO: Suspicion, that this is slow - use Eigen
2493 
2494  FREE_CMATRIX_40(sol);
2495  return OK;
2496 
2497 bad : {
2498  FREE_CMATRIX_40(sol);
2499  return FAIL;
2500 
2501  }
2502 }
2503 
2504 //=============================================================================================================
2505 
2506 void FwdBemModel::fwd_bem_lin_field_calc(float *rd, float *Q, FwdCoilSet *coils, FwdBemModel *m, float *B)
2507 /*
2508  * Calculate the magnetic field in a set of coils
2509  */
2510 {
2511  float *v0;
2512  int s,k,p,np;
2513  FwdCoil* coil;
2514  float mult;
2515  float **rr;
2516  float my_rd[3],my_Q[3];
2517  FwdBemSolution* sol = (FwdBemSolution*)coils->user_data;
2518  /*
2519  * Infinite-medium potentials
2520  */
2521  if (!m->v0)
2522  m->v0 = MALLOC_40(m->nsol,float);
2523  v0 = m->v0;
2524  /*
2525  * The dipole location and orientation must be transformed
2526  */
2527  VEC_COPY_40(my_rd,rd);
2528  VEC_COPY_40(my_Q,Q);
2529  if (m->head_mri_t) {
2530  FiffCoordTransOld::fiff_coord_trans(my_rd,m->head_mri_t,FIFFV_MOVE);
2531  FiffCoordTransOld::fiff_coord_trans(my_Q,m->head_mri_t,FIFFV_NO_MOVE);
2532  }
2533  /*
2534  * Compute the inifinite-medium potentials at the vertices
2535  */
2536  for (s = 0, p = 0; s < m->nsurf; s++) {
2537  np = m->surfs[s]->np;
2538  rr = m->surfs[s]->rr;
2539  mult = m->source_mult[s];
2540  for (k = 0; k < np; k++)
2541  v0[p++] = mult*fwd_bem_inf_pot(my_rd,my_Q,rr[k]);
2542  }
2543  /*
2544  * Primary current contribution
2545  * (can be calculated in the coil/dipole coordinates)
2546  */
2547  for (k = 0; k < coils->ncoil; k++) {
2548  coil = coils->coils[k];
2549  B[k] = 0.0;
2550  for (p = 0; p < coil->np; p++)
2551  B[k] = B[k] + coil->w[p]*fwd_bem_inf_field(rd,Q,coil->rmag[p],coil->cosmag[p]);
2552  }
2553  /*
2554  * Volume current contribution
2555  */
2556  for (k = 0; k < coils->ncoil; k++)
2557  B[k] = B[k] + mne_dot_vectors_40(sol->solution[k],v0,m->nsol);
2558  /*
2559  * Scale correctly
2560  */
2561  for (k = 0; k < coils->ncoil; k++)
2562  B[k] = MAG_FACTOR*B[k];
2563  return;
2564 }
2565 
2566 //=============================================================================================================
2567 
2568 void FwdBemModel::fwd_bem_field_calc(float *rd, float *Q, FwdCoilSet *coils, FwdBemModel *m, float *B)
2569 /*
2570  * Calculate the magnetic field in a set of coils
2571  */
2572 {
2573  float *v0;
2574  int s,k,p,ntri;
2575  FwdCoil* coil;
2576  MneTriangle* tri;
2577  float mult;
2578  float my_rd[3],my_Q[3];
2579  FwdBemSolution* sol = (FwdBemSolution*)coils->user_data;
2580  /*
2581  * Infinite-medium potentials
2582  */
2583  if (!m->v0)
2584  m->v0 = MALLOC_40(m->nsol,float);
2585  v0 = m->v0;
2586  /*
2587  * The dipole location and orientation must be transformed
2588  */
2589  VEC_COPY_40(my_rd,rd);
2590  VEC_COPY_40(my_Q,Q);
2591  if (m->head_mri_t) {
2592  FiffCoordTransOld::fiff_coord_trans(my_rd,m->head_mri_t,FIFFV_MOVE);
2593  FiffCoordTransOld::fiff_coord_trans(my_Q,m->head_mri_t,FIFFV_NO_MOVE);
2594  }
2595  /*
2596  * Compute the inifinite-medium potentials at the centers of the triangles
2597  */
2598  for (s = 0, p = 0; s < m->nsurf; s++) {
2599  ntri = m->surfs[s]->ntri;
2600  tri = m->surfs[s]->tris;
2601  mult = m->source_mult[s];
2602  for (k = 0; k < ntri; k++, tri++)
2603  v0[p++] = mult*fwd_bem_inf_pot(my_rd,my_Q,tri->cent);
2604  }
2605  /*
2606  * Primary current contribution
2607  * (can be calculated in the coil/dipole coordinates)
2608  */
2609  for (k = 0; k < coils->ncoil; k++) {
2610  coil = coils->coils[k];
2611  B[k] = 0.0;
2612  for (p = 0; p < coil->np; p++)
2613  B[k] = B[k] + coil->w[p]*fwd_bem_inf_field(rd,Q,coil->rmag[p],coil->cosmag[p]);
2614  }
2615  /*
2616  * Volume current contribution
2617  */
2618  for (k = 0; k < coils->ncoil; k++)
2619  B[k] = B[k] + mne_dot_vectors_40(sol->solution[k],v0,m->nsol);
2620  /*
2621  * Scale correctly
2622  */
2623  for (k = 0; k < coils->ncoil; k++)
2624  B[k] = MAG_FACTOR*B[k];
2625  return;
2626 }
2627 
2628 //=============================================================================================================
2629 
2630 void FwdBemModel::fwd_bem_field_grad_calc(float *rd, float *Q, FwdCoilSet* coils, FwdBemModel* m, float *xgrad, float *ygrad, float *zgrad)
2631 /*
2632  * Calculate the magnetic field in a set of coils
2633  */
2634 {
2635  FwdBemSolution* sol = (FwdBemSolution*)coils->user_data;
2636  float *v0;
2637  int s,k,p,ntri,pp;
2638  FwdCoil* coil;
2639  MneTriangle* tri;
2640  float mult;
2641  float *grads[3],ee[3],mri_ee[3],mri_rd[3],mri_Q[3];
2642  float *grad;
2643 
2644  grads[0] = xgrad;
2645  grads[1] = ygrad;
2646  grads[2] = zgrad;
2647  /*
2648  * Infinite-medium potentials
2649  */
2650  if (!m->v0)
2651  m->v0 = MALLOC_40(m->nsol,float);
2652  v0 = m->v0;
2653  /*
2654  * The dipole location and orientation must be transformed
2655  */
2656  VEC_COPY_40(mri_rd,rd);
2657  VEC_COPY_40(mri_Q,Q);
2658  if (m->head_mri_t) {
2659  FiffCoordTransOld::fiff_coord_trans(mri_rd,m->head_mri_t,FIFFV_MOVE);
2660  FiffCoordTransOld::fiff_coord_trans(mri_Q,m->head_mri_t,FIFFV_NO_MOVE);
2661  }
2662  for (pp = 0; pp < 3; pp++) {
2663  grad = grads[pp];
2664  /*
2665  * Select the correct gradient component
2666  */
2667  for (p = 0; p < 3; p++)
2668  ee[p] = p == pp ? 1.0 : 0.0;
2669  VEC_COPY_40(mri_ee,ee);
2670  if (m->head_mri_t)
2671  FiffCoordTransOld::fiff_coord_trans(mri_ee,m->head_mri_t,FIFFV_NO_MOVE);
2672  /*
2673  * Compute the inifinite-medium potential derivatives at the centers of the triangles
2674  */
2675  for (s = 0, p = 0; s < m->nsurf; s++) {
2676  ntri = m->surfs[s]->ntri;
2677  tri = m->surfs[s]->tris;
2678  mult = m->source_mult[s];
2679  for (k = 0; k < ntri; k++, tri++) {
2680  v0[p++] = mult*fwd_bem_inf_pot_der(mri_rd,mri_Q,tri->cent,mri_ee);
2681  }
2682  }
2683  /*
2684  * Primary current contribution
2685  * (can be calculated in the coil/dipole coordinates)
2686  */
2687  for (k = 0; k < coils->ncoil; k++) {
2688  coil = coils->coils[k];
2689  grad[k] = 0.0;
2690  for (p = 0; p < coil->np; p++)
2691  grad[k] = grad[k] + coil->w[p]*fwd_bem_inf_field_der(rd,Q,coil->rmag[p],coil->cosmag[p],ee);
2692  }
2693  /*
2694  * Volume current contribution
2695  */
2696  for (k = 0; k < coils->ncoil; k++)
2697  grad[k] = grad[k] + mne_dot_vectors_40(sol->solution[k],v0,m->nsol);
2698  /*
2699  * Scale correctly
2700  */
2701  for (k = 0; k < coils->ncoil; k++)
2702  grad[k] = MAG_FACTOR*grad[k];
2703  }
2704  return;
2705 }
2706 
2707 //=============================================================================================================
2708 
2709 float FwdBemModel::fwd_bem_inf_field_der(float *rd, float *Q, float *rp, float *dir, float *comp) /* Which gradient component */
2710 /*
2711  * Derivative of the infinite-medium magnetic field with respect to
2712  * one of the dipole position coordinates (without \mu_0/4\pi)
2713  */
2714 {
2715  float diff[3],diff2,diff3,diff5,cross[3],crossn[3],res;
2716 
2717  VEC_DIFF_40(rd,rp,diff);
2718  diff2 = VEC_DOT_40(diff,diff);
2719  diff3 = sqrt(diff2)*diff2;
2720  diff5 = diff3*diff2;
2721  CROSS_PRODUCT_40(Q,diff,cross);
2722  CROSS_PRODUCT_40(dir,Q,crossn);
2723 
2724  res = 3*VEC_DOT_40(cross,dir)*VEC_DOT_40(comp,diff)/diff5 - VEC_DOT_40(comp,crossn)/diff3;
2725  return res;
2726 }
2727 
2728 //=============================================================================================================
2729 
2730 float FwdBemModel::fwd_bem_inf_pot_der(float *rd, float *Q, float *rp, float *comp) /* Which gradient component */
2731 /*
2732  * Derivative of the infinite-medium potential with respect to one of
2733  * the dipole position coordinates
2734  */
2735 {
2736  float diff[3];
2737  float diff2,diff5,diff3;
2738  float res;
2739 
2740  VEC_DIFF_40(rd,rp,diff);
2741  diff2 = VEC_DOT_40(diff,diff);
2742  diff3 = sqrt(diff2)*diff2;
2743  diff5 = diff3*diff2;
2744 
2745  res = 3*VEC_DOT_40(Q,diff)*VEC_DOT_40(comp,diff)/diff5 - VEC_DOT_40(comp,Q)/diff3;
2746  return (res/(4.0*M_PI));
2747 }
2748 
2749 //=============================================================================================================
2750 
2751 void FwdBemModel::fwd_bem_lin_field_grad_calc(float *rd, float *Q, FwdCoilSet *coils, FwdBemModel *m, float *xgrad, float *ygrad, float *zgrad)
2752 /*
2753  * Calculate the gradient with respect to dipole position coordinates in a set of coils
2754  */
2755 {
2756  FwdBemSolution* sol = (FwdBemSolution*)coils->user_data;
2757 
2758  float *v0;
2759  int s,k,p,np,pp;
2760  FwdCoil *coil;
2761  float mult;
2762  float **rr,ee[3],mri_ee[3],mri_rd[3],mri_Q[3];
2763  float *grads[3];
2764  float *grad;
2765 
2766  grads[0] = xgrad;
2767  grads[1] = ygrad;
2768  grads[2] = zgrad;
2769  /*
2770  * Space for infinite-medium potentials
2771  */
2772  if (!m->v0)
2773  m->v0 = MALLOC_40(m->nsol,float);
2774  v0 = m->v0;
2775  /*
2776  * The dipole location and orientation must be transformed
2777  */
2778  VEC_COPY_40(mri_rd,rd);
2779  VEC_COPY_40(mri_Q,Q);
2780  if (m->head_mri_t) {
2781  FiffCoordTransOld::fiff_coord_trans(mri_rd,m->head_mri_t,FIFFV_MOVE);
2782  FiffCoordTransOld::fiff_coord_trans(mri_Q,m->head_mri_t,FIFFV_NO_MOVE);
2783  }
2784  for (pp = 0; pp < 3; pp++) {
2785  grad = grads[pp];
2786  /*
2787  * Select the correct gradient component
2788  */
2789  for (p = 0; p < 3; p++)
2790  ee[p] = p == pp ? 1.0 : 0.0;
2791  VEC_COPY_40(mri_ee,ee);
2792  if (m->head_mri_t)
2793  FiffCoordTransOld::fiff_coord_trans(mri_ee,m->head_mri_t,FIFFV_NO_MOVE);
2794  /*
2795  * Compute the inifinite-medium potentials at the vertices
2796  */
2797  for (s = 0, p = 0; s < m->nsurf; s++) {
2798  np = m->surfs[s]->np;
2799  rr = m->surfs[s]->rr;
2800  mult = m->source_mult[s];
2801 
2802  for (k = 0; k < np; k++)
2803  v0[p++] = mult*fwd_bem_inf_pot_der(mri_rd,mri_Q,rr[k],mri_ee);
2804  }
2805  /*
2806  * Primary current contribution
2807  * (can be calculated in the coil/dipole coordinates)
2808  */
2809  for (k = 0; k < coils->ncoil; k++) {
2810  coil = coils->coils[k];
2811  grad[k] = 0.0;
2812  for (p = 0; p < coil->np; p++)
2813  grad[k] = grad[k] + coil->w[p]*fwd_bem_inf_field_der(rd,Q,coil->rmag[p],coil->cosmag[p],ee);
2814  }
2815  /*
2816  * Volume current contribution
2817  */
2818  for (k = 0; k < coils->ncoil; k++)
2819  grad[k] = grad[k] + mne_dot_vectors_40(sol->solution[k],v0,m->nsol);
2820  /*
2821  * Scale correctly
2822  */
2823  for (k = 0; k < coils->ncoil; k++)
2824  grad[k] = MAG_FACTOR*grad[k];
2825  }
2826  return;
2827 }
2828 
2829 //=============================================================================================================
2830 
2831 int FwdBemModel::fwd_bem_field(float *rd, float *Q, FwdCoilSet *coils, float *B, void *client) /* The model */
2832 /*
2833  * This version calculates the magnetic field in a set of coils
2834  * Call fwd_bem_specify_coils first to establish the coil-specific
2835  * solution matrix
2836  */
2837 {
2838  FwdBemModel* m = (FwdBemModel*)client;
2839  FwdBemSolution* sol = (FwdBemSolution*)coils->user_data;
2840 
2841  if (!m) {
2842  printf("No BEM model specified to fwd_bem_field");
2843  return FAIL;
2844  }
2845  if (!sol || !sol->solution || sol->ncoil != coils->ncoil) {
2846  printf("No appropriate coil-specific data available in fwd_bem_field");
2847  return FAIL;
2848  }
2849  if (m->bem_method == FWD_BEM_CONSTANT_COLL)
2850  fwd_bem_field_calc(rd,Q,coils,m,B);
2851  else if (m->bem_method == FWD_BEM_LINEAR_COLL)
2852  fwd_bem_lin_field_calc(rd,Q,coils,m,B);
2853  else {
2854  printf("Unknown BEM method : %d",m->bem_method);
2855  return FAIL;
2856  }
2857  return OK;
2858 }
2859 
2860 //=============================================================================================================
2861 
2862 int FwdBemModel::fwd_bem_field_grad(float *rd,
2863  float Q[],
2864  FwdCoilSet *coils,
2865  float Bval[],
2866  float xgrad[],
2867  float ygrad[],
2868  float zgrad[],
2869  void *client) /* Client data to be passed to some foward modelling routines */
2870 {
2871  FwdBemModel* m = (FwdBemModel*)client;
2872  FwdBemSolution* sol = (FwdBemSolution*)coils->user_data;
2873 
2874  if (!m) {
2875  qCritical("No BEM model specified to fwd_bem_field");
2876  return FAIL;
2877  }
2878 
2879  if (!sol || !sol->solution || sol->ncoil != coils->ncoil) {
2880  qCritical("No appropriate coil-specific data available in fwd_bem_field");
2881  return FAIL;
2882  }
2883 
2884  if (m->bem_method == FWD_BEM_CONSTANT_COLL) {
2885  if (Bval) {
2886  fwd_bem_field_calc(rd,
2887  Q,
2888  coils,
2889  m,
2890  Bval);
2891  }
2892 
2893  fwd_bem_field_grad_calc(rd,
2894  Q,
2895  coils,
2896  m,
2897  xgrad,
2898  ygrad,
2899  zgrad);
2900  } else if (m->bem_method == FWD_BEM_LINEAR_COLL) {
2901  if (Bval) {
2902  fwd_bem_lin_field_calc(rd,
2903  Q,
2904  coils,
2905  m,
2906  Bval);
2907  }
2908 
2909  fwd_bem_lin_field_grad_calc(rd,
2910  Q,
2911  coils,
2912  m,
2913  xgrad,
2914  ygrad,
2915  zgrad);
2916  } else {
2917  qCritical("Unknown BEM method : %d",m->bem_method);
2918  return FAIL;
2919  }
2920 
2921  return OK;
2922 }
2923 
2924 //=============================================================================================================
2925 
2926 void *FwdBemModel::meg_eeg_fwd_one_source_space(void *arg)
2927 /*
2928  * Compute the MEG or EEG forward solution for one source space
2929  * and possibly for only one source component
2930  */
2931 {
2932  FwdThreadArg* a = (FwdThreadArg*)arg;
2933  MneSourceSpaceOld* s = a->s;
2934  int j,p,q;
2935  float *xyz[3];
2936 
2937  p = a->off;
2938  q = 3*a->off;
2939  if (a->fixed_ori) { /* The normal source component only */
2940  if (a->field_pot_grad && a->res_grad) { /* Gradient requested? */
2941  for (j = 0; j < s->np; j++) {
2942  if (s->inuse[j]) {
2943  if (a->field_pot_grad(s->rr[j],
2944  s->nn[j],
2945  a->coils_els,
2946  a->res[p],
2947  a->res_grad[q],
2948  a->res_grad[q+1],
2949  a->res_grad[q+2],
2950  a->client) != OK) {
2951  goto bad;
2952  }
2953  q = q + 3;
2954  p++;
2955  }
2956  }
2957  } else {
2958  for (j = 0; j < s->np; j++)
2959  if (s->inuse[j])
2960  if (a->field_pot(s->rr[j],
2961  s->nn[j],
2962  a->coils_els,
2963  a->res[p++],
2964  a->client) != OK)
2965  goto bad;
2966  }
2967  }
2968  else { /* All source components */
2969  if (a->field_pot_grad && a->res_grad) { /* Gradient requested? */
2970  for (j = 0; j < s->np; j++) {
2971  if (s->inuse[j]) {
2972  if (a->comp < 0) { /* Compute all components */
2973  if (a->field_pot_grad(s->rr[j],
2974  Qx,
2975  a->coils_els,
2976  a->res[p],
2977  a->res_grad[q],
2978  a->res_grad[q+1],
2979  a->res_grad[q+2],
2980  a->client) != OK)
2981  goto bad;
2982  q = q + 3; p++;
2983  if (a->field_pot_grad(s->rr[j],
2984  Qy,
2985  a->coils_els,
2986  a->res[p],
2987  a->res_grad[q],
2988  a->res_grad[q+1],
2989  a->res_grad[q+2],
2990  a->client) != OK)
2991  goto bad;
2992  q = q + 3; p++;
2993  if (a->field_pot_grad(s->rr[j],
2994  Qz,
2995  a->coils_els,
2996  a->res[p],
2997  a->res_grad[q],
2998  a->res_grad[q+1],
2999  a->res_grad[q+2],
3000  a->client) != OK)
3001  goto bad;
3002  q = q + 3; p++;
3003  }
3004  else if (a->comp == 0) { /* Compute x component */
3005  if (a->field_pot_grad(s->rr[j],
3006  Qx,
3007  a->coils_els,
3008  a->res[p],
3009  a->res_grad[q],
3010  a->res_grad[q+1],
3011  a->res_grad[q+2],
3012  a->client) != OK)
3013  goto bad;
3014  q = q + 3; p++;
3015  q = q + 3; p++;
3016  q = q + 3; p++;
3017  }
3018  else if (a->comp == 1) { /* Compute y component */
3019  q = q + 3; p++;
3020  if (a->field_pot_grad(s->rr[j],
3021  Qy,
3022  a->coils_els,
3023  a->res[p],
3024  a->res_grad[q],
3025  a->res_grad[q+1],
3026  a->res_grad[q+2],
3027  a->client) != OK)
3028  goto bad;
3029  q = q + 3; p++;
3030  q = q + 3; p++;
3031  }
3032  else if (a->comp == 2) { /* Compute z component */
3033  q = q + 3; p++;
3034  q = q + 3; p++;
3035  if (a->field_pot_grad(s->rr[j],
3036  Qz,
3037  a->coils_els,
3038  a->res[p],
3039  a->res_grad[q],
3040  a->res_grad[q+1],
3041  a->res_grad[q+2],
3042  a->client) != OK)
3043  goto bad;
3044  q = q + 3; p++;
3045  }
3046  }
3047  }
3048  }
3049  else {
3050  for (j = 0; j < s->np; j++) {
3051  if (s->inuse[j]) {
3052  if (a->vec_field_pot) {
3053  xyz[0] = a->res[p++];
3054  xyz[1] = a->res[p++];
3055  xyz[2] = a->res[p++];
3056  if (a->vec_field_pot(s->rr[j],a->coils_els,xyz,a->client) != OK)
3057  goto bad;
3058  }
3059  else {
3060  if (a->comp < 0) { /* Compute all components here */
3061  if (a->field_pot(s->rr[j],Qx,a->coils_els,a->res[p++],a->client) != OK)
3062  goto bad;
3063  if (a->field_pot(s->rr[j],Qy,a->coils_els,a->res[p++],a->client) != OK)
3064  goto bad;
3065  if (a->field_pot(s->rr[j],Qz,a->coils_els,a->res[p++],a->client) != OK)
3066  goto bad;
3067  }
3068  else if (a->comp == 0) { /* Compute x component */
3069  if (a->field_pot(s->rr[j],Qx,a->coils_els,a->res[p++],a->client) != OK)
3070  goto bad;
3071  p++; p++;
3072  }
3073  else if (a->comp == 1) { /* Compute y component */
3074  p++;
3075  if (a->field_pot(s->rr[j],Qy,a->coils_els,a->res[p++],a->client) != OK)
3076  goto bad;
3077  p++;
3078  }
3079  else if (a->comp == 2) { /* Compute z component */
3080  p++; p++;
3081  if (a->field_pot(s->rr[j],Qz,a->coils_els,a->res[p++],a->client) != OK)
3082  goto bad;
3083  }
3084  }
3085  }
3086  }
3087  }
3088  }
3089  a->stat = OK;
3090  return NULL;
3091 
3092 bad : {
3093  a->stat = FAIL;
3094  return NULL;
3095  }
3096 }
3097 
3098 //=============================================================================================================
3099 
3101  int nspace,
3102  FwdCoilSet *coils,
3103  FwdCoilSet *comp_coils,
3104  MneCTFCompDataSet *comp_data,
3105  bool fixed_ori,
3106  FwdBemModel *bem_model,
3107  Vector3f *r0,
3108  bool use_threads,
3109  FiffNamedMatrix& resp,
3110  FiffNamedMatrix& resp_grad,
3111  bool bDoGrad)
3112 /*
3113  * Compute the MEG forward solution
3114  * Use either the sphere model or BEM in the calculations
3115  */
3116 {
3117  float **res = NULL; /* The forward solution matrix */
3118  float **res_grad = NULL; /* The gradient with respect to the dipole position */
3119  MatrixXd matRes;
3120  MatrixXd matResGrad;
3121  int nrow = 0;
3122  FwdCompData *comp = NULL;
3123  fwdFieldFunc field; /* Computes the field for one dipole orientation */
3124  fwdVecFieldFunc vec_field; /* Computes the field for all dipole orientations */
3125  fwdFieldGradFunc field_grad; /* Computes the field and gradient with respect to dipole position
3126  * for one dipole orientation */
3127  int nmeg = coils->ncoil;/* Number of channels */
3128  int nsource; /* Total number of sources */
3129  int k,p,q,off;
3130  QStringList names; /* Channel names */
3131  void *client;
3132  FwdThreadArg* one_arg = NULL;
3133  int nproc = QThread::idealThreadCount();
3134  QStringList emptyList;
3135 
3136  if (bem_model) {
3137  /*
3138  * Use the new compensated field computation
3139  * It works the same way independent of whether or not the compensation is in effect
3140  */
3141 #ifdef TEST
3142  printf("Using differences.\n");
3143  comp = FwdCompData::fwd_make_comp_data(comp_data,
3144  coils,comp_coils,
3145  FwdBemModel::fwd_bem_field,
3146  NULL,
3147  my_bem_field_grad,
3148  bem_model,
3149  NULL);
3150 #else
3151  comp = FwdCompData::fwd_make_comp_data(comp_data,
3152  coils,
3153  comp_coils,
3154  FwdBemModel::fwd_bem_field,
3155  NULL,
3156  FwdBemModel::fwd_bem_field_grad,
3157  bem_model,
3158  NULL);
3159 #endif
3160  if (!comp)
3161  goto bad;
3162  /*
3163  * Field computation matrices...
3164  */
3165  qDebug() << "!!!TODO Speed the following with Eigen up!";
3166  printf("Composing the field computation matrix...");
3167  if (fwd_bem_specify_coils(bem_model,coils) == FAIL)
3168  goto bad;
3169  printf("[done]\n");
3170 
3171  if (comp->set && comp->set->current) { /* Test just to specify confusing output */
3172  printf("Composing the field computation matrix (compensation coils)...");
3173  if (fwd_bem_specify_coils(bem_model,comp->comp_coils) == FAIL)
3174  goto bad;
3175  printf("[done]\n");
3176  }
3177  field = FwdCompData::fwd_comp_field;
3178  vec_field = NULL;
3179  field_grad = FwdCompData::fwd_comp_field_grad;
3180  client = comp;
3181  }
3182  else {
3183  /*
3184  * Use the new compensated field computation
3185  * It works the same way independent of whether or not the compensation is in effect
3186  */
3187 #ifdef TEST
3188  printf("Using differences.\n");
3189  comp = FwdCompData::fwd_make_comp_data(comp_data,coils,comp_coils,
3190  fwd_sphere_field,
3191  fwd_sphere_field_vec,
3192  my_sphere_field_grad,
3193  r0,NULL);
3194 #else
3195  comp = FwdCompData::fwd_make_comp_data(comp_data,coils,comp_coils,
3196  fwd_sphere_field,
3197  fwd_sphere_field_vec,
3198  fwd_sphere_field_grad,
3199  r0,NULL);
3200 #endif
3201  if (!comp)
3202  goto bad;
3203  field = FwdCompData::fwd_comp_field;
3204  vec_field = FwdCompData::fwd_comp_field_vec;
3205  field_grad = FwdCompData::fwd_comp_field_grad;
3206  client = comp;
3207  }
3208  /*
3209  * Count the sources
3210  */
3211  for (k = 0, nsource = 0; k < nspace; k++)
3212  nsource += spaces[k]->nuse;
3213  /*
3214  * Allocate space for the solution
3215  */
3216  if (fixed_ori)
3217  res = ALLOC_CMATRIX_40(nsource,nmeg);
3218  else
3219  res = ALLOC_CMATRIX_40(3*nsource,nmeg);
3220  if (bDoGrad) {
3221  if (fixed_ori)
3222  res_grad = ALLOC_CMATRIX_40(3*nsource,nmeg);
3223  else
3224  res_grad = ALLOC_CMATRIX_40(3*3*nsource,nmeg);
3225  }
3226  /*
3227  * Set up the argument for the field computation
3228  */
3229  one_arg = new FwdThreadArg();
3230  one_arg->res = res;
3231  one_arg->res_grad = res_grad;
3232  one_arg->off = 0;
3233  one_arg->coils_els = coils;
3234  one_arg->client = client;
3235  one_arg->s = NULL;
3236  one_arg->fixed_ori = fixed_ori;
3237  one_arg->field_pot = field;
3238  one_arg->vec_field_pot = vec_field;
3239  one_arg->field_pot_grad = field_grad;
3240 
3241  if (nproc < 2)
3242  use_threads = false;
3243 
3244  if (use_threads) {
3245  int nthread = (fixed_ori || vec_field || nproc < 6) ? nspace : 3*nspace;
3246  QList <FwdThreadArg*> args; //fwdThreadArg *args = MALLOC_40(nthread,fwdThreadArg);
3247  int stat;
3248  /*
3249  * We need copies to allocate separate workspace for each thread
3250  */
3251  if (fixed_ori || vec_field || nproc < 6) {
3252  for (k = 0, off = 0; k < nthread; k++) {
3253  FwdThreadArg* t_arg = FwdThreadArg::create_meg_multi_thread_duplicate(one_arg,bem_model != NULL);
3254  t_arg->s = spaces[k];
3255  t_arg->off = off;
3256  off = fixed_ori ? off + spaces[k]->nuse : off + 3*spaces[k]->nuse;
3257  args.append(t_arg);
3258  }
3259  printf("%d processors. I will use one thread for each of the %d source spaces.\n",
3260  nproc,nspace);
3261  }
3262  else {
3263  for (k = 0, off = 0, q = 0; k < nspace; k++) {
3264  for (p = 0; p < 3; p++,q++) {
3265  FwdThreadArg* t_arg = FwdThreadArg::create_meg_multi_thread_duplicate(one_arg,bem_model != NULL);
3266  t_arg->s = spaces[k];
3267  t_arg->off = off;
3268  t_arg->comp = p;
3269  args.append(t_arg);
3270  }
3271  off = fixed_ori ? off + spaces[k]->nuse : off + 3*spaces[k]->nuse;
3272  }
3273  printf("%d processors. I will use %d threads : %d source spaces x 3 source components.\n",
3274  nproc,nthread,nspace);
3275  }
3276  printf("Computing MEG at %d source locations (%s orientations)...",
3277  nsource,fixed_ori ? "fixed" : "free");
3278  /*
3279  * Ready to start the threads & Wait for them to complete
3280  */
3281  QtConcurrent::blockingMap(args, meg_eeg_fwd_one_source_space);
3282  /*
3283  * Check the results
3284  */
3285  for (k = 0, stat = OK; k < nthread; k++)
3286  if (args[k]->stat != OK) {
3287  stat = FAIL;
3288  break;
3289  }
3290  for (k = 0; k < args.size(); k++)//nthread == args.size()
3291  FwdThreadArg::free_meg_multi_thread_duplicate(args[k],bem_model != NULL);
3292  if (stat != OK)
3293  goto bad;
3294  }
3295  else {
3296  printf("Computing MEG at %d source locations (%s orientations, no threads)...",
3297  nsource,fixed_ori ? "fixed" : "free");
3298  for (k = 0, off = 0; k < nspace; k++) {
3299  one_arg->s = spaces[k];
3300  one_arg->off = off;
3301  meg_eeg_fwd_one_source_space(one_arg);
3302  if (one_arg->stat != OK)
3303  goto bad;
3304  off = fixed_ori ? off + one_arg->s->nuse : off + 3*one_arg->s->nuse;
3305  }
3306  }
3307  printf("done.\n");
3308  {
3309  QStringList orig_names;
3310  for (k = 0; k < nmeg; k++)
3311  orig_names.append(coils->coils[k]->chname);
3312  names = orig_names;
3313  }
3314  if(one_arg)
3315  delete one_arg;
3316  if(comp)
3317  delete comp;
3318 
3319  // Store solution in fiff named matrix
3320  nrow = fixed_ori ? nsource : 3*nsource;
3321  matRes.conservativeResize(nrow,nmeg);
3322  for(int i = 0; i < nrow; i++) {
3323  for(int j = 0; j < nmeg; j++) {
3324  matRes(i,j) = res[i][j];
3325  }
3326  }
3327  resp.nrow = nrow;
3328  resp.ncol = nmeg;
3329  resp.row_names = emptyList;
3330  resp.col_names = names;
3331  resp.data = matRes;
3332  resp.transpose_named_matrix();
3333 
3334  if (bDoGrad && res_grad) {
3335  nrow = fixed_ori ? 3*nsource : 3*3*nsource;
3336  matResGrad = MatrixXd::Zero(nrow,nmeg);
3337  for(int i = 0; i < nrow; i++) {
3338  for(int j = 0; j < nmeg; j++) {
3339  matResGrad(i,j) = res_grad[i][j];
3340  }
3341  }
3342  resp_grad.nrow = nrow;
3343  resp_grad.ncol = nmeg;
3344  resp_grad.row_names = emptyList;
3345  resp_grad.col_names = names;
3346  resp_grad.data = matResGrad;
3347  resp_grad.transpose_named_matrix();
3348  }
3349  return OK;
3350 
3351 bad : {
3352  if(one_arg)
3353  delete one_arg;
3354  if(comp)
3355  delete comp;
3356  FREE_CMATRIX_40(res);
3357  FREE_CMATRIX_40(res_grad);
3358  return FAIL;
3359  }
3360 }
3361 
3362 //=============================================================================================================
3363 
3365  int nspace,
3366  FwdCoilSet *els,
3367  bool fixed_ori,
3368  FwdBemModel *bem_model,
3369  FwdEegSphereModel *m,
3370  bool use_threads,
3371  FiffNamedMatrix& resp,
3372  FiffNamedMatrix& resp_grad,
3373  bool bDoGrad)
3374 /*
3375  * Compute the EEG forward solution
3376  * Use either the sphere model or BEM in the calculations
3377  */
3378 {
3379  float **res = NULL; /* The forward solution matrix */
3380  float **res_grad = NULL; /* The gradient with respect to the dipole position */
3381  MatrixXd matRes;
3382  MatrixXd matResGrad;
3383  int nrow = 0;
3384  fwdFieldFunc pot; /* Computes the potentials for one dipole orientation */
3385  fwdVecFieldFunc vec_pot; /* Computes the potentials for all dipole orientations */
3386  fwdFieldGradFunc pot_grad; /* Computes the potential and gradient with respect to dipole position
3387  * for one dipole orientation */
3388  int nsource; /* Total number of sources */
3389  int neeg = els->ncoil; /* Number of channels */
3390  int k,p,q,off;
3391  QStringList names; /* Channel names */
3392  void *client;
3393  FwdThreadArg* one_arg = NULL;
3394  int nproc = QThread::idealThreadCount();
3395  QStringList emptyList;
3396  /*
3397  * Count the sources
3398  */
3399  for (k = 0, nsource = 0; k < nspace; k++)
3400  nsource += spaces[k]->nuse;
3401 
3402  if (bem_model) {
3403  if (fwd_bem_specify_els(bem_model,els) == FAIL)
3404  goto bad;
3405  client = bem_model;
3406  pot = fwd_bem_pot_els;
3407  vec_pot = NULL;
3408 #ifdef TEST
3409  printf("Using differences.\n");
3410  pot_grad = my_bem_pot_grad;
3411 #else
3412  pot_grad = fwd_bem_pot_grad_els;
3413 #endif
3414  }
3415  else {
3416  if (m->nfit == 0) {
3417  printf("Using the standard series expansion for a multilayer sphere model for EEG\n");
3418  pot = FwdEegSphereModel::fwd_eeg_multi_spherepot_coil1;
3419  vec_pot = NULL;
3420  pot_grad = NULL;
3421  }
3422  else {
3423  printf("Using the equivalent source approach in the homogeneous sphere for EEG\n");
3426  pot_grad = FwdEegSphereModel::fwd_eeg_spherepot_grad_coil;
3427  }
3428  client = m;
3429  }
3430  /*
3431  * Allocate space for the solution
3432  */
3433  if (fixed_ori)
3434  res = ALLOC_CMATRIX_40(nsource,neeg);
3435  else
3436  res = ALLOC_CMATRIX_40(3*nsource,neeg);
3437  if (bDoGrad) {
3438  if (!pot_grad) {
3439  qCritical("EEG gradient calculation function not available");
3440  goto bad;
3441  }
3442  if (fixed_ori)
3443  res_grad = ALLOC_CMATRIX_40(3*nsource,neeg);
3444  else
3445  res_grad = ALLOC_CMATRIX_40(3*3*nsource,neeg);
3446  }
3447  /*
3448  * Set up the argument for the field computation
3449  */
3450  one_arg = new FwdThreadArg();
3451  one_arg->res = res;
3452  one_arg->res_grad = res_grad;
3453  one_arg->off = 0;
3454  one_arg->coils_els = els;
3455  one_arg->client = client;
3456  one_arg->s = NULL;
3457  one_arg->fixed_ori = fixed_ori;
3458  one_arg->field_pot = pot;
3459  one_arg->vec_field_pot = vec_pot;
3460  one_arg->field_pot_grad = pot_grad;
3461 
3462  if (nproc < 2)
3463  use_threads = false;
3464 
3465  if (use_threads) {
3466  int nthread = (fixed_ori || vec_pot || nproc < 6) ? nspace : 3*nspace;
3467  QList <FwdThreadArg*> args; //FwdThreadArg* *args = MALLOC_40(nthread,FwdThreadArg*);
3468  int stat;
3469  /*
3470  * We need copies to allocate separate workspace for each thread
3471  */
3472  if (fixed_ori || vec_pot || nproc < 6) {
3473  for (k = 0, off = 0; k < nthread; k++) {
3474  FwdThreadArg* t_arg = FwdThreadArg::create_eeg_multi_thread_duplicate(one_arg,bem_model != NULL);
3475  t_arg->s = spaces[k];
3476  t_arg->off = off;
3477  off = fixed_ori ? off + spaces[k]->nuse : off + 3*spaces[k]->nuse;
3478  args.append(t_arg);
3479  }
3480  printf("%d processors. I will use one thread for each of the %d source spaces.\n",nproc,nspace);
3481  }
3482  else {
3483  for (k = 0, off = 0, q = 0; k < nspace; k++) {
3484  for (p = 0; p < 3; p++,q++) {
3485  FwdThreadArg* t_arg = FwdThreadArg::create_eeg_multi_thread_duplicate(one_arg,bem_model != NULL);
3486  t_arg->s = spaces[k];
3487  t_arg->off = off;
3488  t_arg->comp = p;
3489  args.append(t_arg);
3490  }
3491  off = fixed_ori ? off + spaces[k]->nuse : off + 3*spaces[k]->nuse;
3492  }
3493  printf("%d processors. I will use %d threads : %d source spaces x 3 source components.\n",nproc,nthread,nspace);
3494  }
3495  printf("Computing EEG at %d source locations (%s orientations)...",
3496  nsource,fixed_ori ? "fixed" : "free");
3497  /*
3498  * Ready to start the threads & Wait for them to complete
3499  */
3500  QtConcurrent::blockingMap(args, meg_eeg_fwd_one_source_space);
3501  /*
3502  * Check the results
3503  */
3504  for (k = 0, stat = OK; k < nthread; k++)
3505  if (args[k]->stat != OK) {
3506  stat = FAIL;
3507  break;
3508  }
3509  for (k = 0; k < args.size(); k++)//nthread == args.size()
3510  FwdThreadArg::free_eeg_multi_thread_duplicate(args[k],bem_model != NULL);
3511  if (stat != OK)
3512  goto bad;
3513  }
3514  else {
3515  printf("Computing EEG at %d source locations (%s orientations, no threads)...",
3516  nsource,fixed_ori ? "fixed" : "free");
3517  for (k = 0, off = 0; k < nspace; k++) {
3518  one_arg->s = spaces[k];
3519  one_arg->off = off;
3520  meg_eeg_fwd_one_source_space(one_arg);
3521  if (one_arg->stat != OK)
3522  goto bad;
3523  off = fixed_ori ? off + one_arg->s->nuse : off + 3*one_arg->s->nuse;
3524  }
3525  }
3526  printf("done.\n");
3527  {
3528  QStringList orig_names;
3529  for (k = 0; k < neeg; k++)
3530  orig_names.append(els->coils[k]->chname);
3531  names = orig_names;
3532  }
3533  if(one_arg)
3534  delete one_arg;
3535 
3536  // Store solution in fiff named matrix
3537  nrow = fixed_ori ? nsource : 3*nsource;
3538  matRes.conservativeResize(nrow,neeg);
3539  for(int i = 0; i < nrow; i++) {
3540  for(int j = 0; j < neeg; j++) {
3541  matRes(i,j) = res[i][j];
3542  }
3543  }
3544  resp.nrow = nrow;
3545  resp.ncol = neeg;
3546  resp.row_names = emptyList;
3547  resp.col_names = names;
3548  resp.data = matRes;
3549  resp.transpose_named_matrix();
3550 
3551  if (bDoGrad && res_grad) {
3552  matResGrad = MatrixXd::Zero(fixed_ori ? 3*nsource : 3*3*nsource,neeg);
3553  for(int i = 0; i < nrow; i++) {
3554  for(int j = 0; j < neeg; j++) {
3555  matResGrad(i,j) = res_grad[i][j];
3556  }
3557  }
3558  resp_grad.nrow = nrow;
3559  resp_grad.ncol = neeg;
3560  resp_grad.row_names = emptyList;
3561  resp_grad.col_names = names;
3562  resp_grad.data = matResGrad;
3563  resp_grad.transpose_named_matrix();
3564  }
3565  return OK;
3566 
3567 bad : {
3568  if(one_arg)
3569  delete one_arg;
3570  FREE_CMATRIX_40(res);
3571  FREE_CMATRIX_40(res_grad);
3572  return FAIL;
3573  }
3574 }
3575 
3576 //=============================================================================================================
3577 
3578 #define EPS 1e-5 /* Points closer to origin than this many
3579  meters are considered to be at the
3580  origin */
3581 #define CEPS 1e-5
3582 
3583 int FwdBemModel::fwd_sphere_field(float *rd, float Q[], FwdCoilSet *coils, float Bval[], void *client) /* Client data will be the sphere model origin */
3584 {
3585  /* This version uses Jukka Sarvas' field computation
3586  for details, see
3587 
3588  Jukka Sarvas:
3589 
3590  Basic mathematical and electromagnetic concepts
3591  of the biomagnetic inverse problem,
3592 
3593  Phys. Med. Biol. 1987, Vol. 32, 1, 11-22
3594 
3595  The formulas have been manipulated for efficient computation
3596  by Matti Hamalainen, February 1990
3597 
3598  */
3599  float *r0 = (float *)client; /* The sphere model origin */
3600  float v[3],a_vec[3];
3601  float a,a2,r,r2;
3602  float ar,ar0,rr0;
3603  float vr,ve,re,r0e;
3604  float F,g0,gr,result,sum;
3605  int j,k,p;
3606  FwdCoil* this_coil;
3607  float *this_pos,*this_dir; /* These point to the coil structure! */
3608  int np;
3609  float myrd[3];
3610  float pos[3];
3611  /*
3612  * Shift to the sphere model coordinates
3613  */
3614  for (p = 0; p < 3; p++)
3615  myrd[p] = rd[p] - r0[p];
3616  rd = myrd;
3617  /*
3618  * Check for a dipole at the origin
3619  */
3620  for (k = 0 ; k < coils->ncoil ; k++)
3621  if (FWD_IS_MEG_COIL(coils->coils[k]->coil_class))
3622  Bval[k] = 0.0;
3623  r = VEC_LEN_40(rd);
3624  if (r > EPS) { /* The hard job */
3625 
3626  CROSS_PRODUCT_40(Q,rd,v);
3627 
3628  for (k = 0; k < coils->ncoil; k++) {
3629  this_coil = coils->coils[k];
3630  if (FWD_IS_MEG_COIL(this_coil->type)) {
3631 
3632  np = this_coil->np;
3633 
3634  for (j = 0, sum = 0.0; j < np; j++) {
3635 
3636  this_pos = this_coil->rmag[j];
3637  this_dir = this_coil->cosmag[j];
3638 
3639  for (p = 0; p < 3; p++)
3640  pos[p] = this_pos[p] - r0[p];
3641  this_pos = pos;
3642  result = 0.0;
3643 
3644  /* Vector from dipole to the field point */
3645 
3646  VEC_DIFF_40(rd,this_pos,a_vec);
3647 
3648  /* Compute the dot products needed */
3649 
3650  a2 = VEC_DOT_40(a_vec,a_vec); a = sqrt(a2);
3651 
3652  if (a > 0.0) {
3653  r2 = VEC_DOT_40(this_pos,this_pos); r = sqrt(r2);
3654  if (r > 0.0) {
3655  rr0 = VEC_DOT_40(this_pos,rd);
3656  ar = (r2-rr0);
3657  if (std::fabs(ar/(a*r)+1.0) > CEPS) { /* There is a problem on the negative 'z' axis if the dipole location
3658  * and the field point are on the same line */
3659  ar0 = ar/a;
3660 
3661  ve = VEC_DOT_40(v,this_dir); vr = VEC_DOT_40(v,this_pos);
3662  re = VEC_DOT_40(this_pos,this_dir); r0e = VEC_DOT_40(rd,this_dir);
3663 
3664  /* The main ingredients */
3665 
3666  F = a*(r*a + ar);
3667  gr = a2/r + ar0 + 2.0*(a+r);
3668  g0 = a + 2*r + ar0;
3669 
3670  /* Mix them together... */
3671 
3672  sum = sum + this_coil->w[j]*(ve*F + vr*(g0*r0e - gr*re))/(F*F);
3673  }
3674  }
3675  }
3676  } /* All points done */
3677  Bval[k] = MAG_FACTOR*sum;
3678  }
3679  }
3680  }
3681  return OK; /* Happy conclusion: this works always */
3682 }
3683 
3684 //=============================================================================================================
3685 
3686 int FwdBemModel::fwd_sphere_field_vec(float *rd, FwdCoilSet *coils, float **Bval, void *client) /* Client data will be the sphere model origin */
3687 {
3688  /* This version uses Jukka Sarvas' field computation
3689  for details, see
3690 
3691  Jukka Sarvas:
3692 
3693  Basic mathematical and electromagnetic concepts
3694  of the biomagnetic inverse problem,
3695 
3696  Phys. Med. Biol. 1987, Vol. 32, 1, 11-22
3697 
3698  The formulas have been manipulated for efficient computation
3699  by Matti Hamalainen, February 1990
3700 
3701  The idea of matrix kernels is from
3702 
3703  Mosher, Leahy, and Lewis: EEG and MEG: Forward Solutions for Inverse Methods
3704 
3705  which has been simplified here using standard vector notation
3706 
3707  */
3708  float *r0 = (float *)client; /* The sphere model origin */
3709  float a_vec[3],v1[3],v2[3];
3710  float a,a2,r,r2;
3711  float ar,ar0,rr0;
3712  float re,r0e;
3713  float F,g0,gr,g,sum[3];
3714  int j,k,p;
3715  FwdCoil* this_coil;
3716  float *this_pos,*this_dir; /* These point to the coil structure! */
3717  int np;
3718  float myrd[3];
3719  float pos[3];
3720  /*
3721  * Shift to the sphere model coordinates
3722  */
3723  for (p = 0; p < 3; p++)
3724  myrd[p] = rd[p] - r0[p];
3725  rd = myrd;
3726  /*
3727  * Check for a dipole at the origin
3728  */
3729  r = VEC_LEN_40(rd);
3730  for (k = 0; k < coils->ncoil; k++) {
3731  this_coil = coils->coils[k];
3732  if (FWD_IS_MEG_COIL(this_coil->coil_class)) {
3733  if (r < EPS) {
3734  Bval[0][k] = Bval[1][k] = Bval[2][k] = 0.0;
3735  }
3736  else { /* The hard job */
3737 
3738  np = this_coil->np;
3739  sum[0] = sum[1] = sum[2] = 0.0;
3740 
3741  for (j = 0; j < np; j++) {
3742 
3743  this_pos = this_coil->rmag[j];
3744  this_dir = this_coil->cosmag[j];
3745 
3746  for (p = 0; p < 3; p++)
3747  pos[p] = this_pos[p] - r0[p];
3748  this_pos = pos;
3749 
3750  /* Vector from dipole to the field point */
3751 
3752  VEC_DIFF_40(rd,this_pos,a_vec);
3753 
3754  /* Compute the dot products needed */
3755 
3756  a2 = VEC_DOT_40(a_vec,a_vec); a = sqrt(a2);
3757 
3758  if (a > 0.0) {
3759  r2 = VEC_DOT_40(this_pos,this_pos); r = sqrt(r2);
3760  if (r > 0.0) {
3761  rr0 = VEC_DOT_40(this_pos,rd);
3762  ar = (r2-rr0);
3763  if (std::fabs(ar/(a*r)+1.0) > CEPS) { /* There is a problem on the negative 'z' axis if the dipole location
3764  * and the field point are on the same line */
3765 
3766  /* The main ingredients */
3767 
3768  ar0 = ar/a;
3769  F = a*(r*a + ar);
3770  gr = a2/r + ar0 + 2.0*(a+r);
3771  g0 = a + 2*r + ar0;
3772 
3773  re = VEC_DOT_40(this_pos,this_dir); r0e = VEC_DOT_40(rd,this_dir);
3774  CROSS_PRODUCT_40(rd,this_dir,v1);
3775  CROSS_PRODUCT_40(rd,this_pos,v2);
3776 
3777  g = (g0*r0e - gr*re)/(F*F);
3778  /*
3779  * Mix them together...
3780  */
3781  for (p = 0; p < 3; p++)
3782  sum[p] = sum[p] + this_coil->w[j]*(v1[p]/F + v2[p]*g);
3783  }
3784  }
3785  }
3786  } /* All points done */
3787  for (p = 0; p < 3; p++)
3788  Bval[p][k] = MAG_FACTOR*sum[p];
3789  }
3790  }
3791  }
3792  return OK; /* Happy conclusion: this works always */
3793 }
3794 
3795 //=============================================================================================================
3796 
3797 int FwdBemModel::fwd_sphere_field_grad(float *rd, float Q[], FwdCoilSet *coils, float Bval[], float xgrad[], float ygrad[], float zgrad[], void *client) /* Client data to be passed to some foward modelling routines */
3798 /*
3799  * Compute the derivatives of the sphere model field with respect to
3800  * dipole coordinates
3801  */
3802 {
3803  /* This version uses Jukka Sarvas' field computation
3804  for details, see
3805 
3806  Jukka Sarvas:
3807 
3808  Basic mathematical and electromagnetic concepts
3809  of the biomagnetic inverse problem,
3810 
3811  Phys. Med. Biol. 1987, Vol. 32, 1, 11-22
3812 
3813  The formulas have been manipulated for efficient computation
3814  by Matti Hamalainen, February 1990
3815 
3816  */
3817 
3818  float v[3],a_vec[3];
3819  float a,a2,r,r2;
3820  float ar,rr0;
3821  float vr,ve,re,r0e;
3822  float F,g0,gr,result,G,F2;
3823 
3824  int j,k,p;
3825  float huu;
3826  float ggr[3],gg0[3]; /* Gradient of gr & g0 */
3827  float ga[3]; /* Grapdient of a */
3828  float gar[3]; /* Gradient of ar */
3829  float gFF[3]; /* Gradient of F divided by F */
3830  float gresult[3];
3831  float eQ[3],rQ[3]; /* e x Q and r x Q */
3832  int do_field = 0;
3833  FwdCoil* this_coil;
3834  float *this_pos,*this_dir;
3835  int np;
3836  float myrd[3];
3837  float pos[3];
3838  float *r0 = (float *)client; /* The sphere model origin */
3839  /*
3840  * Shift to the sphere model coordinates
3841  */
3842  for (p = 0; p < 3; p++)
3843  myrd[p] = rd[p] - r0[p];
3844  rd = myrd;
3845 
3846  if (Bval)
3847  do_field = 1;
3848 
3849  /* Check for a dipole at the origin */
3850 
3851  r = VEC_LEN_40(rd);
3852  for (k = 0; k < coils->ncoil ; k++) {
3853  if (FWD_IS_MEG_COIL(coils->coils[k]->coil_class)) {
3854  if (do_field)
3855  Bval[k] = 0.0;
3856  xgrad[k] = 0.0;
3857  ygrad[k] = 0.0;
3858  zgrad[k] = 0.0;
3859  }
3860  }
3861  if (r > EPS) { /* The hard job */
3862 
3863  v[X_40] = Q[Y_40]*rd[Z_40] - Q[Z_40]*rd[Y_40];
3864  v[Y_40] = -Q[X_40]*rd[Z_40] + Q[Z_40]*rd[X_40];
3865  v[Z_40] = Q[X_40]*rd[Y_40] - Q[Y_40]*rd[X_40];
3866 
3867  for (k = 0 ; k < coils->ncoil ; k++) {
3868 
3869  this_coil = coils->coils[k];
3870 
3871  if (FWD_IS_MEG_COIL(this_coil->type)) {
3872 
3873  np = this_coil->np;
3874 
3875  for (j = 0; j < np; j++) {
3876 
3877  this_pos = this_coil->rmag[j];
3878  /*
3879  * Shift to the sphere model coordinates
3880  */
3881  for (p = 0; p < 3; p++)
3882  pos[p] = this_pos[p] - r0[p];
3883  this_pos = pos;
3884 
3885  this_dir = this_coil->cosmag[j];
3886 
3887  /* Vector from dipole to the field point */
3888 
3889  a_vec[X_40] = this_pos[X_40] - rd[X_40];
3890  a_vec[Y_40] = this_pos[Y_40] - rd[Y_40];
3891  a_vec[Z_40] = this_pos[Z_40] - rd[Z_40];
3892 
3893  /* Compute the dot and cross products needed */
3894 
3895  a2 = VEC_DOT_40(a_vec,a_vec); a = sqrt(a2);
3896  r2 = VEC_DOT_40(this_pos,this_pos); r = sqrt(r2);
3897  rr0 = VEC_DOT_40(this_pos,rd);
3898  ar = (r2 - rr0)/a;
3899 
3900  ve = VEC_DOT_40(v,this_dir); vr = VEC_DOT_40(v,this_pos);
3901  re = VEC_DOT_40(this_pos,this_dir); r0e = VEC_DOT_40(rd,this_dir);
3902 
3903  /* eQ = this_dir x Q */
3904 
3905  eQ[X_40] = this_dir[Y_40]*Q[Z_40] - this_dir[Z_40]*Q[Y_40];
3906  eQ[Y_40] = -this_dir[X_40]*Q[Z_40] + this_dir[Z_40]*Q[X_40];
3907  eQ[Z_40] = this_dir[X_40]*Q[Y_40] - this_dir[Y_40]*Q[X_40];
3908 
3909  /* rQ = this_pos x Q */
3910 
3911  rQ[X_40] = this_pos[Y_40]*Q[Z_40] - this_pos[Z_40]*Q[Y_40];
3912  rQ[Y_40] = -this_pos[X_40]*Q[Z_40] + this_pos[Z_40]*Q[X_40];
3913  rQ[Z_40] = this_pos[X_40]*Q[Y_40] - this_pos[Y_40]*Q[X_40];
3914 
3915  /* The main ingredients */
3916 
3917  F = a*(r*a + r2 - rr0);
3918  F2 = F*F;
3919  gr = a2/r + ar + 2.0*(a+r);
3920  g0 = a + 2.0*r + ar;
3921  G = g0*r0e - gr*re;
3922 
3923  /* Mix them together... */
3924 
3925  result = (ve*F + vr*G)/F2;
3926 
3927  /* The computation of the gradient... */
3928 
3929  huu = 2.0 + 2.0*a/r;
3930  for (p = X_40; p <= Z_40; p++) {
3931  ga[p] = -a_vec[p]/a;
3932  gar[p] = -(ga[p]*ar + this_pos[p])/a;
3933  gg0[p] = ga[p] + gar[p];
3934  ggr[p] = huu*ga[p] + gar[p];
3935  gFF[p] = ga[p]/a - (r*a_vec[p] + a*this_pos[p])/F;
3936  gresult[p] = -2.0*result*gFF[p] + (eQ[p]+gFF[p]*ve)/F +
3937  (rQ[p]*G + vr*(gg0[p]*r0e + g0*this_dir[p] - ggr[p]*re))/F2;
3938  }
3939 
3940  if (do_field)
3941  Bval[k] = Bval[k] + this_coil->w[j]*result;
3942  xgrad[k] = xgrad[k] + this_coil->w[j]*gresult[X_40];
3943  ygrad[k] = ygrad[k] + this_coil->w[j]*gresult[Y_40];
3944  zgrad[k] = zgrad[k] + this_coil->w[j]*gresult[Z_40];
3945  }
3946  if (do_field)
3947  Bval[k] = MAG_FACTOR*Bval[k];
3948  xgrad[k] = MAG_FACTOR*xgrad[k];
3949  ygrad[k] = MAG_FACTOR*ygrad[k];
3950  zgrad[k] = MAG_FACTOR*zgrad[k];
3951  }
3952  }
3953  }
3954  return OK; /* Happy conclusion: this works always */
3955 }
3956 
3957 //=============================================================================================================
3958 
3959 int FwdBemModel::fwd_mag_dipole_field(float *rm, float M[], FwdCoilSet *coils, float Bval[], void *client) /* Client data will be the sphere model origin */
3960 /*
3961  * This is for a specific dipole component
3962  */
3963 {
3964  int j,k,np;
3965  FwdCoil* this_coil;
3966  float sum,diff[3],dist,dist2,dist5,*dir;
3967 
3968  for (k = 0; k < coils->ncoil; k++) {
3969  this_coil = coils->coils[k];
3970  if (FWD_IS_MEG_COIL(this_coil->type)) {
3971  np = this_coil->np;
3972  /*
3973  * Go through all points
3974  */
3975  for (j = 0, sum = 0.0; j < np; j++) {
3976  dir = this_coil->cosmag[j];
3977  VEC_DIFF_40(rm,this_coil->rmag[j],diff);
3978  dist = VEC_LEN_40(diff);
3979  if (dist > EPS) {
3980  dist2 = dist*dist;
3981  dist5 = dist2*dist2*dist;
3982  sum = sum + this_coil->w[j]*(3*VEC_DOT_40(M,diff)*VEC_DOT_40(diff,dir) - dist2*VEC_DOT_40(M,dir))/dist5;
3983  }
3984  } /* All points done */
3985  Bval[k] = MAG_FACTOR*sum;
3986  }
3987  else if (this_coil->type == FWD_COILC_EEG)
3988  Bval[k] = 0.0;
3989  }
3990  return OK;
3991 }
3992 
3993 //=============================================================================================================
3994 
3995 int FwdBemModel::fwd_mag_dipole_field_vec(float *rm, FwdCoilSet *coils, float **Bval, void *client) /* Client data will be the sphere model origin */
3996 /*
3997  * This is for all dipole components
3998  * For EEG this produces a zero result
3999  */
4000 {
4001  int j,k,p,np;
4002  FwdCoil* this_coil;
4003  float sum[3],diff[3],dist,dist2,dist5,*dir;
4004 
4005  for (k = 0; k < coils->ncoil; k++) {
4006  this_coil = coils->coils[k];
4007  if (FWD_IS_MEG_COIL(this_coil->type)) {
4008  np = this_coil->np;
4009  sum[0] = sum[1] = sum[2] = 0.0;
4010  /*
4011  * Go through all points
4012  */
4013  for (j = 0; j < np; j++) {
4014  dir = this_coil->cosmag[j];
4015  VEC_DIFF_40(rm,this_coil->rmag[j],diff);
4016  dist = VEC_LEN_40(diff);
4017  if (dist > EPS) {
4018  dist2 = dist*dist;
4019  dist5 = dist2*dist2*dist;
4020  for (p = 0; p < 3; p++)
4021  sum[p] = sum[p] + this_coil->w[j]*(3*diff[p]*VEC_DOT_40(diff,dir) - dist2*dir[p])/dist5;
4022  }
4023  } /* All points done */
4024  for (p = 0; p < 3; p++)
4025  Bval[p][k] = MAG_FACTOR*sum[p];
4026  }
4027  else if (this_coil->type == FWD_COILC_EEG) {
4028  for (p = 0; p < 3; p++)
4029  Bval[p][k] = 0.0;
4030  }
4031  }
4032  return OK;
4033 }
#define FIFFV_MNE_COORD_MRI_VOXEL
#define FIFFV_BEM_APPROX_CONST
Definition: fiff_file.h:781
FiffStream class declaration.
Fwd Thread Argument (FwdThreadArg) class declaration.
float ** cosmag
Definition: fwd_coil.h:181
Filter Thread Argument Description.
#define FIFFV_MNE_COORD_CTF_DEVICE
#define FIFFV_BEM_APPROX_LINEAR
Definition: fiff_file.h:782
float ** rmag
Definition: fwd_coil.h:180
QSharedPointer< FiffDirNode > SPtr
Definition: fiff_dir_node.h:77
#define FIFF_BEM_APPROX
Definition: fiff_file.h:745
FwdCoilSet description.
Definition: fwd_coil_set.h:75
Forward BEM Solution (FwdBemSolution) class declaration.
#define FIFFV_MNE_COORD_FS_TAL_LTZ
Holds the BEM model definition.
QSharedPointer< FiffTag > SPtr
Definition: fiff_tag.h:152
Triangle data.
Definition: mne_triangle.h:75
static QString fwd_bem_make_bem_sol_name(const QString &name)
FwdBemModel class declaration.
This defines a surface.
FwdEegSphereModel class declaration.
This defines a source space.
FwdCompData class declaration.
One MNE CTF Compensation Data Set description.
#define FIFFB_BEM
Definition: fiff_file.h:403
FIFF File I/O routines.
Definition: fiff_stream.h:104
Coordinate transformation descriptor.
#define FIFF_BEM_POT_SOLUTION
Definition: fiff_file.h:744
FiffNamedMatrix class declaration.
static int fwd_eeg_spherepot_coil_vec(float *rd, FwdCoilSet *els, float **Vval_vec, void *client)
FwdCoil description.
Definition: fwd_coil.h:88
MneTriangle class declaration.
QString chname
Definition: fwd_coil.h:167
This structure is used in the compensated field calculations.
Definition: fwd_comp_data.h:87
#define FIFFV_MNE_COORD_RAS
FwdCoilSet * dup_coil_set(const FIFFLIB::FiffCoordTransOld *t) const
Mapping from infinite medium potentials to a particular set of coils or electrodes.
QSharedPointer< FiffStream > SPtr
Definition: fiff_stream.h:107
#define FIFFV_MNE_COORD_CTF_HEAD
#define FIFFV_MNE_COORD_MNI_TAL
static int compute_forward_eeg(MNELIB::MneSourceSpaceOld **spaces, int nspace, FwdCoilSet *els, bool fixed_ori, FwdBemModel *bem_model, FwdEegSphereModel *m, bool use_threads, FIFFLIB::FiffNamedMatrix &resp, FIFFLIB::FiffNamedMatrix &resp_grad, bool bDoGrad)
#define FIFFV_MNE_COORD_FS_TAL_GTZ
Electric Current Dipole description.
static int compute_forward_meg(MNELIB::MneSourceSpaceOld **spaces, int nspace, FwdCoilSet *coils, FwdCoilSet *comp_coils, MNELIB::MneCTFCompDataSet *comp_data, bool fixed_ori, FwdBemModel *bem_model, Eigen::Vector3f *r0, bool use_threads, FIFFLIB::FiffNamedMatrix &resp, FIFFLIB::FiffNamedMatrix &resp_grad, bool bDoGRad)
MneSurfaceOld class declaration.
MneSourceSpaceOld class declaration.
static int fwd_eeg_spherepot_coil(float *rd, float *Q, FwdCoilSet *els, float *Vval, void *client)