MNE-CPP  0.1.9
A Framework for Electrophysiology
dipole_fit_data.cpp
1 
2 #include <fwd/fwd_types.h>
3 
4 #include "dipole_fit_data.h"
5 #include "guess_data.h"
6 #include "../c/mne_meas_data.h"
7 #include "../c/mne_meas_data_set.h"
8 #include <mne/c/mne_proj_item.h>
9 #include <mne/c/mne_cov_matrix.h>
10 #include "ecd.h"
11 
12 #include <fiff/fiff_stream.h>
13 #include <fwd/fwd_bem_model.h>
14 #include <mne/c/mne_surface_old.h>
15 
16 #include <fwd/fwd_comp_data.h>
17 
18 #include <Eigen/Dense>
19 
20 #include <QFile>
21 #include <QCoreApplication>
22 #include <QDebug>
23 
24 #define _USE_MATH_DEFINES
25 #include <math.h>
26 
27 using namespace Eigen;
28 using namespace FIFFLIB;
29 using namespace MNELIB;
30 using namespace FWDLIB;
31 using namespace INVERSELIB;
32 
33 //============================= ctf_types.h =============================
34 
35 #ifndef FIFFV_COIL_CTF_GRAD
36 #define FIFFV_COIL_CTF_GRAD 5001
37 #endif
38 
39 #ifndef FIFFV_COIL_CTF_REF_MAG
40 #define FIFFV_COIL_CTF_REF_MAG 5002
41 #endif
42 
43 #ifndef FIFFV_COIL_CTF_REF_GRAD
44 #define FIFFV_COIL_CTF_REF_GRAD 5003
45 #endif
46 
47 #ifndef FIFFV_COIL_CTF_OFFDIAG_REF_GRAD
48 #define FIFFV_COIL_CTF_OFFDIAG_REF_GRAD 5004
49 #endif
50 
51 #ifndef TRUE
52 #define TRUE 1
53 #endif
54 
55 #ifndef FALSE
56 #define FALSE 0
57 #endif
58 
59 #ifndef FAIL
60 #define FAIL -1
61 #endif
62 
63 #ifndef OK
64 #define OK 0
65 #endif
66 
67 #define X_3 0
68 #define Y_3 1
69 #define Z_3 2
70 
71 #define MALLOC_3(x,t) (t *)malloc((x)*sizeof(t))
72 #define REALLOC_3(x,y,t) (t *)((x == NULL) ? malloc((y)*sizeof(t)) : realloc((x),(y)*sizeof(t)))
73 #define FREE_3(x) if ((char *)(x) != NULL) free((char *)(x))
74 
75 /*
76  * Float, double, and int arrays
77  */
78 #define ALLOC_FLOAT_3(x) MALLOC_3(x,float)
79 
80 #define ALLOC_DCMATRIX_3(x,y) mne_dmatrix_3((x),(y))
81 #define ALLOC_CMATRIX_3(x,y) mne_cmatrix_3((x),(y))
82 #define FREE_CMATRIX_3(m) mne_free_cmatrix_3((m))
83 #define FREE_DCMATRIX_3(m) mne_free_dcmatrix_3((m))
84 
85 static void matrix_error_3(int kind, int nr, int nc)
86 
87 {
88  if (kind == 1)
89  printf("Failed to allocate memory pointers for a %d x %d matrix\n",nr,nc);
90  else if (kind == 2)
91  printf("Failed to allocate memory for a %d x %d matrix\n",nr,nc);
92  else
93  printf("Allocation error for a %d x %d matrix\n",nr,nc);
94  if (sizeof(void *) == 4) {
95  printf("This is probably because you seem to be using a computer with 32-bit architecture.\n");
96  printf("Please consider moving to a 64-bit platform.");
97  }
98  printf("Cannot continue. Sorry.\n");
99  exit(1);
100 }
101 
102 float **mne_cmatrix_3 (int nr,int nc)
103 
104 {
105  int i;
106  float **m;
107  float *whole;
108 
109  m = MALLOC_3(nr,float *);
110  if (!m) matrix_error_3(1,nr,nc);
111  whole = MALLOC_3(nr*nc,float);
112  if (!whole) matrix_error_3(2,nr,nc);
113 
114  for(i=0;i<nr;i++)
115  m[i] = whole + i*nc;
116  return m;
117 }
118 
119 double **mne_dmatrix_3(int nr, int nc)
120 
121 {
122  int i;
123  double **m;
124  double *whole;
125 
126  m = MALLOC_3(nr,double *);
127  if (!m) matrix_error_3(1,nr,nc);
128  whole = MALLOC_3(nr*nc,double);
129  if (!whole) matrix_error_3(2,nr,nc);
130 
131  for(i=0;i<nr;i++)
132  m[i] = whole + i*nc;
133  return m;
134 }
135 
136 void mne_free_dcmatrix_3 (double **m)
137 
138 {
139  if (m) {
140  FREE_3(*m);
141  FREE_3(m);
142  }
143 }
144 
145 /*
146  * Dot product and length
147  */
148 #define VEC_DOT_3(x,y) ((x)[X_3]*(y)[X_3] + (x)[Y_3]*(y)[Y_3] + (x)[Z_3]*(y)[Z_3])
149 #define VEC_LEN_3(x) sqrt(VEC_DOT_3(x,x))
150 
151 /*
152  * Others...
153  */
154 
155 #define VEC_DIFF_3(from,to,diff) {\
156  (diff)[X_3] = (to)[X_3] - (from)[X_3];\
157  (diff)[Y_3] = (to)[Y_3] - (from)[Y_3];\
158  (diff)[Z_3] = (to)[Z_3] - (from)[Z_3];\
159  }
160 
161 #define VEC_COPY_3(to,from) {\
162  (to)[X_3] = (from)[X_3];\
163  (to)[Y_3] = (from)[Y_3];\
164  (to)[Z_3] = (from)[Z_3];\
165  }
166 
167 #define CROSS_PRODUCT_3(x,y,xy) {\
168  (xy)[X_3] = (x)[Y_3]*(y)[Z_3]-(y)[Y_3]*(x)[Z_3];\
169  (xy)[Y_3] = -((x)[X_3]*(y)[Z_3]-(y)[X_3]*(x)[Z_3]);\
170  (xy)[Z_3] = (x)[X_3]*(y)[Y_3]-(y)[X_3]*(x)[Y_3];\
171  }
172 
173 //============================= mne_matop.c =============================
174 
175 #define MIN_3(a,b) ((a) < (b) ? (a) : (b))
176 
177 void mne_transpose_square_3(float **mat, int n)
178 /*
179  * In-place transpose of a square matrix
180  */
181 {
182  int j,k;
183  float val;
184 
185  for (j = 1; j < n; j++)
186  for (k = 0; k < j; k++) {
187  val = mat[j][k];
188  mat[j][k] = mat[k][j];
189  mat[k][j] = val;
190  }
191  return;
192 }
193 
194 float mne_dot_vectors_3(float *v1,
195  float *v2,
196  int nn)
197 
198 {
199 #ifdef BLAS
200  int one = 1;
201  float res = sdot(&nn,v1,&one,v2,&one);
202  return res;
203 #else
204  float res = 0.0;
205  int k;
206 
207  for (k = 0; k < nn; k++)
208  res = res + v1[k]*v2[k];
209  return res;
210 #endif
211 }
212 
213 void mne_add_scaled_vector_to_3(float *v1,float scale, float *v2,int nn)
214 
215 {
216 #ifdef BLAS
217  float fscale = scale;
218  int one = 1;
219  saxpy(&nn,&fscale,v1,&one,v2,&one);
220 #else
221  int k;
222  for (k = 0; k < nn; k++)
223  v2[k] = v2[k] + scale*v1[k];
224 #endif
225  return;
226 }
227 
228 double **mne_dmatt_dmat_mult2_3 (double **m1,double **m2, int d1,int d2,int d3)
229 /* Matrix multiplication
230  * result(d1 x d3) = m1(d2 x d1)^T * m2(d2 x d3) */
231 
232 {
233  double **result = ALLOC_DCMATRIX_3(d1,d3);
234 #ifdef BLAS
235  char *transa = "N";
236  char *transb = "T";
237  double zero = 0.0;
238  double one = 1.0;
239 
240  dgemm (transa,transb,&d3,&d1,&d2,
241  &one,m2[0],&d3,m1[0],&d1,&zero,result[0],&d3);
242 
243  return result;
244 #else
245  int j,k,p;
246  double sum;
247 
248  for (j = 0; j < d1; j++)
249  for (k = 0; k < d3; k++) {
250  sum = 0.0;
251  for (p = 0; p < d2; p++)
252  sum = sum + m1[p][j]*m2[p][k];
253  result[j][k] = sum;
254  }
255  return result;
256 #endif
257 }
258 
259 float **mne_mat_mat_mult_3 (float **m1,float **m2,int d1,int d2,int d3)
260 /* Matrix multiplication
261  * result(d1 x d3) = m1(d1 x d2) * m2(d2 x d3) */
262 
263 {
264 #ifdef BLAS
265  float **result = ALLOC_CMATRIX_3(d1,d3);
266  char *transa = "N";
267  char *transb = "N";
268  float zero = 0.0;
269  float one = 1.0;
270  sgemm (transa,transb,&d3,&d1,&d2,
271  &one,m2[0],&d3,m1[0],&d2,&zero,result[0],&d3);
272  return (result);
273 #else
274  float **result = ALLOC_CMATRIX_3(d1,d3);
275  int j,k,p;
276  float sum;
277 
278  for (j = 0; j < d1; j++)
279  for (k = 0; k < d3; k++) {
280  sum = 0.0;
281  for (p = 0; p < d2; p++)
282  sum = sum + m1[j][p]*m2[p][k];
283  result[j][k] = sum;
284  }
285  return (result);
286 #endif
287 }
288 
289 void mne_mat_vec_mult2_3(float **m,float *v,float *result, int d1,int d2)
290 /*
291  * Matrix multiplication
292  * result(d1) = m(d1 x d2) * v(d2)
293  */
294 
295 {
296  int j;
297 
298  for (j = 0; j < d1; j++)
299  result[j] = mne_dot_vectors_3(m[j],v,d2);
300  return;
301 }
302 
303 //============================= mne_named_matrix.c =============================
304 
305 QString mne_name_list_to_string_3(const QStringList& list)
306 /*
307  * Convert a string array to a colon-separated string
308  */
309 {
310  int nlist = list.size();
311  QString res;
312  if (nlist == 0 || list.isEmpty())
313  return res;
314 // res[0] = '\0';
315  for (int k = 0; k < nlist-1; k++) {
316  res += list[k];
317  res += ":";
318  }
319  res += list[nlist-1];
320  return res;
321 }
322 
323 QString mne_channel_names_to_string_3(const QList<FiffChInfo>& chs, int nch)
324 /*
325  * Make a colon-separated string out of channel names
326  */
327 {
328  QStringList names;
329  QString res;
330  int k;
331 
332  if (nch <= 0)
333  return NULL;
334  for (k = 0; k < nch; k++)
335  names.append(chs[k].ch_name);
336  res = mne_name_list_to_string_3(names);
337  names.clear();
338  return res;
339 }
340 
341 void mne_string_to_name_list_3(const QString& s, QStringList& listp,int &nlistp)
342 /*
343  * Convert a colon-separated list into a string array
344  */
345 {
346  QStringList list;
347 
348  if (!s.isEmpty() && s.size() > 0) {
350  //list = s.split(":");
351  }
352  listp = list;
353  nlistp = list.size();
354  return;
355 }
356 
357 //============================= mne_sparse_matop.c =============================
358 
359 FiffSparseMatrix* mne_convert_to_sparse_3(float **dense, /* The dense matrix to be converted */
360  int nrow, /* Number of rows in the dense matrix */
361  int ncol, /* Number of columns in the dense matrix */
362  int stor_type, /* Either FIFFTS_MC_CCS or FIFFTS_MC_RCS */
363  float small) /* How small elements should be ignored? */
364 /*
365  * Create the compressed row or column storage sparse matrix representation
366  * including a vector containing the nonzero matrix element values,
367  * the row or column pointer vector and the appropriate index vector(s).
368  */
369 {
370  int j,k;
371  int nz;
372  int ptr;
373  FiffSparseMatrix* sparse = NULL;
374  int size;
375 
376  if (small < 0) { /* Automatic scaling */
377  float maxval = 0.0;
378  float val;
379 
380  for (j = 0; j < nrow; j++)
381  for (k = 0; k < ncol; k++) {
382  val = std::fabs(dense[j][k]);
383  if (val > maxval)
384  maxval = val;
385  }
386  if (maxval > 0)
387  small = maxval * std::fabs(small);
388  else
389  small = std::fabs(small);
390  }
391  for (j = 0, nz = 0; j < nrow; j++)
392  for (k = 0; k < ncol; k++) {
393  if (std::fabs(dense[j][k]) > small)
394  nz++;
395  }
396 
397  if (nz <= 0) {
398  printf("No nonzero elements found.");
399  return NULL;
400  }
401  if (stor_type == FIFFTS_MC_CCS) {
402  size = nz*(sizeof(fiff_float_t) + sizeof(fiff_int_t)) +
403  (ncol+1)*(sizeof(fiff_int_t));
404  }
405  else if (stor_type == FIFFTS_MC_RCS) {
406  size = nz*(sizeof(fiff_float_t) + sizeof(fiff_int_t)) +
407  (nrow+1)*(sizeof(fiff_int_t));
408  }
409  else {
410  printf("Unknown sparse matrix storage type: %d",stor_type);
411  return NULL;
412  }
413  sparse = new FiffSparseMatrix;
414  sparse->coding = stor_type;
415  sparse->m = nrow;
416  sparse->n = ncol;
417  sparse->nz = nz;
418  sparse->data = (float *)malloc(size);
419  sparse->inds = (int *)(sparse->data+nz);
420  sparse->ptrs = sparse->inds+nz;
421 
422  if (stor_type == FIFFTS_MC_RCS) {
423  for (j = 0, nz = 0; j < nrow; j++) {
424  ptr = -1;
425  for (k = 0; k < ncol; k++)
426  if (std::fabs(dense[j][k]) > small) {
427  sparse->data[nz] = dense[j][k];
428  if (ptr < 0)
429  ptr = nz;
430  sparse->inds[nz++] = k;
431  }
432  sparse->ptrs[j] = ptr;
433  }
434  sparse->ptrs[nrow] = nz;
435  for (j = nrow - 1; j >= 0; j--) /* Take care of the empty rows */
436  if (sparse->ptrs[j] < 0)
437  sparse->ptrs[j] = sparse->ptrs[j+1];
438  }
439  else if (stor_type == FIFFTS_MC_CCS) {
440  for (k = 0, nz = 0; k < ncol; k++) {
441  ptr = -1;
442  for (j = 0; j < nrow; j++)
443  if (std::fabs(dense[j][k]) > small) {
444  sparse->data[nz] = dense[j][k];
445  if (ptr < 0)
446  ptr = nz;
447  sparse->inds[nz++] = j;
448  }
449  sparse->ptrs[k] = ptr;
450  }
451  sparse->ptrs[ncol] = nz;
452  for (k = ncol-1; k >= 0; k--) /* Take care of the empty columns */
453  if (sparse->ptrs[k] < 0)
454  sparse->ptrs[k] = sparse->ptrs[k+1];
455  }
456  return sparse;
457 }
458 
459 namespace MNELIB
460 {
461 
462 typedef struct {
463  float limit;
464  int report_dim;
465  float *B;
466  double B2;
467  DipoleForward* fwd;
469 
470 }
471 
472 int mne_is_diag_cov_3(MneCovMatrix* c)
473 
474 {
475  return c->cov_diag != NULL;
476 }
477 
478 void mne_scale_vector_3 (double scale,float *v,int nn)
479 
480 {
481 #ifdef BLAS
482  float fscale = scale;
483  int one = 1;
484  sscal(&nn,&fscale,v,&one);
485 #else
486  int k;
487  for (k = 0; k < nn; k++)
488  v[k] = v[k]*scale;
489 #endif
490 }
491 
492 //float
493 Eigen::MatrixXf toFloatEigenMatrix_3(float **mat, const int m, const int n)
494 {
495  Eigen::MatrixXf eigen_mat(m,n);
496 
497  for ( int i = 0; i < m; ++i)
498  for ( int j = 0; j < n; ++j)
499  eigen_mat(i,j) = mat[i][j];
500 
501  return eigen_mat;
502 }
503 
504 void fromFloatEigenMatrix_3(const Eigen::MatrixXf& from_mat, float **& to_mat, const int m, const int n)
505 {
506  for ( int i = 0; i < m; ++i)
507  for ( int j = 0; j < n; ++j)
508  to_mat[i][j] = from_mat(i,j);
509 }
510 
511 void fromFloatEigenMatrix_3(const Eigen::MatrixXf& from_mat, float **& to_mat)
512 {
513  fromFloatEigenMatrix_3(from_mat, to_mat, from_mat.rows(), from_mat.cols());
514 }
515 
516 void fromFloatEigenVector_3(const Eigen::VectorXf& from_vec, float *to_vec, const int n)
517 {
518  for ( int i = 0; i < n; ++i)
519  to_vec[i] = from_vec[i];
520 }
521 
522 float **mne_lu_invert_3(float **mat,int dim)
523 /*
524  * Invert a matrix using the LU decomposition from
525  * LAPACK
526  */
527 {
528  Eigen::MatrixXf eigen_mat = toFloatEigenMatrix_3(mat, dim, dim);
529  Eigen::MatrixXf eigen_mat_inv = eigen_mat.inverse();
530  fromFloatEigenMatrix_3(eigen_mat_inv, mat);
531  return mat;
532 }
533 
534 void mne_free_cmatrix_3 (float **m)
535 {
536  if (m) {
537  FREE_3(*m);
538  FREE_3(m);
539  }
540 }
541 
542 int mne_svd_3(float **mat, /* The matrix */
543  int m,int n, /* m rows n columns */
544  float *sing, /* Singular values (must have size
545  * MIN(m,n)+1 */
546  float **uu, /* Left eigenvectors */
547  float **vv) /* Right eigenvectors */
548 /*
549  * Compute the SVD of mat.
550  * The singular vector calculations depend on whether
551  * or not u and v are given.
552  *
553  * The allocations should be done as follows
554  *
555  * mat = ALLOC_CMATRIX_3(m,n);
556  * vv = ALLOC_CMATRIX_3(MIN(m,n),n);
557  * uu = ALLOC_CMATRIX_3(MIN(m,n),m);
558  * sing = MALLOC_3(MIN(m,n),float);
559  *
560  * mat is modified by this operation
561  *
562  * This simply allocates the workspace and calls the
563  * LAPACK Fortran routine
564  */
565 
566 {
567  int udim = MIN_3(m,n);
568 
569  Eigen::MatrixXf eigen_mat = toFloatEigenMatrix_3(mat, m, n);
570 
571  //ToDo Optimize computation depending of whether uu or vv are defined
572  Eigen::JacobiSVD< Eigen::MatrixXf > svd(eigen_mat ,Eigen::ComputeFullU | Eigen::ComputeFullV);
573 
574  fromFloatEigenVector_3(svd.singularValues(), sing, svd.singularValues().size());
575 
576  if (uu != NULL)
577  fromFloatEigenMatrix_3(svd.matrixU().transpose(), uu, udim, m);
578 
579  if (vv != NULL)
580  fromFloatEigenMatrix_3(svd.matrixV().transpose(), vv, m, n);
581 
582  return 0;
583  // return info;
584 }
585 
586 void mne_free_cov_3(MneCovMatrix* c)
587 /*
588  * Free a covariance matrix and all its data
589  */
590 {
591  if (c == NULL)
592  return;
593  FREE_3(c->cov);
594  FREE_3(c->cov_diag);
595  if(c->cov_sparse)
596  delete c->cov_sparse;
597  c->names.clear();
598  FREE_CMATRIX_3(c->eigen);
599  FREE_3(c->lambda);
600  FREE_3(c->inv_lambda);
601  FREE_3(c->chol);
602  FREE_3(c->ch_class);
603  if(c->proj)
604  delete c->proj;
605  if (c->sss)
606  delete c->sss;
607  c->bads.clear();
608  FREE_3(c);
609  return;
610 }
611 
612 #define EPS_3 0.05
613 
614 int mne_get_values_from_data_3 (float time, /* Interesting time point */
615  float integ, /* Time integration */
616  float **data, /* The data values (time by time) */
617  int nsamp, /* How many time points? */
618  int nch, /* How many channels */
619  float tmin, /* Time of first sample */
620  float sfreq, /* Sampling frequency */
621  int use_abs, /* Use absolute values */
622  float *value) /* The picked values */
623 /*
624  * Pick a signal value using linear interpolation
625  */
626 {
627  int n1,n2,k;
628  float s1,s2;
629  float f1,f2;
630  float sum;
631  float width;
632  int ch;
633 
634  for (ch = 0; ch < nch; ch++) {
635  /*
636  * Find out the correct samples
637  */
638  if (std::fabs(sfreq*integ) < EPS_3) { /* This is the single-sample case */
639  s1 = sfreq*(time - tmin);
640  n1 = floor(s1);
641  f1 = 1.0 + n1 - s1;
642  if (n1 < 0 || n1 > nsamp-1) {
643  printf("Sample value out of range %d (0..%d)",n1,nsamp-1);
644  return(-1);
645  }
646  /*
647  * Avoid rounding error
648  */
649  if (n1 == nsamp-1) {
650  if (std::fabs(f1 - 1.0) < 1e-3)
651  f1 = 1.0;
652  }
653  if (f1 < 1.0 && n1 > nsamp-2) {
654  printf("Sample value out of range %d (0..%d) %.4f",n1,nsamp-1,f1);
655  return(-1);
656  }
657  if (f1 < 1.0) {
658  if (use_abs)
659  sum = f1*std::fabs(data[n1][ch]) + (1.0-f1)*std::fabs(data[n1+1][ch]);
660  else
661  sum = f1*data[n1][ch] + (1.0-f1)*data[n1+1][ch];
662  }
663  else {
664  if (use_abs)
665  sum = std::fabs(data[n1][ch]);
666  else
667  sum = data[n1][ch];
668  }
669  }
670  else { /* Multiple samples */
671  s1 = sfreq*(time - 0.5*integ - tmin);
672  s2 = sfreq*(time + 0.5*integ - tmin);
673  n1 = ceil(s1); n2 = floor(s2);
674  if (n2 < n1) { /* We are within one sample interval */
675  n1 = floor(s1);
676  if (n1 < 0 || n1 > nsamp-2)
677  return (-1);
678  f1 = s1 - n1;
679  f2 = s2 - n1;
680  if (use_abs)
681  sum = 0.5*((f1+f2)*std::fabs(data[n1+1][ch]) + (2.0-f1-f2)*std::fabs(data[n1][ch]));
682  else
683  sum = 0.5*((f1+f2)*data[n1+1][ch] + (2.0-f1-f2)*data[n1][ch]);
684  }
685  else {
686  f1 = n1 - s1;
687  f2 = s2 - n2;
688  if (n1 < 0 || n1 > nsamp-1) {
689  printf("Sample value out of range %d (0..%d)",n1,nsamp-1);
690  return(-1);
691  }
692  if (n2 < 0 || n2 > nsamp-1) {
693  printf("Sample value out of range %d (0..%d)",n2,nsamp-1);
694  return(-1);
695  }
696  if (f1 != 0.0 && n1 < 1)
697  return(-1);
698  if (f2 != 0.0 && n2 > nsamp-2)
699  return(-1);
700  sum = 0.0;
701  width = 0.0;
702  if (n2 > n1) { /* Do the whole intervals */
703  if (use_abs) {
704  sum = 0.5 * std::fabs(data[n1][ch]);
705  for (k = n1+1; k < n2; k++)
706  sum = sum + std::fabs(data[k][ch]);
707  sum = sum + 0.5 * std::fabs(data[n2][ch]);
708  }
709  else {
710  sum = 0.5*data[n1][ch];
711  for (k = n1+1; k < n2; k++)
712  sum = sum + data[k][ch];
713  sum = sum + 0.5*data[n2][ch];
714  }
715  width = n2 - n1;
716  }
717  /*
718  * Add tails
719  */
720  if (use_abs) {
721  if (f1 != 0.0)
722  sum = sum + 0.5*f1*(f1*std::fabs(data[n1-1][ch]) + (2.0-f1)*std::fabs(data[n1][ch]));
723  if (f2 != 0.0)
724  sum = sum + 0.5*f2*(f2*std::fabs(data[n2+1][ch]) + (2.0-f2)*std::fabs(data[n2][ch]));
725  }
726  else {
727  if (f1 != 0.0)
728  sum = sum + 0.5*f1*(f1*data[n1-1][ch] + (2.0-f1)*data[n1][ch]);
729  if (f2 != 0.0)
730  sum = sum + 0.5*f2*(f2*data[n2+1][ch] + (2.0-f2)*data[n2][ch]);
731  }
732  width = width + f1 + f2;
733  sum = sum/width;
734  }
735  }
736  value[ch] = sum;
737  }
738  return (0);
739 }
740 
741 //============================= mne_decompose.c =============================
742 
743 int mne_decompose_eigen_3(double *mat,
744  double *lambda,
745  float **vectors, /* Eigenvectors fit into floats easily */
746  int dim)
747 /*
748  * Compute the eigenvalue decomposition of
749  * a symmetric matrix using the LAPACK routines
750  *
751  * 'mat' contains the lower triangle of the matrix
752  */
753 {
754  int np = dim*(dim+1)/2;
755  double *w = MALLOC_3(dim,double);
756  double *z = MALLOC_3(dim*dim,double);
757  double *work = MALLOC_3(3*dim,double);
758  double *dmat = MALLOC_3(np,double);
759  float *vecp = vectors[0];
760 
761  int info,k;
762  int maxi;
763  double scale;
764 
765 // maxi = idamax(&np,mat,&one);
766 // idamax workaround begin
767  maxi = 0;
768  for(int i = 0; i < np; ++i)
769  if (std::fabs(mat[i]) > std::fabs(mat[maxi]))
770  maxi = i;
771 // idamax workaround end
772 
773  scale = 1.0/mat[maxi];//scale = 1.0/mat[maxi-1];
774 
775  for (k = 0; k < np; k++)
776  dmat[k] = mat[k]*scale;
777 // dspev(compz,uplo,&dim,dmat,w,z,&dim,work,&info);
778 
779 // dspev workaround begin
780  MatrixXd dmat_tmp = MatrixXd::Zero(dim,dim);
781  int idx = 0;
782  for (int i = 0; i < dim; ++i) {
783  for(int j = 0; j <= i; ++j) {
784  dmat_tmp(i,j) = dmat[idx];
785  dmat_tmp(j,i) = dmat[idx];
786  ++idx;
787  }
788  }
789  SelfAdjointEigenSolver<MatrixXd> es;
790  es.compute(dmat_tmp);
791  for ( int i = 0; i < dim; ++i )
792  w[i] = es.eigenvalues()[i];
793 
794  idx = 0;
795  for ( int j = 0; j < dim; ++j ) {
796  for( int i = 0; i < dim; ++i ) {
797  z[idx] = es.eigenvectors()(i,j);// Column Major
798  ++idx;
799  }
800  }
801 // dspev workaround end
802 
803  info = 0;
804 
805  qDebug() << "!!!DEBUG ToDo: dspev(compz,uplo,&dim,dmat,w,z,&dim,work,&info);";
806 
807  FREE_3(work);
808  if (info != 0)
809  printf("Eigenvalue decomposition failed (LAPACK info = %d)",info);
810  else {
811  scale = 1.0/scale;
812  for (k = 0; k < dim; k++)
813  lambda[k] = scale*w[k];
814  for (k = 0; k < dim*dim; k++)
815  vecp[k] = z[k];
816  }
817  FREE_3(w);
818  FREE_3(z);
819  if (info == 0)
820  return 0;
821  else
822  return -1;
823 }
824 
825 //============================= mne_cov_matrix.c =============================
826 
827 /*
828  * Routines for handling the covariance matrices
829  */
830 
831 static int mne_lt_packed_index_3(int j, int k)
832 
833 {
834  if (j >= k)
835  return k + j*(j+1)/2;
836  else
837  return j + k*(k+1)/2;
838 }
839 
840 /*
841  * Handle the linear projection operators
842  */
843 
844 int mne_add_inv_cov_3(MneCovMatrix* c)
845 /*
846  * Calculate the inverse square roots for whitening
847  */
848 {
849  double *src = c->lambda ? c->lambda : c->cov_diag;
850  int k;
851 
852  if (src == NULL) {
853  qCritical("Covariance matrix is not diagonal or not decomposed.");
854  return FAIL;
855  }
856  c->inv_lambda = REALLOC_3(c->inv_lambda,c->ncov,double);
857  for (k = 0; k < c->ncov; k++) {
858  if (src[k] <= 0.0)
859  c->inv_lambda[k] = 0.0;
860  else
861  c->inv_lambda[k] = 1.0/sqrt(src[k]);
862  }
863  return OK;
864 }
865 
866 static int condition_cov_3(MneCovMatrix* c, float rank_threshold, int use_rank)
867 
868 {
869  double *scale = NULL;
870  double *cov = NULL;
871  double *lambda = NULL;
872  float **eigen = NULL;
873  double **data1 = NULL;
874  double **data2 = NULL;
875  double magscale,gradscale,eegscale;
876  int nmag,ngrad,neeg,nok;
877  int j,k;
878  int res = FAIL;
879 
880  if (c->cov_diag)
881  return OK;
882  if (!c->ch_class) {
883  qCritical("Channels not classified. Rank cannot be determined.");
884  return FAIL;
885  }
886  magscale = gradscale = eegscale = 0.0;
887  nmag = ngrad = neeg = 0;
888  for (k = 0; k < c->ncov; k++) {
889  if (c->ch_class[k] == MNE_COV_CH_MEG_MAG) {
890  magscale += c->cov[mne_lt_packed_index_3(k,k)]; nmag++;
891  }
892  else if (c->ch_class[k] == MNE_COV_CH_MEG_GRAD) {
893  gradscale += c->cov[mne_lt_packed_index_3(k,k)]; ngrad++;
894  }
895  else if (c->ch_class[k] == MNE_COV_CH_EEG) {
896  eegscale += c->cov[mne_lt_packed_index_3(k,k)]; neeg++;
897  }
898 #ifdef DEBUG
899  fprintf(stdout,"%d ",c->ch_class[k]);
900 #endif
901  }
902 #ifdef DEBUG
903  fprintf(stdout,"\n");
904 #endif
905  if (nmag > 0)
906  magscale = magscale > 0.0 ? sqrt(nmag/magscale) : 0.0;
907  if (ngrad > 0)
908  gradscale = gradscale > 0.0 ? sqrt(ngrad/gradscale) : 0.0;
909  if (neeg > 0)
910  eegscale = eegscale > 0.0 ? sqrt(neeg/eegscale) : 0.0;
911 #ifdef DEBUG
912  fprintf(stdout,"%d %g\n",nmag,magscale);
913  fprintf(stdout,"%d %g\n",ngrad,gradscale);
914  fprintf(stdout,"%d %g\n",neeg,eegscale);
915 #endif
916  scale = MALLOC_3(c->ncov,double);
917  for (k = 0; k < c->ncov; k++) {
918  if (c->ch_class[k] == MNE_COV_CH_MEG_MAG)
919  scale[k] = magscale;
920  else if (c->ch_class[k] == MNE_COV_CH_MEG_GRAD)
921  scale[k] = gradscale;
922  else if (c->ch_class[k] == MNE_COV_CH_EEG)
923  scale[k] = eegscale;
924  else
925  scale[k] = 1.0;
926  }
927  cov = MALLOC_3(c->ncov*(c->ncov+1)/2.0,double);
928  lambda = MALLOC_3(c->ncov,double);
929  eigen = ALLOC_CMATRIX_3(c->ncov,c->ncov);
930  for (j = 0; j < c->ncov; j++)
931  for (k = 0; k <= j; k++)
932  cov[mne_lt_packed_index_3(j,k)] = c->cov[mne_lt_packed_index_3(j,k)]*scale[j]*scale[k];
933  if (mne_decompose_eigen_3(cov,lambda,eigen,c->ncov) == 0) {
934 #ifdef DEBUG
935  for (k = 0; k < c->ncov; k++)
936  fprintf(stdout,"%g ",lambda[k]/lambda[c->ncov-1]);
937  fprintf(stdout,"\n");
938 #endif
939  nok = 0;
940  for (k = c->ncov-1; k >= 0; k--) {
941  if (lambda[k] >= rank_threshold*lambda[c->ncov-1])
942  nok++;
943  else
944  break;
945  }
946  printf("\n\tEstimated covariance matrix rank = %d (%g)\n",nok,lambda[c->ncov-nok]/lambda[c->ncov-1]);
947  if (use_rank > 0 && use_rank < nok) {
948  nok = use_rank;
949  printf("\tUser-selected covariance matrix rank = %d (%g)\n",nok,lambda[c->ncov-nok]/lambda[c->ncov-1]);
950  }
951  /*
952  * Put it back together
953  */
954  for (j = 0; j < c->ncov-nok; j++)
955  lambda[j] = 0.0;
956  data1 = ALLOC_DCMATRIX_3(c->ncov,c->ncov);
957  for (j = 0; j < c->ncov; j++) {
958 #ifdef DEBUG
959  mne_print_vector(stdout,NULL,eigen[j],c->ncov);
960 #endif
961  for (k = 0; k < c->ncov; k++)
962  data1[j][k] = sqrt(lambda[j])*eigen[j][k];
963  }
964  data2 = mne_dmatt_dmat_mult2_3(data1,data1,c->ncov,c->ncov,c->ncov);
965 #ifdef DEBUG
966  printf(">>>\n");
967  for (j = 0; j < c->ncov; j++)
968  mne_print_dvector(stdout,NULL,data2[j],c->ncov);
969  printf(">>>\n");
970 #endif
971  /*
972  * Scale back
973  */
974  for (k = 0; k < c->ncov; k++)
975  if (scale[k] > 0.0)
976  scale[k] = 1.0/scale[k];
977  for (j = 0; j < c->ncov; j++)
978  for (k = 0; k <= j; k++)
979  if (c->cov[mne_lt_packed_index_3(j,k)] != 0.0)
980  c->cov[mne_lt_packed_index_3(j,k)] = scale[j]*scale[k]*data2[j][k];
981  res = nok;
982  }
983  FREE_3(cov);
984  FREE_3(scale);
985  FREE_3(lambda);
986  FREE_CMATRIX_3(eigen);
987  FREE_DCMATRIX_3(data1);
988  FREE_DCMATRIX_3(data2);
989  return res;
990 }
991 
992 static int check_cov_data(double *vals, int nval)
993 
994 {
995  int k;
996  double sum;
997 
998  for (k = 0, sum = 0.0; k < nval; k++)
999  sum += vals[k];
1000  if (sum == 0.0) {
1001  qCritical("Sum of covariance matrix elements is zero!");
1002  return FAIL;
1003  }
1004  return OK;
1005 }
1006 
1007 int mne_classify_channels_cov(MneCovMatrix* cov,
1008  const QList<FiffChInfo>& chs,
1009  int nchan)
1010 /*
1011  * Assign channel classes in a covariance matrix with help of channel infos
1012  */
1013 {
1014  int k,p;
1015  FiffChInfo ch;
1016 
1017  if (chs.isEmpty()) {
1018  qCritical("Channel information not available in mne_classify_channels_cov");
1019  goto bad;
1020  }
1021  cov->ch_class = REALLOC_3(cov->ch_class,cov->ncov,int);
1022  for (k = 0; k < cov->ncov; k++) {
1023  cov->ch_class[k] = MNE_COV_CH_UNKNOWN;
1024  for (p = 0; p < nchan; p++) {
1025  if (QString::compare(chs[p].ch_name,cov->names[k]) == 0) {
1026  ch = chs[p];
1027  if (ch.kind == FIFFV_MEG_CH) {
1028  if (ch.unit == FIFF_UNIT_T)
1029  cov->ch_class[k] = MNE_COV_CH_MEG_MAG;
1030  else
1031  cov->ch_class[k] = MNE_COV_CH_MEG_GRAD;
1032  }
1033  else if (ch.kind == FIFFV_EEG_CH)
1034  cov->ch_class[k] = MNE_COV_CH_EEG;
1035  break;
1036  }
1037  }
1038 // if (!ch) {
1039 // printf("Could not find channel info for channel %s in mne_classify_channels_cov",cov->names[k].toUtf8().constData());
1040 // goto bad;
1041 // }
1042  }
1043  return OK;
1044 
1045 bad : {
1046  FREE_3(cov->ch_class);
1047  cov->ch_class = NULL;
1048  return FAIL;
1049  }
1050 }
1051 
1052 static int mne_decompose_eigen_cov_small_3(MneCovMatrix* c,float small, int use_rank)
1053 /*
1054  * Do the eigenvalue decomposition
1055  */
1056 {
1057  int np,k,p,rank;
1058  float rank_threshold = 1e-6;
1059 
1060  if (small < 0)
1061  small = 1.0;
1062 
1063  if (!c)
1064  return OK;
1065  if (c->cov_diag)
1066  return mne_add_inv_cov_3(c);
1067  if (c->lambda && c->eigen) {
1068  printf("\n\tEigenvalue decomposition had been precomputed.\n");
1069  c->nzero = 0;
1070  for (k = 0; k < c->ncov; k++, c->nzero++)
1071  if (c->lambda[k] > 0)
1072  break;
1073  }
1074  else {
1075  FREE_3(c->lambda); c->lambda = NULL;
1076  FREE_CMATRIX_3(c->eigen); c->eigen = NULL;
1077 
1078  if ((rank = condition_cov_3(c,rank_threshold,use_rank)) < 0)
1079  return FAIL;
1080 
1081  np = c->ncov*(c->ncov+1)/2;
1082  c->lambda = MALLOC_3(c->ncov,double);
1083  c->eigen = ALLOC_CMATRIX_3(c->ncov,c->ncov);
1084  if (mne_decompose_eigen_3(c->cov,c->lambda,c->eigen,c->ncov) != 0)
1085  goto bad;
1086  c->nzero = c->ncov - rank;
1087  for (k = 0; k < c->nzero; k++)
1088  c->lambda[k] = 0.0;
1089  /*
1090  * Find which eigenvectors correspond to EEG/MEG
1091  */
1092  {
1093  float meglike,eeglike;
1094  int nmeg,neeg;
1095 
1096  nmeg = neeg = 0;
1097  for (k = c->nzero; k < c->ncov; k++) {
1098  meglike = eeglike = 0.0;
1099  for (p = 0; p < c->ncov; p++) {
1100  if (c->ch_class[p] == MNE_COV_CH_EEG)
1101  eeglike += std::fabs(c->eigen[k][p]);
1102  else if (c->ch_class[p] == MNE_COV_CH_MEG_MAG || c->ch_class[p] == MNE_COV_CH_MEG_GRAD)
1103  meglike += std::fabs(c->eigen[k][p]);
1104  }
1105  if (meglike > eeglike)
1106  nmeg++;
1107  else
1108  neeg++;
1109  }
1110  printf("\t%d MEG and %d EEG-like channels remain in the whitened data\n",nmeg,neeg);
1111  }
1112  }
1113  return mne_add_inv_cov_3(c);
1114 
1115 bad : {
1116  FREE_3(c->lambda); c->lambda = NULL;
1117  FREE_CMATRIX_3(c->eigen); c->eigen = NULL;
1118  return FAIL;
1119  }
1120 }
1121 
1122 int mne_decompose_eigen_cov_3(MneCovMatrix* c)
1123 
1124 {
1125  return mne_decompose_eigen_cov_small_3(c,-1.0,-1);
1126 }
1127 
1128 //============================= mne_whiten.c =============================
1129 
1130 int mne_whiten_data(float **data, float **whitened_data, int np, int nchan, MneCovMatrix* C)
1131 /*
1132  * Apply the whitening operation
1133  */
1134 {
1135  int j,k;
1136  float *one = NULL,*orig,*white;
1137  double *inv;
1138 
1139  if (data == NULL || np <= 0)
1140  return OK;
1141 
1142  if (C->ncov != nchan) {
1143  printf("Incompatible covariance matrix. Cannot whiten the data.");
1144  return FAIL;
1145  }
1146  inv = C->inv_lambda;
1147  if (mne_is_diag_cov_3(C)) {
1148  // printf("<DEBUG> Performing Diag\n");
1149  for (j = 0; j < np; j++) {
1150  orig = data[j];
1151  white = whitened_data[j];
1152  for (k = 0; k < nchan; k++)
1153  white[k] = orig[k]*inv[k];
1154  }
1155  }
1156  else {
1157  /*
1158  * This is arranged so that whitened_data can be the same matrix as the original
1159  */
1160  one = MALLOC_3(nchan,float);
1161  for (j = 0; j < np; j++) {
1162  orig = data[j];
1163  white = whitened_data[j];
1164  for (k = C->nzero; k < nchan; k++)
1165  one[k] = mne_dot_vectors_3(C->eigen[k],orig,nchan);
1166  for (k = 0; k < C->nzero; k++)
1167  white[k] = 0.0;
1168  for (k = C->nzero; k < nchan; k++)
1169  white[k] = one[k]*inv[k];
1170  }
1171  FREE_3(one);
1172  }
1173  return OK;
1174 }
1175 
1176 int mne_whiten_one_data(float *data, float *whitened_data, int nchan, MneCovMatrix* C)
1177 
1178 {
1179  float *datap[1];
1180  float *whitened_datap[1];
1181 
1182  datap[0] = data;
1183  whitened_datap[0] = whitened_data;
1184 
1185  return mne_whiten_data(datap,whitened_datap,1,nchan,C);
1186 }
1187 
1188 //=============================================================================================================
1189 //============================= dipole_fit_setup.c =============================
1190 static void free_dipole_fit_funcs(dipoleFitFuncs f)
1191 
1192 {
1193  if (!f)
1194  return;
1195 
1196  if (f->meg_client_free && f->meg_client)
1197  f->meg_client_free(f->meg_client);
1198  if (f->eeg_client_free && f->eeg_client)
1199  f->eeg_client_free(f->eeg_client);
1200 
1201  FREE_3(f);
1202  return;
1203 }
1204 
1205 //static void regularize_cov(MneCovMatrix* c, /* The matrix to regularize */
1206 // float *regs, /* Regularization values to apply (fractions of the
1207 // * average diagonal values for each class */
1208 // int *active) /* Which of the channels are 'active' */
1210 // * Regularize different parts of the noise covariance matrix differently
1211 // */
1212 //{
1213 // int j;
1214 // float sums[3],nn[3];
1215 // int nkind = 3;
1216 
1217 // if (!c->cov || !c->ch_class)
1218 // return;
1219 
1220 // for (j = 0; j < nkind; j++) {
1221 // sums[j] = 0.0;
1222 // nn[j] = 0;
1223 // }
1224 // /*
1225 // * Compute the averages over the diagonal elements for each class
1226 // */
1227 // for (j = 0; j < c->ncov; j++) {
1228 // if (c->ch_class[j] >= 0) {
1229 // if (!active || active[j]) {
1230 // sums[c->ch_class[j]] += c->cov[mne_lt_packed_index_3(j,j)];
1231 // nn[c->ch_class[j]]++;
1232 // }
1233 // }
1234 // }
1235 // printf("Average noise-covariance matrix diagonals:\n");
1236 // for (j = 0; j < nkind; j++) {
1237 // if (nn[j] > 0) {
1238 // sums[j] = sums[j]/nn[j];
1239 // if (j == MNE_COV_CH_MEG_MAG)
1240 // printf("\tMagnetometers : %-7.2f fT reg = %-6.2f\n",1e15*sqrt(sums[j]),regs[j]);
1241 // else if (j == MNE_COV_CH_MEG_GRAD)
1242 // printf("\tPlanar gradiometers : %-7.2f fT/cm reg = %-6.2f\n",1e13*sqrt(sums[j]),regs[j]);
1243 // else
1244 // printf("\tEEG : %-7.2f uV reg = %-6.2f\n",1e6*sqrt(sums[j]),regs[j]);
1245 // sums[j] = regs[j]*sums[j];
1246 // }
1247 // }
1248 // /*
1249 // * Add thee proper amount to the diagonal
1250 // */
1251 // for (j = 0; j < c->ncov; j++)
1252 // if (c->ch_class[j] >= 0)
1253 // c->cov[mne_lt_packed_index_3(j,j)] += sums[c->ch_class[j]];
1254 
1255 // printf("Noise-covariance regularized as requested.\n");
1256 // return;
1257 //}
1258 
1259 void mne_regularize_cov(MneCovMatrix* c, /* The matrix to regularize */
1260  float *regs) /* Regularization values to apply (fractions of the
1261  * average diagonal values for each class */
1262 /*
1263  * Regularize different parts of the noise covariance matrix differently
1264  */
1265 {
1266  int j;
1267  float sums[3],nn[3];
1268  int nkind = 3;
1269 
1270  if (!c->cov || !c->ch_class)
1271  return;
1272 
1273  for (j = 0; j < nkind; j++) {
1274  sums[j] = 0.0;
1275  nn[j] = 0;
1276  }
1277  /*
1278  * Compute the averages over the diagonal elements for each class
1279  */
1280  for (j = 0; j < c->ncov; j++) {
1281  if (c->ch_class[j] >= 0) {
1282  sums[c->ch_class[j]] += c->cov[mne_lt_packed_index_3(j,j)];
1283  nn[c->ch_class[j]]++;
1284  }
1285  }
1286  printf("Average noise-covariance matrix diagonals:\n");
1287  for (j = 0; j < nkind; j++) {
1288  if (nn[j] > 0) {
1289  sums[j] = sums[j]/nn[j];
1290  if (j == MNE_COV_CH_MEG_MAG)
1291  printf("\tMagnetometers : %-7.2f fT reg = %-6.2f\n",1e15*sqrt(sums[j]),regs[j]);
1292  else if (j == MNE_COV_CH_MEG_GRAD)
1293  printf("\tPlanar gradiometers : %-7.2f fT/cm reg = %-6.2f\n",1e13*sqrt(sums[j]),regs[j]);
1294  else
1295  printf("\tEEG : %-7.2f uV reg = %-6.2f\n",1e6*sqrt(sums[j]),regs[j]);
1296  sums[j] = regs[j]*sums[j];
1297  }
1298  }
1299  /*
1300  * Add thee proper amount to the diagonal
1301  */
1302  for (j = 0; j < c->ncov; j++)
1303  if (c->ch_class[j] >= 0)
1304  c->cov[mne_lt_packed_index_3(j,j)] += sums[c->ch_class[j]];
1305 
1306  printf("Noise-covariance regularized as requested.\n");
1307  return;
1308 }
1309 
1310 //============================= dipole_fit_setup.c =============================
1311 
1312 static dipoleFitFuncs new_dipole_fit_funcs()
1313 
1314 {
1315  dipoleFitFuncs f = MALLOC_3(1,dipoleFitFuncsRec);
1316 
1317  f->meg_field = NULL;
1318  f->eeg_pot = NULL;
1319  f->meg_vec_field = NULL;
1320  f->eeg_vec_pot = NULL;
1321  f->meg_client = NULL;
1322  f->meg_client_free = NULL;
1323  f->eeg_client = NULL;
1324  f->eeg_client_free = NULL;
1325 
1326  return f;
1327 }
1328 
1329 //============================= mne_simplex_fit.c =============================
1330 
1331 /*
1332  * This routine comes from Numerical recipes
1333  */
1334 
1335 #define ALPHA 1.0
1336 #define BETA 0.5
1337 #define GAMMA 2.0
1338 
1339 static float tryit (float **p,
1340  float *y,
1341  float *psum,
1342  int ndim,
1343  float (*func)(float *x,int npar,void *user_data), /* The function to be evaluated */
1344  void *user_data, /* Data to be passed to the above function in each evaluation */
1345  int ihi,
1346  int *neval,
1347  float fac)
1348 
1349 {
1350  int j;
1351  float fac1,fac2,ytry,*ptry;
1352 
1353  ptry = ALLOC_FLOAT_3(ndim);
1354  fac1 = (1.0-fac)/ndim;
1355  fac2 = fac1-fac;
1356  for (j = 0; j < ndim; j++)
1357  ptry[j] = psum[j]*fac1-p[ihi][j]*fac2;
1358  ytry = (*func)(ptry,ndim,user_data);
1359  ++(*neval);
1360  if (ytry < y[ihi]) {
1361  y[ihi] = ytry;
1362  for (j = 0; j < ndim; j++) {
1363  psum[j] += ptry[j]-p[ihi][j];
1364  p[ihi][j] = ptry[j];
1365  }
1366  }
1367  FREE_3(ptry);
1368  return ytry;
1369 }
1370 
1371 int mne_simplex_minimize(float **p, /* The initial simplex */
1372  float *y, /* Function values at the vertices */
1373  int ndim, /* Number of variables */
1374  float ftol, /* Relative convergence tolerance */
1375  float (*func)(float *x,int npar,void *user_data), /* The function to be evaluated */
1376  void *user_data, /* Data to be passed to the above function in each evaluation */
1377  int max_eval, /* Maximum number of function evaluations */
1378  int *neval, /* Number of function evaluations */
1379  int report, /* How often to report (-1 = no_reporting) */
1380  int (*report_func)(int loop,
1381  float *fitpar, int npar,
1382  double fval)) /* The function to be called when reporting */
1383 
1384 /*
1385  * Minimization with the simplex algorithm
1386  * Modified from Numerical recipes
1387  */
1388 
1389 {
1390  int i,j,ilo,ihi,inhi;
1391  int mpts = ndim+1;
1392  float ytry,ysave,sum,rtol,*psum;
1393  int result = 0;
1394  int count = 0;
1395  int loop = 1;
1396 
1397  psum = ALLOC_FLOAT_3(ndim);
1398  *neval = 0;
1399  for (j = 0; j < ndim; j++) {
1400  for (i = 0,sum = 0.0; i<mpts; i++)
1401  sum += p[i][j];
1402  psum[j] = sum;
1403  }
1404  if (report_func != NULL && report > 0)
1405  (void)report_func (0,p[0],ndim,-1.0);
1406 
1407  for (;;count++,loop++) {
1408  ilo = 1;
1409  ihi = y[1]>y[2] ? (inhi = 2,1) : (inhi = 1,2);
1410  for (i = 0; i < mpts; i++) {
1411  if (y[i] < y[ilo]) ilo = i;
1412  if (y[i] > y[ihi]) {
1413  inhi = ihi;
1414  ihi = i;
1415  } else if (y[i] > y[inhi])
1416  if (i != ihi) inhi = i;
1417  }
1418  rtol = 2.0*std::fabs(y[ihi]-y[ilo])/(std::fabs(y[ihi])+std::fabs(y[ilo]));
1419  /*
1420  * Report that we are proceeding...
1421  */
1422  if (count == report && report_func != NULL) {
1423  if (report_func (loop,p[ilo],ndim,y[ilo])) {
1424  qWarning("Interation interrupted.");
1425  result = -1;
1426  break;
1427  }
1428  count = 0;
1429  }
1430  if (rtol < ftol) break;
1431  if (*neval >= max_eval) {
1432  qWarning("Maximum number of evaluations exceeded.");
1433  result = -1;
1434  break;
1435  }
1436  ytry = tryit(p,y,psum,ndim,func,user_data,ihi,neval,-ALPHA);
1437  if (ytry <= y[ilo])
1438  ytry = tryit(p,y,psum,ndim,func,user_data,ihi,neval,GAMMA);
1439  else if (ytry >= y[inhi]) {
1440  ysave = y[ihi];
1441  ytry = tryit(p,y,psum,ndim,func,user_data,ihi,neval,BETA);
1442  if (ytry >= ysave) {
1443  for (i = 0; i < mpts; i++) {
1444  if (i != ilo) {
1445  for (j = 0; j < ndim; j++) {
1446  psum[j] = 0.5*(p[i][j]+p[ilo][j]);
1447  p[i][j] = psum[j];
1448  }
1449  y[i] = (*func)(psum,ndim,user_data);
1450  }
1451  }
1452  *neval += ndim;
1453  for (j = 0; j < ndim; j++) {
1454  for (i = 0,sum = 0.0; i < mpts; i++)
1455  sum += p[i][j];
1456  psum[j] = sum;
1457  }
1458  }
1459  }
1460  }
1461  FREE_3(psum);
1462  return (result);
1463 }
1464 
1465 #undef ALPHA
1466 #undef BETA
1467 #undef GAMMA
1468 
1469 //============================= fit_sphere.c =============================
1470 
1471 namespace MNELIB
1472 {
1473 
1474 typedef struct {
1475  float **rr;
1476  int np;
1477  int report;
1479 
1480 }
1481 
1482 static int report_func(int loop,
1483  float *fitpar,
1484  int npar,
1485  double fval)
1486 /*
1487  * Report periodically
1488  */
1489 {
1490  float *r0 = fitpar;
1491  (void) npar; // avoid unused variable warning.
1492 
1493  printf("loop %d r0 %7.1f %7.1f %7.1f fval %g\n",
1494  loop,1000*r0[0],1000*r0[1],1000*r0[2],fval);
1495 
1496  return OK;
1497 }
1498 
1499 static float fit_sphere_eval(float *fitpar,
1500  int npar,
1501  void *user_data)
1502 /*
1503  * Calculate the cost function value
1504  * Optimize for the radius inside here
1505  */
1506 {
1507  fitSphereUser user = (fitSphereUser)user_data;
1508  float *r0 = fitpar;
1509  float diff[3];
1510  int k;
1511  float sum,sum2,one,F;
1512  (void) npar; // avoid unused variable warning.
1513 
1514  for (k = 0, sum = sum2 = 0.0; k < user->np; k++) {
1515  VEC_DIFF_3(r0,user->rr[k],diff);
1516  one = VEC_LEN_3(diff);
1517  sum += one;
1518  sum2 += one*one;
1519  }
1520  F = sum2 - sum*sum/user->np;
1521 
1522  if (user->report)
1523  printf("r0 %7.1f %7.1f %7.1f R %7.1f fval %g\n",
1524  1000*r0[0],1000*r0[1],1000*r0[2],1000*sum/user->np,F);
1525 
1526  return F;
1527 }
1528 
1529 static float opt_rad(float *r0,fitSphereUser user)
1530 
1531 {
1532  float sum, diff[3], one;
1533  int k;
1534 
1535  for (k = 0, sum = 0.0; k < user->np; k++) {
1536  VEC_DIFF_3(r0,user->rr[k],diff);
1537  one = VEC_LEN_3(diff);
1538  sum += one;
1539  }
1540  return sum/user->np;
1541 }
1542 
1543 static void calculate_cm_ave_dist(float **rr, int np, float *cm, float *avep)
1544 
1545 {
1546  int k,q;
1547  float ave,diff[3];
1548 
1549  for (q = 0; q < 3; q++)
1550  cm[q] = 0.0;
1551 
1552  for (k = 0; k < np; k++)
1553  for (q = 0; q < 3; q++)
1554  cm[q] += rr[k][q];
1555 
1556  if (np > 0) {
1557  for (q = 0; q < 3; q++)
1558  cm[q] = cm[q]/np;
1559 
1560  for (k = 0, ave = 0.0; k < np; k++) {
1561  for (q = 0; q < 3; q++)
1562  diff[q] = rr[k][q] - cm[q];
1563  ave += VEC_LEN_3(diff);
1564  }
1565  *avep = ave/np;
1566  }
1567  return;
1568 }
1569 
1570 static float **make_initial_simplex(float *pars,
1571  int npar,
1572  float size)
1573 /*
1574  * Make the initial tetrahedron
1575  */
1576 {
1577  float **simplex = ALLOC_CMATRIX_3(npar+1,npar);
1578  int k;
1579 
1580  for (k = 0; k < npar+1; k++)
1581  memcpy (simplex[k],pars,npar*sizeof(float));
1582 
1583  for (k = 1; k < npar+1; k++)
1584  simplex[k][k-1] = simplex[k][k-1] + size;
1585  return (simplex);
1586 }
1587 
1588 int fit_sphere_to_points(float **rr,
1589  int np,
1590  float simplex_size,
1591  float *r0,
1592  float *R)
1593 /*
1594  * Find the optimal sphere origin
1595  */
1596 {
1597  fitSphereUserRec user;
1598  float ftol = 1e-5;
1599  int max_eval = 500;
1600  int report_interval = -1;
1601  int neval = 0;
1602  float **init_simplex = NULL;
1603  float *init_vals = NULL;
1604 
1605  float cm[3],R0;
1606  int k;
1607 
1608  int res = FAIL;
1609 
1610  user.rr = rr;
1611  user.np = np;
1612 
1613  calculate_cm_ave_dist(rr,np,cm,&R0);
1614 
1615 #ifdef DEBUG
1616  printf("cm %7.1f %7.1f %7.1f R %7.1f\n",
1617  1000*cm[0],1000*cm[1],1000*cm[2],1000*R0);
1618 #endif
1619 
1620  init_simplex = make_initial_simplex(cm,3,simplex_size);
1621 
1622  init_vals = MALLOC_3(4,float);
1623 
1624 #ifdef DEBUG
1625  user.report = TRUE;
1626 #else
1627  user.report = FALSE;
1628 #endif
1629 
1630  for (k = 0; k < 4; k++)
1631  init_vals[k] = fit_sphere_eval(init_simplex[k],3,&user);
1632 
1633  if (mne_simplex_minimize(init_simplex, /* The initial simplex */
1634  init_vals, /* Function values at the vertices */
1635  3, /* Number of variables */
1636  ftol, /* Relative convergence tolerance */
1637  fit_sphere_eval, /* The function to be evaluated */
1638  &user, /* Data to be passed to the above function in each evaluation */
1639  max_eval, /* Maximum number of function evaluations */
1640  &neval, /* Number of function evaluations */
1641  report_interval, /* How often to report (-1 = no_reporting) */
1642  report_func) != OK) /* The function to be called when reporting */
1643  goto out;
1644 
1645  r0[X_3] = init_simplex[0][X_3];
1646  r0[Y_3] = init_simplex[0][Y_3];
1647  r0[Z_3] = init_simplex[0][Z_3];
1648  *R = opt_rad(r0,&user);
1649 
1650  res = OK;
1651  goto out;
1652 
1653 out : {
1654  FREE_3(init_vals);
1655  FREE_CMATRIX_3(init_simplex);
1656  return res;
1657  }
1658 }
1659 
1660 //============================= mne_lin_proj.c =============================
1661 
1662 void mne_proj_op_report_data_3(FILE *out,const char *tag, MneProjOp* op, int list_data,
1663  char **exclude, int nexclude)
1664 /*
1665  * Output info about the projection operator
1666  */
1667 {
1668  int j,k,p,q;
1669  MneProjItem* it;
1670  MneNamedMatrix* vecs;
1671  int found;
1672 
1673  if (out == NULL)
1674  return;
1675  if (op == NULL)
1676  return;
1677  if (op->nitems <= 0) {
1678  fprintf(out,"Empty operator\n");
1679  return;
1680  }
1681 
1682  for (k = 0; k < op->nitems; k++) {
1683  it = op->items[k];
1684  if (list_data && tag)
1685  fprintf(out,"%s\n",tag);
1686  if (tag)
1687  fprintf(out,"%s",tag);
1688  fprintf(out,"# %d : %s : %d vecs : %d chs %s %s\n",
1689  k+1,it->desc.toUtf8().constData(),it->nvec,it->vecs->ncol,
1690  it->has_meg ? "MEG" : "EEG",
1691  it->active ? "active" : "idle");
1692  if (list_data && tag)
1693  fprintf(out,"%s\n",tag);
1694  if (list_data) {
1695  vecs = op->items[k]->vecs;
1696 
1697  for (q = 0; q < vecs->ncol; q++) {
1698  fprintf(out,"%-10s",vecs->collist[q].toUtf8().constData());
1699  fprintf(out,q < vecs->ncol-1 ? " " : "\n");
1700  }
1701  for (p = 0; p < vecs->nrow; p++)
1702  for (q = 0; q < vecs->ncol; q++) {
1703  for (j = 0, found = 0; j < nexclude; j++) {
1704  if (QString::compare(exclude[j],vecs->collist[q]) == 0) {
1705  found = 1;
1706  break;
1707  }
1708  }
1709  fprintf(out,"%10.5g ",found ? 0.0 : vecs->data[p][q]);
1710  fprintf(out,q < vecs->ncol-1 ? " " : "\n");
1711  }
1712  if (list_data && tag)
1713  fprintf(out,"%s\n",tag);
1714  }
1715  }
1716  return;
1717 }
1718 
1719 void mne_proj_op_report_3(FILE *out,const char *tag, MneProjOp* op)
1720 {
1721  mne_proj_op_report_data_3(out,tag,op, FALSE, NULL, 0);
1722 }
1723 
1724 //============================= mne_named_vector.c =============================
1725 
1726 int mne_pick_from_named_vector_3(mneNamedVector vec, const QStringList& names, int nnames, int require_all, float *res)
1727 /*
1728  * Pick the desired elements from the named vector
1729  */
1730 {
1731  int found;
1732  int k,p;
1733 
1734  if (vec->names.size() == 0) {
1735  qCritical("No names present in vector. Cannot pick.");
1736  return FAIL;
1737  }
1738 
1739  for (k = 0; k < nnames; k++)
1740  res[k] = 0.0;
1741 
1742  for (k = 0; k < nnames; k++) {
1743  found = 0;
1744  for (p = 0; p < vec->nvec; p++) {
1745  if (QString::compare(vec->names[p],names[k]) == 0) {
1746  res[k] = vec->data[p];
1747  found = TRUE;
1748  break;
1749  }
1750  }
1751  if (!found && require_all) {
1752  qCritical("All required elements not found in named vector.");
1753  return FAIL;
1754  }
1755  }
1756  return OK;
1757 }
1758 
1759 //============================= mne_lin_proj_io.c =============================
1760 
1761 MneProjOp* mne_read_proj_op_from_node_3(//fiffFile in,
1762  FiffStream::SPtr& stream,
1763  const FiffDirNode::SPtr& start)
1764 /*
1765  * Load all the linear projection data
1766  */
1767 {
1768  MneProjOp* op = NULL;
1769  QList<FiffDirNode::SPtr> proj;
1770  FiffDirNode::SPtr start_node;
1771  QList<FiffDirNode::SPtr> items;
1772  FiffDirNode::SPtr node;
1773  int k;
1774  QString item_desc, desc_tag;
1775  int global_nchan,item_nchan,nlist;
1776  QStringList item_names;
1777  int item_kind;
1778  float **item_vectors = NULL;
1779  int item_nvec;
1780  int item_active;
1781  MneNamedMatrix* item;
1782  FiffTag::SPtr t_pTag;
1783  QStringList emptyList;
1784  int pos;
1785 
1786  if (!stream) {
1787  qCritical("File not open mne_read_proj_op_from_node");
1788  goto bad;
1789  }
1790 
1791  if (!start || start->isEmpty())
1792  start_node = stream->dirtree();
1793  else
1794  start_node = start;
1795 
1796  op = new MneProjOp();
1797  proj = start_node->dir_tree_find(FIFFB_PROJ);
1798  if (proj.size() == 0 || proj[0]->isEmpty()) /* The caller must recognize an empty projection */
1799  goto out;
1800  /*
1801  * Only the first projection block is recognized
1802  */
1803  items = proj[0]->dir_tree_find(FIFFB_PROJ_ITEM);
1804  if (items.size() == 0 || items[0]->isEmpty()) /* The caller must recognize an empty projection */
1805  goto out;
1806  /*
1807  * Get a common number of channels
1808  */
1809  node = proj[0];
1810  if(!node->find_tag(stream, FIFF_NCHAN, t_pTag))
1811  global_nchan = 0;
1812  else {
1813  global_nchan = *t_pTag->toInt();
1814 // TAG_FREE_3(tag);
1815  }
1816  /*
1817  * Proceess each item
1818  */
1819  for (k = 0; k < items.size(); k++) {
1820  node = items[k];
1821  /*
1822  * Complicated procedure for getting the description
1823  */
1824  item_desc.clear();
1825 
1826  if (node->find_tag(stream, FIFF_NAME, t_pTag)) {
1827  item_desc += t_pTag->toString();
1828  }
1829 
1830  /*
1831  * Take the first line of description if it exists
1832  */
1833  if (node->find_tag(stream, FIFF_DESCRIPTION, t_pTag)) {
1834  desc_tag = t_pTag->toString();
1835  if((pos = desc_tag.indexOf("\n")) >= 0)
1836  desc_tag.truncate(pos);
1837  if (!item_desc.isEmpty())
1838  item_desc += " ";
1839  item_desc += desc_tag;
1840  }
1841  /*
1842  * Possibility to override number of channels here
1843  */
1844  if (!node->find_tag(stream, FIFF_NCHAN, t_pTag)) {
1845  item_nchan = global_nchan;
1846  }
1847  else {
1848  item_nchan = *t_pTag->toInt();
1849  }
1850  if (item_nchan <= 0) {
1851  qCritical("Number of channels incorrectly specified for one of the projection items.");
1852  goto bad;
1853  }
1854  /*
1855  * Take care of the channel names
1856  */
1857  if (!node->find_tag(stream, FIFF_PROJ_ITEM_CH_NAME_LIST, t_pTag))
1858  goto bad;
1859 
1860  mne_string_to_name_list_3(t_pTag->toString(),item_names,nlist);
1861  if (nlist != item_nchan) {
1862  printf("Channel name list incorrectly specified for proj item # %d",k+1);
1863  item_names.clear();
1864  goto bad;
1865  }
1866  /*
1867  * Kind of item
1868  */
1869  if (!node->find_tag(stream, FIFF_PROJ_ITEM_KIND, t_pTag))
1870  goto bad;
1871  item_kind = *t_pTag->toInt();
1872  /*
1873  * How many vectors
1874  */
1875  if (!node->find_tag(stream,FIFF_PROJ_ITEM_NVEC, t_pTag))
1876  goto bad;
1877  item_nvec = *t_pTag->toInt();
1878  /*
1879  * The projection data
1880  */
1881  if (!node->find_tag(stream,FIFF_PROJ_ITEM_VECTORS, t_pTag))
1882  goto bad;
1883 
1884  MatrixXf tmp_item_vectors = t_pTag->toFloatMatrix().transpose();
1885  item_vectors = ALLOC_CMATRIX_3(tmp_item_vectors.rows(),tmp_item_vectors.cols());
1886  fromFloatEigenMatrix_3(tmp_item_vectors, item_vectors);
1887 
1888  /*
1889  * Is this item active?
1890  */
1891  if (node->find_tag(stream, FIFF_MNE_PROJ_ITEM_ACTIVE, t_pTag)) {
1892  item_active = *t_pTag->toInt();
1893  }
1894  else
1895  item_active = FALSE;
1896  /*
1897  * Ready to add
1898  */
1899  item = MneNamedMatrix::build_named_matrix(item_nvec,item_nchan,emptyList,item_names,item_vectors);
1900  MneProjOp::mne_proj_op_add_item_act(op,item,item_kind,item_desc,item_active);
1901  delete item;
1902  op->items[op->nitems-1]->active_file = item_active;
1903  }
1904 
1905 out :
1906  return op;
1907 
1908 bad : {
1909  if(op)
1910  delete op;
1911  return NULL;
1912  }
1913 }
1914 
1915 MneProjOp* mne_read_proj_op_3(const QString& name)
1916 
1917 {
1918  QFile file(name);
1919  FiffStream::SPtr stream(new FiffStream(&file));
1920 
1921  if(!stream->open())
1922  return NULL;
1923 
1924  MneProjOp* res = NULL;
1925 
1926  FiffDirNode::SPtr t_default;
1927  res = mne_read_proj_op_from_node_3(stream,t_default);
1928 
1929  stream->close();
1930 
1931  return res;
1932 }
1933 
1934 int mne_proj_op_chs_3(MneProjOp* op, const QStringList& list, int nlist)
1935 
1936 {
1937  if (op == NULL)
1938  return OK;
1939 
1940  MneProjOp::mne_free_proj_op_proj(op); /* These data are not valid any more */
1941 
1942  if (nlist == 0)
1943  return OK;
1944 
1945  op->names = list;
1946  op->nch = nlist;
1947 
1948  return OK;
1949 }
1950 
1951 static void clear_these(float *data, const QStringList& names, int nnames, const QString& start)
1952 
1953 {
1954  int k;
1955  for (k = 0; k < nnames; k++)
1956  if (names[k].contains(start))//strstr(names[k],start) == names[k])
1957  data[k] = 0.0;
1958 }
1959 
1960 #define USE_LIMIT 1e-5
1961 #define SMALL_VALUE 1e-4
1962 
1963 int mne_proj_op_make_proj_bad(MneProjOp* op, char **bad, int nbad)
1964 /*
1965  * Do the channel picking and SVD
1966  * Include a bad list at this phase
1967  * Input to the projection can include the bad channels
1968  * but they are not affected
1969  */
1970 {
1971  int k,p,q,r,nvec;
1972  float **vv_meg = NULL;
1973  float *sing_meg = NULL;
1974  float **vv_eeg = NULL;
1975  float *sing_eeg = NULL;
1976  float **mat_meg = NULL;
1977  float **mat_eeg = NULL;
1978  int nvec_meg;
1979  int nvec_eeg;
1980  mneNamedVectorRec vec;
1981  float size;
1982  int nzero;
1983 #ifdef DEBUG
1984  char name[20];
1985 #endif
1986 
1987  FREE_CMATRIX_3(op->proj_data);
1988  op->proj_data = NULL;
1989  op->nvec = 0;
1990 
1991  if (op->nch <= 0)
1992  return OK;
1993  if (op->nitems <= 0)
1994  return OK;
1995 
1996  nvec = MneProjOp::mne_proj_op_affect(op,op->names,op->nch);
1997  if (nvec == 0)
1998  return OK;
1999 
2000  mat_meg = ALLOC_CMATRIX_3(nvec, op->nch);
2001  mat_eeg = ALLOC_CMATRIX_3(nvec, op->nch);
2002 
2003 #ifdef DEBUG
2004  fprintf(stdout,"mne_proj_op_make_proj_bad\n");
2005 #endif
2006  for (k = 0, nvec_meg = nvec_eeg = 0; k < op->nitems; k++) {
2007  if (op->items[k]->active && MneProjItem::mne_proj_item_affect(op->items[k],op->names,op->nch)) {
2008  vec.nvec = op->items[k]->vecs->ncol;
2009  vec.names = op->items[k]->vecs->collist;
2010  if (op->items[k]->has_meg) {
2011  for (p = 0; p < op->items[k]->nvec; p++, nvec_meg++) {
2012  vec.data = op->items[k]->vecs->data[p];
2013  if (mne_pick_from_named_vector_3(&vec,op->names,op->nch,FALSE,mat_meg[nvec_meg]) == FAIL)
2014  goto bad;
2015 #ifdef DEBUG
2016  fprintf(stdout,"Original MEG:\n");
2017  mne_print_vector(stdout,op->items[k]->desc,mat_meg[nvec_meg],op->nch);
2018  fflush(stdout);
2019 #endif
2020  }
2021  }
2022  else if (op->items[k]->has_eeg) {
2023  for (p = 0; p < op->items[k]->nvec; p++, nvec_eeg++) {
2024  vec.data = op->items[k]->vecs->data[p];
2025  if (mne_pick_from_named_vector_3(&vec,op->names,op->nch,FALSE,mat_eeg[nvec_eeg]) == FAIL)
2026  goto bad;
2027 #ifdef DEBUG
2028  fprintf (stdout,"Original EEG:\n");
2029  mne_print_vector(stdout,op->items[k]->desc,mat_eeg[nvec_eeg],op->nch);
2030  fflush(stdout);
2031 #endif
2032  }
2033  }
2034  }
2035  }
2036  /*
2037  * Replace bad channel entries with zeroes
2038  */
2039  for (q = 0; q < nbad; q++)
2040  for (r = 0; r < op->nch; r++)
2041  if (QString::compare(op->names[r],bad[q]) == 0) {
2042  for (p = 0; p < nvec_meg; p++)
2043  mat_meg[p][r] = 0.0;
2044  for (p = 0; p < nvec_eeg; p++)
2045  mat_eeg[p][r] = 0.0;
2046  }
2047  /*
2048  * Scale the rows so that detection of linear dependence becomes easy
2049  */
2050  for (p = 0, nzero = 0; p < nvec_meg; p++) {
2051  size = sqrt(mne_dot_vectors_3(mat_meg[p],mat_meg[p],op->nch));
2052  if (size > 0) {
2053  for (k = 0; k < op->nch; k++)
2054  mat_meg[p][k] = mat_meg[p][k]/size;
2055  }
2056  else
2057  nzero++;
2058  }
2059  if (nzero == nvec_meg) {
2060  FREE_CMATRIX_3(mat_meg); mat_meg = NULL; nvec_meg = 0;
2061  }
2062  for (p = 0, nzero = 0; p < nvec_eeg; p++) {
2063  size = sqrt(mne_dot_vectors_3(mat_eeg[p],mat_eeg[p],op->nch));
2064  if (size > 0) {
2065  for (k = 0; k < op->nch; k++)
2066  mat_eeg[p][k] = mat_eeg[p][k]/size;
2067  }
2068  else
2069  nzero++;
2070  }
2071  if (nzero == nvec_eeg) {
2072  FREE_CMATRIX_3(mat_eeg); mat_eeg = NULL; nvec_eeg = 0;
2073  }
2074  if (nvec_meg + nvec_eeg == 0) {
2075  printf("No projection remains after excluding bad channels. Omitting projection.\n");
2076  return OK;
2077  }
2078  /*
2079  * Proceed to SVD
2080  */
2081 #ifdef DEBUG
2082  fprintf(stdout,"Before SVD:\n");
2083 #endif
2084  if (nvec_meg > 0) {
2085 #ifdef DEBUG
2086  fprintf(stdout,"---->>\n");
2087  for (p = 0; p < nvec_meg; p++) {
2088  sprintf(name,"MEG %02d",p+1);
2089  mne_print_vector(stdout,name,mat_meg[p],op->nch);
2090  }
2091  fprintf(stdout,"---->>\n");
2092 #endif
2093  sing_meg = MALLOC_3(nvec_meg+1,float);
2094  vv_meg = ALLOC_CMATRIX_3(nvec_meg,op->nch);
2095  if (mne_svd_3(mat_meg,nvec_meg,op->nch,sing_meg,NULL,vv_meg) != OK)
2096  goto bad;
2097  }
2098  if (nvec_eeg > 0) {
2099 #ifdef DEBUG
2100  fprintf(stdout,"---->>\n");
2101  for (p = 0; p < nvec_eeg; p++) {
2102  sprintf(name,"EEG %02d",p+1);
2103  mne_print_vector(stdout,name,mat_eeg[p],op->nch);
2104  }
2105  fprintf(stdout,"---->>\n");
2106 #endif
2107  sing_eeg = MALLOC_3(nvec_eeg+1,float);
2108  vv_eeg = ALLOC_CMATRIX_3(nvec_eeg,op->nch);
2109  if (mne_svd_3(mat_eeg,nvec_eeg,op->nch,sing_eeg,NULL,vv_eeg) != OK)
2110  goto bad;
2111  }
2112  /*
2113  * Check for linearly dependent vectors
2114  */
2115  for (p = 0, op->nvec = 0; p < nvec_meg; p++, op->nvec++)
2116  if (sing_meg[p]/sing_meg[0] < USE_LIMIT)
2117  break;
2118  for (p = 0; p < nvec_eeg; p++, op->nvec++)
2119  if (sing_eeg[p]/sing_eeg[0] < USE_LIMIT)
2120  break;
2121 #ifdef DEBUG
2122  printf("Number of linearly independent vectors = %d\n",op->nvec);
2123 #endif
2124  op->proj_data = ALLOC_CMATRIX_3(op->nvec,op->nch);
2125 #ifdef DEBUG
2126  fprintf(stdout,"Final projection data:\n");
2127 #endif
2128  for (p = 0, op->nvec = 0; p < nvec_meg; p++, op->nvec++) {
2129  if (sing_meg[p]/sing_meg[0] < USE_LIMIT)
2130  break;
2131  for (k = 0; k < op->nch; k++) {
2132  /*
2133  * Avoid crosstalk between MEG/EEG
2134  */
2135  if (std::fabs(vv_meg[p][k]) < SMALL_VALUE)
2136  op->proj_data[op->nvec][k] = 0.0;
2137  else
2138  op->proj_data[op->nvec][k] = vv_meg[p][k];
2139 
2140  /*
2141  * If the above did not work, this will (provided that EEG channels are called EEG*)
2142  */
2143  clear_these(op->proj_data[op->nvec],op->names,op->nch,"EEG");
2144  }
2145 #ifdef DEBUG
2146  sprintf(name,"MEG %02d",p+1);
2147  mne_print_vector(stdout,name,op->proj_data[op->nvec],op->nch);
2148 #endif
2149  }
2150  for (p = 0; p < nvec_eeg; p++, op->nvec++) {
2151  if (sing_eeg[p]/sing_eeg[0] < USE_LIMIT)
2152  break;
2153  for (k = 0; k < op->nch; k++) {
2154  /*
2155  * Avoid crosstalk between MEG/EEG
2156  */
2157  if (std::fabs(vv_eeg[p][k]) < SMALL_VALUE)
2158  op->proj_data[op->nvec][k] = 0.0;
2159  else
2160  op->proj_data[op->nvec][k] = vv_eeg[p][k];
2161  /*
2162  * If the above did not work, this will (provided that MEG channels are called MEG*)
2163  */
2164  clear_these(op->proj_data[op->nvec],op->names,op->nch,"MEG");
2165  }
2166 #ifdef DEBUG
2167  sprintf(name,"EEG %02d",p+1);
2168  mne_print_vector(stdout,name,op->proj_data[op->nvec],op->nch);
2169  fflush(stdout);
2170 #endif
2171  }
2172  FREE_3(sing_meg);
2173  FREE_CMATRIX_3(vv_meg);
2174  FREE_CMATRIX_3(mat_meg);
2175  FREE_3(sing_eeg);
2176  FREE_CMATRIX_3(vv_eeg);
2177  FREE_CMATRIX_3(mat_eeg);
2178  /*
2179  * Make sure that the stimulus channels are not modified
2180  */
2181  for (k = 0; k < op->nch; k++)
2182  if (op->names[k].contains("STI")){ //strstr(op->names[k],"STI") == op->names[k]) {
2183  for (p = 0; p < op->nvec; p++)
2184  op->proj_data[p][k] = 0.0;
2185  }
2186 
2187  return OK;
2188 
2189 bad : {
2190  FREE_3(sing_meg);
2191  FREE_CMATRIX_3(vv_meg);
2192  FREE_CMATRIX_3(mat_meg);
2193  FREE_3(sing_eeg);
2194  FREE_CMATRIX_3(vv_eeg);
2195  FREE_CMATRIX_3(mat_eeg);
2196  return FAIL;
2197  }
2198 }
2199 
2200 int mne_proj_op_make_proj(MneProjOp* op)
2201 /*
2202  * Do the channel picking and SVD
2203  */
2204 {
2205  return mne_proj_op_make_proj_bad(op,NULL,0);
2206 }
2207 
2208 //============================= mne_read_forward_solution.c =============================
2209 
2210 int mne_read_meg_comp_eeg_ch_info_3(const QString& name,
2211  QList<FiffChInfo>& megp, /* MEG channels */
2212  int *nmegp,
2213  QList<FiffChInfo>& meg_compp,
2214  int *nmeg_compp,
2215  QList<FiffChInfo>& eegp, /* EEG channels */
2216  int *neegp,
2217  FiffCoordTransOld** meg_head_t,
2218  fiffId *idp) /* The measurement ID */
2219 /*
2220  * Read the channel information and split it into three arrays,
2221  * one for MEG, one for MEG compensation channels, and one for EEG
2222  */
2223 {
2224  QFile file(name);
2225  FiffStream::SPtr stream(new FiffStream(&file));
2226 
2227  QList<FiffChInfo> chs;
2228  int nchan = 0;
2229  QList<FiffChInfo> meg;
2230  int nmeg = 0;
2231  QList<FiffChInfo> meg_comp;
2232  int nmeg_comp = 0;
2233  QList<FiffChInfo> eeg;
2234  int neeg = 0;
2235  fiffId id = NULL;
2236  QList<FiffDirNode::SPtr> nodes;
2237  FiffDirNode::SPtr info;
2238  FiffTag::SPtr t_pTag;
2239  FiffChInfo this_ch;
2240  FiffCoordTransOld* t = NULL;
2241  fiff_int_t kind, pos;
2242  int j,k,to_find;
2243 
2244  if(!stream->open())
2245  goto bad;
2246 
2247  nodes = stream->dirtree()->dir_tree_find(FIFFB_MNE_PARENT_MEAS_FILE);
2248 
2249  if (nodes.size() == 0) {
2250  nodes = stream->dirtree()->dir_tree_find(FIFFB_MEAS_INFO);
2251  if (nodes.size() == 0) {
2252  qCritical ("Could not find the channel information.");
2253  goto bad;
2254  }
2255  }
2256  info = nodes[0];
2257  to_find = 0;
2258  for (k = 0; k < info->nent(); k++) {
2259  kind = info->dir[k]->kind;
2260  pos = info->dir[k]->pos;
2261  switch (kind) {
2262  case FIFF_NCHAN :
2263  if (!stream->read_tag(t_pTag,pos))
2264  goto bad;
2265  nchan = *t_pTag->toInt();
2266 
2267  for (j = 0; j < nchan; j++)
2268  chs.append(FiffChInfo());
2269  chs[j].scanNo = -1;
2270  to_find = nchan;
2271  break;
2272 
2273  case FIFF_PARENT_BLOCK_ID :
2274  if(!stream->read_tag(t_pTag, pos))
2275  goto bad;
2276  *id = *(fiffId)t_pTag->data();
2277  break;
2278 
2279  case FIFF_COORD_TRANS :
2280  if(!stream->read_tag(t_pTag, pos))
2281  goto bad;
2282  t = FiffCoordTransOld::read_helper( t_pTag );
2283  if (t && (t->from != FIFFV_COORD_DEVICE || t->to != FIFFV_COORD_HEAD))
2284  {
2285  delete t;
2286  t = NULL;
2287  }
2288  break;
2289 
2290  case FIFF_CH_INFO : /* Information about one channel */
2291  if(!stream->read_tag(t_pTag, pos))
2292  goto bad;
2293  this_ch = t_pTag->toChInfo();
2294  if (this_ch.scanNo <= 0 || this_ch.scanNo > nchan) {
2295  printf ("FIFF_CH_INFO : scan # out of range %d (%d)!",this_ch.scanNo,nchan);
2296  goto bad;
2297  }
2298  else
2299  chs[this_ch.scanNo-1] = this_ch;
2300  to_find--;
2301  break;
2302  }
2303  }
2304  if (to_find != 0) {
2305  qCritical("Some of the channel information was missing.");
2306  goto bad;
2307  }
2308  if (t == NULL && meg_head_t != NULL) {
2309  /*
2310  * Try again in a more general fashion
2311  */
2312  if ((t = FiffCoordTransOld::mne_read_meas_transform(name)) == NULL) {
2313  qCritical("MEG -> head coordinate transformation not found.");
2314  goto bad;
2315  }
2316  }
2317  /*
2318  * Sort out the channels
2319  */
2320  for (k = 0; k < nchan; k++) {
2321  if (chs[k].kind == FIFFV_MEG_CH) {
2322  meg.append(chs[k]);
2323  nmeg++;
2324  } else if (chs[k].kind == FIFFV_REF_MEG_CH) {
2325  meg_comp.append(chs[k]);
2326  nmeg_comp++;
2327  } else if (chs[k].kind == FIFFV_EEG_CH) {
2328  eeg.append(chs[k]);
2329  neeg++;
2330  }
2331  }
2332 
2333  stream->close();
2334 
2335  megp = meg;
2336  if(nmegp) {
2337  *nmegp = nmeg;
2338  }
2339 
2340  meg_compp = meg_comp;
2341  *nmeg_compp = nmeg_comp;
2342 
2343  eegp = eeg;
2344  *neegp = neeg;
2345 
2346  if (idp == NULL) {
2347  FREE_3(id);
2348  }
2349  else
2350  *idp = id;
2351  if (meg_head_t == NULL) {
2352  delete t;
2353  }
2354  else
2355  *meg_head_t = t;
2356 
2357  return FIFF_OK;
2358 
2359 bad : {
2360 
2361  stream->close();
2362  FREE_3(id);
2363 
2364  delete t;
2365  return FIFF_FAIL;
2366  }
2367 }
2368 
2369 //============================= data.c =============================
2370 
2371 #if defined(_WIN32) || defined(_WIN64)
2372 #define snprintf _snprintf
2373 #define vsnprintf _vsnprintf
2374 #define strcasecmp _stricmp
2375 #define strncasecmp _strnicmp
2376 #endif
2377 
2378 int is_selected_in_data(mshMegEegData d, const QString& ch_name)
2379 /*
2380  * Is this channel selected in data
2381  */
2382 {
2383  int issel = FALSE;
2384  int k;
2385 
2386  for (k = 0; k < d->meas->nchan; k++)
2387  if (QString::compare(ch_name,d->meas->chs[k].ch_name,Qt::CaseInsensitive) == 0) {
2388  issel = d->sels[k];
2389  break;
2390  }
2391  return issel;
2392 }
2393 
2394 //============================= mne_process_bads.c =============================
2395 
2396 static int whitespace_3(char *text)
2397 
2398 {
2399  if (text == NULL || strlen(text) == 0)
2400  return TRUE;
2401  if (strspn(text," \t\n\r") == strlen(text))
2402  return TRUE;
2403  return FALSE;
2404 }
2405 
2406 static char *next_line_3(char *line, int n, FILE *in)
2407 {
2408  char *res;
2409 
2410  for (res = fgets(line,n,in); res != NULL; res = fgets(line,n,in))
2411  if (!whitespace_3(res))
2412  if (res[0] != '#')
2413  break;
2414  return res;
2415 }
2416 
2417 #define MAXLINE 500
2418 
2419 int mne_read_bad_channels_3(const QString& name, QStringList& listp, int& nlistp)
2420 /*
2421  * Read bad channel names
2422  */
2423 {
2424  FILE *in = NULL;
2425  QStringList list;
2426  char line[MAXLINE+1];
2427  char *next;
2428 
2429  if (name.isEmpty())
2430  return OK;
2431 
2432  if ((in = fopen(name.toUtf8().data(),"r")) == NULL) {
2433  qCritical() << name;
2434  goto bad;
2435  }
2436  while ((next = next_line_3(line,MAXLINE,in)) != NULL) {
2437  if (strlen(next) > 0) {
2438  if (next[strlen(next)-1] == '\n')
2439  next[strlen(next)-1] = '\0';
2440  list.append(next);
2441  }
2442  }
2443  if (ferror(in))
2444  goto bad;
2445 
2446  listp = list;
2447  nlistp = list.size();
2448 
2449  return OK;
2450 
2451 bad : {
2452  list.clear();
2453  if (in != NULL)
2454  fclose(in);
2455  return FAIL;
2456  }
2457 }
2458 
2459 int mne_read_bad_channel_list_from_node_3(FiffStream::SPtr& stream,
2460  const FiffDirNode::SPtr& pNode,
2461  QStringList& listp,
2462  int& nlistp)
2463 {
2464  FiffDirNode::SPtr node,bad;
2465  QList<FiffDirNode::SPtr> temp;
2466  QStringList list;
2467  int nlist = 0;
2468  FiffTag::SPtr t_pTag;
2469  QString names;
2470 
2471  if (pNode->isEmpty())
2472  node = stream->dirtree();
2473  else
2474  node = pNode;
2475 
2476  temp = node->dir_tree_find(FIFFB_MNE_BAD_CHANNELS);
2477  if (temp.size() > 0) {
2478  bad = temp[0];
2479 
2480  bad->find_tag(stream, FIFF_MNE_CH_NAME_LIST, t_pTag);
2481  if (t_pTag) {
2482  names = t_pTag->toString();
2483  mne_string_to_name_list_3(names,list,nlist);
2484  }
2485  }
2486  listp = list;
2487  nlistp = list.size();
2488  return OK;
2489 }
2490 
2491 int mne_read_bad_channel_list_3(const QString& name, QStringList& listp, int& nlistp)
2492 
2493 {
2494  QFile file(name);
2495  FiffStream::SPtr stream(new FiffStream(&file));
2496 
2497  int res;
2498 
2499  if(!stream->open())
2500  return FAIL;
2501 
2502  res = mne_read_bad_channel_list_from_node_3(stream,stream->dirtree(),listp,nlistp);
2503 
2504  stream->close();
2505 
2506  return res;
2507 }
2508 
2509 MneCovMatrix* mne_read_cov(const QString& name,int kind)
2510 /*
2511  * Read a covariance matrix from a fiff
2512  */
2513 {
2514  QFile file(name);
2515  FiffStream::SPtr stream(new FiffStream(&file));
2516 
2517  FiffTag::SPtr t_pTag;
2518  QList<FiffDirNode::SPtr> nodes;
2519  FiffDirNode::SPtr covnode;
2520 
2521  QStringList names; /* Optional channel name list */
2522  int nnames = 0;
2523  double *cov = NULL;
2524  double *cov_diag = NULL;
2525  FiffSparseMatrix* cov_sparse = NULL;
2526  double *lambda = NULL;
2527  float **eigen = NULL;
2528  MatrixXf tmp_eigen;
2529  QStringList bads;
2530  int nbad = 0;
2531  int ncov = 0;
2532  int nfree = 1;
2533  MneCovMatrix* res = NULL;
2534 
2535  int k,p,nn;
2536  float *f;
2537  double *d;
2538  MneProjOp* op = NULL;
2539  MneSssData* sss = NULL;
2540 
2541  if(!stream->open())
2542  goto out;
2543 
2544  nodes = stream->dirtree()->dir_tree_find(FIFFB_MNE_COV);
2545 
2546  if (nodes.size() == 0) {
2547  printf("No covariance matrix available in %s",name.toUtf8().data());
2548  goto out;
2549  }
2550  /*
2551  * Locate the desired matrix
2552  */
2553  for (k = 0 ; k < nodes.size(); ++k) {
2554  if (!nodes[k]->find_tag(stream, FIFF_MNE_COV_KIND, t_pTag))
2555  continue;
2556 
2557  if (*t_pTag->toInt() == kind) {
2558  covnode = nodes[k];
2559  break;
2560  }
2561  }
2562  if (covnode->isEmpty()) {
2563  printf("Desired covariance matrix not found from %s",name.toUtf8().data());
2564  goto out;
2565  }
2566  /*
2567  * Read the data
2568  */
2569  if (!nodes[k]->find_tag(stream, FIFF_MNE_COV_DIM, t_pTag))
2570  goto out;
2571  ncov = *t_pTag->toInt();
2572 
2573  if (nodes[k]->find_tag(stream, FIFF_MNE_COV_NFREE, t_pTag)) {
2574  nfree = *t_pTag->toInt();
2575  }
2576  if (covnode->find_tag(stream, FIFF_MNE_ROW_NAMES, t_pTag)) {
2577  mne_string_to_name_list_3(t_pTag->toString(),names,nnames);
2578  if (nnames != ncov) {
2579  qCritical("Incorrect number of channel names for a covariance matrix");
2580  goto out;
2581  }
2582  }
2583  if (!nodes[k]->find_tag(stream, FIFF_MNE_COV, t_pTag)) {
2584  if (!nodes[k]->find_tag(stream, FIFF_MNE_COV_DIAG, t_pTag))
2585  goto out;
2586  else {
2587  if (t_pTag->getType() == FIFFT_DOUBLE) {
2588  cov_diag = MALLOC_3(ncov,double);
2589  qDebug() << "ToDo: Check whether cov_diag contains the right stuff!!! - use VectorXd instead";
2590  d = t_pTag->toDouble();
2591  for (p = 0; p < ncov; p++)
2592  cov_diag[p] = d[p];
2593  if (check_cov_data(cov_diag,ncov) != OK)
2594  goto out;
2595  }
2596  else if (t_pTag->getType() == FIFFT_FLOAT) {
2597  cov_diag = MALLOC_3(ncov,double);
2598  qDebug() << "ToDo: Check whether f contains the right stuff!!! - use VectorXf instead";
2599  f = t_pTag->toFloat();
2600  for (p = 0; p < ncov; p++)
2601  cov_diag[p] = f[p];
2602  }
2603  else {
2604  printf("Illegal data type for covariance matrix");
2605  goto out;
2606  }
2607  }
2608  }
2609  else {
2610  nn = ncov*(ncov+1)/2;
2611  if (t_pTag->getType() == FIFFT_DOUBLE) {
2612  qDebug() << "ToDo: Check whether cov contains the right stuff!!! - use VectorXd instead";
2613 // qDebug() << t_pTag->size() / sizeof(double);
2614  cov = MALLOC_3(nn,double);
2615  d = t_pTag->toDouble();
2616  for (p = 0; p < nn; p++)
2617  cov[p] = d[p];
2618  if (check_cov_data(cov,nn) != OK)
2619  goto out;
2620  }
2621  else if (t_pTag->getType() == FIFFT_FLOAT) {
2622  cov = MALLOC_3(nn,double);
2623  f = t_pTag->toFloat();
2624  for (p = 0; p < nn; p++)
2625  cov[p] = f[p];
2626  }
2627  else if ((cov_sparse = FiffSparseMatrix::fiff_get_float_sparse_matrix(t_pTag)) == NULL) {
2628  goto out;
2629  }
2630 
2631  if (nodes[k]->find_tag(stream, FIFF_MNE_COV_EIGENVALUES, t_pTag)) {
2632  lambda = (double *)t_pTag->toDouble();
2633  if (nodes[k]->find_tag(stream, FIFF_MNE_COV_EIGENVECTORS, t_pTag))
2634  goto out;
2635 
2636  tmp_eigen = t_pTag->toFloatMatrix().transpose();
2637  eigen = ALLOC_CMATRIX_3(tmp_eigen.rows(),tmp_eigen.cols());
2638  fromFloatEigenMatrix_3(tmp_eigen, eigen);
2639  }
2640  /*
2641  * Read the optional projection operator
2642  */
2643  if ((op = mne_read_proj_op_from_node_3(stream,nodes[k])) == NULL)
2644  goto out;
2645  /*
2646  * Read the optional SSS data
2647  */
2648  if ((sss = MneSssData::read_sss_data_from_node(stream,nodes[k])) == NULL)
2649  goto out;
2650  /*
2651  * Read the optional bad channel list
2652  */
2653  if (mne_read_bad_channel_list_from_node_3(stream,nodes[k],bads,nbad) == FAIL)
2654  goto out;
2655  }
2656  if (cov_sparse)
2657  res = MneCovMatrix::mne_new_cov_sparse(kind,ncov,names,cov_sparse);
2658  else if (cov)
2659  res = MneCovMatrix::mne_new_cov_dense(kind,ncov,names,cov);
2660  else if (cov_diag)
2661  res = MneCovMatrix::mne_new_cov_diag(kind,ncov,names,cov_diag);
2662  else {
2663  qCritical("mne_read_cov : covariance matrix data is not defined. How come?");
2664  goto out;
2665  }
2666  cov = NULL;
2667  cov_diag = NULL;
2668  cov_sparse = NULL;
2669  res->eigen = eigen;
2670  res->lambda = lambda;
2671  res->nfree = nfree;
2672  res->bads = bads;
2673  res->nbad = nbad;
2674  /*
2675  * Count the non-zero eigenvalues
2676  */
2677  if (res->lambda) {
2678  res->nzero = 0;
2679  for (k = 0; k < res->ncov; k++, res->nzero++)
2680  if (res->lambda[k] > 0)
2681  break;
2682  }
2683 
2684  if (op && op->nitems > 0) {
2685  res->proj = op;
2686  op = NULL;
2687  }
2688  if (sss && sss->ncomp > 0 && sss->job != FIFFV_SSS_JOB_NOTHING) {
2689  res->sss = sss;
2690  sss = NULL;
2691  }
2692 
2693 out : {
2694  stream->close();
2695  if(op)
2696  delete op;
2697  if(sss)
2698  delete sss;
2699 
2700  if (!res) {
2701  names.clear();
2702  bads.clear();
2703  FREE_3(cov);
2704  FREE_3(cov_diag);
2705  if(cov_sparse)
2706  delete cov_sparse;
2707  }
2708  return res;
2709  }
2710 }
2711 
2712 //============================= mne_coord_transforms.c =============================
2713 
2714 namespace MNELIB
2715 {
2716 
2717 typedef struct {
2718  int frame;
2719  const char *name;
2720 } frameNameRec_3;
2721 
2722 }
2723 
2724 const char *mne_coord_frame_name_3(int frame)
2725 
2726 {
2727  static frameNameRec_3 frames[] = {
2728  {FIFFV_COORD_UNKNOWN,"unknown"},
2729  {FIFFV_COORD_DEVICE,"MEG device"},
2730  {FIFFV_COORD_ISOTRAK,"isotrak"},
2731  {FIFFV_COORD_HPI,"hpi"},
2732  {FIFFV_COORD_HEAD,"head"},
2733  {FIFFV_COORD_MRI,"MRI (surface RAS)"},
2734  {FIFFV_MNE_COORD_MRI_VOXEL, "MRI voxel"},
2735  {FIFFV_COORD_MRI_SLICE,"MRI slice"},
2736  {FIFFV_COORD_MRI_DISPLAY,"MRI display"},
2737  {FIFFV_MNE_COORD_CTF_DEVICE,"CTF MEG device"},
2738  {FIFFV_MNE_COORD_CTF_HEAD,"CTF/4D/KIT head"},
2739  {FIFFV_MNE_COORD_RAS,"RAS (non-zero origin)"},
2740  {FIFFV_MNE_COORD_MNI_TAL,"MNI Talairach"},
2741  {FIFFV_MNE_COORD_FS_TAL_GTZ,"Talairach (MNI z > 0)"},
2742  {FIFFV_MNE_COORD_FS_TAL_LTZ,"Talairach (MNI z < 0)"},
2743  {-1,"unknown"}
2744  };
2745  int k;
2746  for (k = 0; frames[k].frame != -1; k++) {
2747  if (frame == frames[k].frame)
2748  return frames[k].name;
2749  }
2750  return frames[k].name;
2751 }
2752 
2753 //============================= mne_read_process_forward_solution.c =============================
2754 
2755 void mne_merge_channels(const QList<FiffChInfo>& chs1,
2756  int nch1,
2757  const QList<FiffChInfo>& chs2,
2758  int nch2,
2759  QList<FiffChInfo>& resp,
2760  int *nresp)
2761 
2762 {
2763  resp.clear();
2764  resp.reserve(nch1+nch2);
2765 
2766  int k;
2767  for (k = 0; k < nch1; k++) {
2768  resp.append(chs1.at(k));
2769  }
2770  for (k = 0; k < nch2; k++) {
2771  resp.append(chs2.at(k));
2772  }
2773 
2774  *nresp = nch1+nch2;
2775  return;
2776 }
2777 
2778 //============================= read_ch_info.c =============================
2779 
2780 static FiffDirNode::SPtr find_meas_3 (const FiffDirNode::SPtr& node)
2781 /*
2782  * Find corresponding meas node
2783  */
2784 {
2785  FiffDirNode::SPtr empty_node;
2786  FiffDirNode::SPtr tmp_node = node;
2787 
2788  while (tmp_node->type != FIFFB_MEAS) {
2789  if (tmp_node->parent == NULL)
2790  return empty_node;//(NULL);
2791  tmp_node = tmp_node->parent;
2792  }
2793  return (tmp_node);
2794 }
2795 
2796 static FiffDirNode::SPtr find_meas_info_3 (const FiffDirNode::SPtr& node)
2797 /*
2798  * Find corresponding meas info node
2799  */
2800 {
2801  int k;
2802  FiffDirNode::SPtr empty_node;
2803  FiffDirNode::SPtr tmp_node = node;
2804 
2805  while (tmp_node->type != FIFFB_MEAS) {
2806  if (tmp_node->parent == NULL)
2807  return empty_node;
2808  tmp_node = tmp_node->parent;
2809  }
2810  for (k = 0; k < tmp_node->nchild(); k++)
2811  if (tmp_node->children[k]->type == FIFFB_MEAS_INFO)
2812  return (tmp_node->children[k]);
2813  return empty_node;
2814 }
2815 
2816 static int get_all_chs (//fiffFile file, /* The file we are reading */
2817  FiffStream::SPtr& stream,
2818  const FiffDirNode::SPtr& p_node, /* The directory node containing our data */
2819  fiffId *id, /* The block id from the nearest FIFFB_MEAS parent */
2820  QList<FiffChInfo>& chp, /* Channel descriptions */
2821  int *nchanp) /* Number of channels */
2822 /*
2823  * Find channel information from
2824  * nearest FIFFB_MEAS_INFO parent of
2825  * node.
2826  */
2827 {
2828  QList<FiffChInfo> ch;
2829  FiffChInfo this_ch;
2830  int j,k;
2831  int nchan = 0;
2832  int to_find = 0;
2833  FiffDirNode::SPtr meas;
2834  FiffDirNode::SPtr meas_info;
2835  fiff_int_t kind, pos;
2836  FiffTag::SPtr t_pTag;
2837 
2838  *id = NULL;
2839  /*
2840  * Find desired parents
2841  */
2842  if (!(meas = find_meas_3(p_node))) {
2843  qCritical("Meas. block not found!");
2844  goto bad;
2845  }
2846 
2847  if (!(meas_info = find_meas_info_3(p_node))) {
2848  qCritical ("Meas. info not found!");
2849  goto bad;
2850  }
2851  /*
2852  * Is there a block id is in the FIFFB_MEAS node?
2853  */
2854  if (!meas->id.isEmpty()) {
2855  *id = MALLOC_3(1,fiffIdRec);
2856  (*id)->version = meas->id.version;
2857  (*id)->machid[0] = meas->id.machid[0];
2858  (*id)->machid[1] = meas->id.machid[1];
2859  (*id)->time = meas->id.time;
2860  }
2861  /*
2862  * Others from FIFFB_MEAS_INFO
2863  */
2864  for (k = 0; k < meas_info->nent(); k++) {
2865  kind = meas_info->dir[k]->kind;
2866  pos = meas_info->dir[k]->pos;
2867  switch (kind) {
2868  case FIFF_NCHAN :
2869  if (!stream->read_tag(t_pTag,pos))
2870  goto bad;
2871  nchan = *t_pTag->toInt();
2872  *nchanp = nchan;
2873 
2874  for (j = 0; j < nchan; j++) {
2875  ch.append(FiffChInfo());
2876  ch[j].scanNo = -1;
2877  }
2878 
2879  to_find += nchan - 1;
2880  break;
2881 
2882  case FIFF_CH_INFO : /* Information about one channel */
2883  if (!stream->read_tag(t_pTag,pos))
2884  goto bad;
2885  this_ch = t_pTag->toChInfo();
2886  if (this_ch.scanNo <= 0 || this_ch.scanNo > nchan) {
2887  qCritical ("FIFF_CH_INFO : scan # out of range!");
2888  goto bad;
2889  }
2890  else
2891  ch[this_ch.scanNo-1] = this_ch;
2892  to_find--;
2893  break;
2894  }
2895  }
2896 
2897  chp = ch;
2898  return FIFF_OK;
2899 
2900 bad : {
2901  return FIFF_FAIL;
2902  }
2903 }
2904 
2905 static int read_ch_info(const QString& name,
2906  QList<FiffChInfo>& chsp,
2907  int *nchanp,
2908  fiffId *idp)
2909 /*
2910  * Read channel information from a measurement file
2911  */
2912 {
2913  QFile file(name);
2914  FiffStream::SPtr stream(new FiffStream(&file));
2915 
2916  QList<FiffChInfo> chs;
2917  int nchan = 0;
2918  fiffId id = NULL;
2919 
2920  QList<FiffDirNode::SPtr> meas;
2921  FiffDirNode::SPtr node;
2922 
2923  if(!stream->open())
2924  goto bad;
2925 
2926  meas = stream->dirtree()->dir_tree_find(FIFFB_MEAS);
2927  if (meas.size() == 0) {
2928  qCritical ("%s : no MEG data available here",name.toUtf8().data());
2929  goto bad;
2930  }
2931  node = meas[0];
2932  if (get_all_chs(stream,
2933  node,
2934  &id,
2935  chs,
2936  &nchan) == FIFF_FAIL)
2937  goto bad;
2938  chsp = chs;
2939  *nchanp = nchan;
2940  *idp = id;
2941  stream->close();
2942  return FIFF_OK;
2943 
2944 bad : {
2945  FREE_3(id);
2946  stream->close();
2947  return FIFF_FAIL;
2948  }
2949 }
2950 
2951 #define TOO_CLOSE 1e-4
2952 
2953 static int at_origin (const Eigen::Vector3f& rr)
2954 {
2955  return (rr.norm() < TOO_CLOSE);
2956 }
2957 
2958 static int is_valid_eeg_ch(const FiffChInfo& ch)
2959 /*
2960  * Is the electrode position information present?
2961  */
2962 {
2963  if (ch.kind == FIFFV_EEG_CH) {
2964  if (at_origin(ch.chpos.r0) ||
2966  return FALSE;
2967  else
2968  return TRUE;
2969  }
2970  return FALSE;
2971 }
2972 
2973 static int accept_ch(const FiffChInfo& ch,
2974  const QStringList& bads,
2975  int nbad)
2976 
2977 {
2978  int k;
2979  for (k = 0; k < nbad; k++)
2980  if (QString::compare(ch.ch_name,bads[k]) == 0)
2981  return FALSE;
2982  return TRUE;
2983 }
2984 
2985 int read_meg_eeg_ch_info(const QString& name, /* Input file */
2986  int do_meg, /* Use MEG */
2987  int do_eeg, /* Use EEG */
2988  const QStringList& bads, /* List of bad channels */
2989  int nbad,
2990  QList<FiffChInfo>& chsp, /* MEG + EEG channels */
2991  int *nmegp, /* Count of each */
2992  int *neegp)
2993 /*
2994  * Read the channel information and split it into two arrays,
2995  * one for MEG and one for EEG
2996  */
2997 {
2998  QList<FiffChInfo> chs;
2999  int nchan = 0;
3000  QList<FiffChInfo> meg;
3001  int nmeg = 0;
3002  QList<FiffChInfo> eeg;
3003  int neeg = 0;
3004  fiffId id = NULL;
3005  int nch;
3006 
3007  int k;
3008 
3009  if (read_ch_info(name,
3010  chs,
3011  &nchan,
3012  &id) != FIFF_OK)
3013  goto bad;
3014  /*
3015  * Sort out the channels
3016  */
3017  for (k = 0; k < nchan; k++) {
3018  if (accept_ch(chs[k],bads,nbad)) {
3019  if (do_meg && chs[k].kind == FIFFV_MEG_CH) {
3020  meg.append(chs[k]);
3021  nmeg++;
3022  } else if (do_eeg && chs[k].kind == FIFFV_EEG_CH && is_valid_eeg_ch(chs[k])) {
3023  eeg.append(chs[k]);
3024  neeg++;
3025  }
3026  }
3027  }
3028 
3029  mne_merge_channels(meg,
3030  nmeg,
3031  eeg,
3032  neeg,
3033  chsp,
3034  &nch);
3035 
3036  if(nmegp) {
3037  *nmegp = nmeg;
3038  }
3039  if(neegp) {
3040  *neegp = neeg;
3041  }
3042  FREE_3(id);
3043  return FIFF_OK;
3044 
3045 bad : {
3046  FREE_3(id);
3047  return FIFF_FAIL;
3048  }
3049 }
3050 
3051 void mne_revert_to_diag_cov(MneCovMatrix* c)
3052 /*
3053  * Pick the diagonal elements of the full covariance matrix
3054  */
3055 {
3056  int k,p;
3057  if (c->cov == NULL)
3058  return;
3059 #define REALLY_REVERT
3060 #ifdef REALLY_REVERT
3061  c->cov_diag = REALLOC_3(c->cov_diag,c->ncov,double);
3062 
3063  for (k = p = 0; k < c->ncov; k++) {
3064  c->cov_diag[k] = c->cov[p];
3065  p = p + k + 2;
3066  }
3067  FREE_3(c->cov); c->cov = NULL;
3068 #else
3069  for (j = 0, p = 0; j < c->ncov; j++)
3070  for (k = 0; k <= j; k++, p++)
3071  if (j != k)
3072  c->cov[p] = 0.0;
3073  else
3074 #endif
3075  FREE_3(c->lambda); c->lambda = NULL;
3076  FREE_CMATRIX_3(c->eigen); c->eigen = NULL;
3077  return;
3078 }
3079 
3080 MneCovMatrix* mne_pick_chs_cov_omit(MneCovMatrix* c,
3081  const QStringList& new_names,
3082  int ncov,
3083  int omit_meg_eeg,
3084  const QList<FiffChInfo>& chs)
3085 /*
3086  * Pick designated channels from a covariance matrix, optionally omit MEG/EEG correlations
3087  */
3088 {
3089  int j,k;
3090  int *pick = NULL;
3091  double *cov = NULL;
3092  double *cov_diag = NULL;
3093  QStringList names;
3094  int *is_meg = NULL;
3095  int from,to;
3096  MneCovMatrix* res;
3097 
3098  if (ncov == 0) {
3099  qCritical("No channels specified for picking in mne_pick_chs_cov_omit");
3100  return NULL;
3101  }
3102  if (c->names.isEmpty()) {
3103  qCritical("No names in covariance matrix. Cannot do picking.");
3104  return NULL;
3105  }
3106  pick = MALLOC_3(ncov,int);
3107  for (j = 0; j < ncov; j++)
3108  pick[j] = -1;
3109  for (j = 0; j < ncov; j++)
3110  for (k = 0; k < c->ncov; k++)
3111  if (QString::compare(c->names[k],new_names[j]) == 0) {
3112  pick[j] = k;
3113  break;
3114  }
3115  for (j = 0; j < ncov; j++) {
3116  if (pick[j] < 0) {
3117  printf("All desired channels not found in the covariance matrix (at least missing %s).", new_names[j].toUtf8().constData());
3118  FREE_3(pick);
3119  return NULL;
3120  }
3121  }
3122  if (omit_meg_eeg) {
3123  is_meg = MALLOC_3(ncov,int);
3124  if (!chs.isEmpty()) {
3125  for (j = 0; j < ncov; j++)
3126  if (chs[j].kind == FIFFV_MEG_CH)
3127  is_meg[j] = TRUE;
3128  else
3129  is_meg[j] = FALSE;
3130  }
3131  else {
3132  for (j = 0; j < ncov; j++)
3133  if (new_names[j].startsWith("MEG"))
3134  is_meg[j] = TRUE;
3135  else
3136  is_meg[j] = FALSE;
3137  }
3138  }
3139  if (c->cov_diag) {
3140  cov_diag = MALLOC_3(ncov,double);
3141  for (j = 0; j < ncov; j++) {
3142  cov_diag[j] = c->cov_diag[pick[j]];
3143  names.append(c->names[pick[j]]);
3144  }
3145  }
3146  else {
3147  cov = MALLOC_3(ncov*(ncov+1)/2,double);
3148  for (j = 0; j < ncov; j++) {
3149  names.append(c->names[pick[j]]);
3150  for (k = 0; k <= j; k++) {
3151  from = mne_lt_packed_index_3(pick[j],pick[k]);
3152  to = mne_lt_packed_index_3(j,k);
3153  if (to < 0 || to > ncov*(ncov+1)/2-1) {
3154  printf("Wrong destination index in mne_pick_chs_cov : %d %d %d\n",j,k,to);
3155  exit(1);
3156  }
3157  if (from < 0 || from > c->ncov*(c->ncov+1)/2-1) {
3158  printf("Wrong source index in mne_pick_chs_cov : %d %d %d\n",pick[j],pick[k],from);
3159  exit(1);
3160  }
3161  cov[to] = c->cov[from];
3162  if (omit_meg_eeg)
3163  if (is_meg[j] != is_meg[k])
3164  cov[to] = 0.0;
3165  }
3166  }
3167  }
3168 
3169  res = MneCovMatrix::mne_new_cov(c->kind,ncov,names,cov,cov_diag);
3170 
3171  res->bads = c->bads;
3172  res->nbad = c->nbad;
3173  res->proj = MneProjOp::mne_dup_proj_op(c->proj);
3174  res->sss = c->sss ? new MneSssData(*(c->sss)) : NULL;
3175 
3176  if (c->ch_class) {
3177  res->ch_class = MALLOC_3(res->ncov,int);
3178  for (k = 0; k < res->ncov; k++)
3179  res->ch_class[k] = c->ch_class[pick[k]];
3180  }
3181  FREE_3(pick);
3182  FREE_3(is_meg);
3183  return res;
3184 }
3185 
3186 int mne_proj_op_proj_dvector(MneProjOp* op, double *vec, int nch, int do_complement)
3187 /*
3188  * Apply projection operator to a vector (doubles)
3189  * Assume that all dimension checking etc. has been done before
3190  */
3191 {
3192  float *pvec;
3193  double w;
3194  int k,p;
3195 
3196  if (op->nvec <= 0)
3197  return OK;
3198 
3199  if (op->nch != nch) {
3200  qCritical("Data vector size does not match projection operator");
3201  return FAIL;
3202  }
3203 
3204  Eigen::VectorXd res = Eigen::VectorXd::Zero(op->nch);
3205 
3206  for (p = 0; p < op->nvec; p++) {
3207  pvec = op->proj_data[p];
3208  w = 0.0;
3209  for (k = 0; k < op->nch; k++)
3210  w += vec[k]*pvec[k];
3211  for (k = 0; k < op->nch; k++)
3212  res[k] = res[k] + w*pvec[k];
3213  }
3214 
3215  if (do_complement) {
3216  for (k = 0; k < op->nch; k++)
3217  vec[k] = vec[k] - res[k];
3218  }
3219  else {
3220  for (k = 0; k < op->nch; k++)
3221  vec[k] = res[k];
3222  }
3223 
3224  return OK;
3225 }
3226 
3227 int mne_name_list_match(const QStringList& list1, int nlist1,
3228  const QStringList& list2, int nlist2)
3229 /*
3230  * Check whether two name lists are identical
3231  */
3232 {
3233  int k;
3234  if (list1.isEmpty() && list2.isEmpty())
3235  return 0;
3236  if (list1.isEmpty() || list2.isEmpty())
3237  return 1;
3238  if (nlist1 != nlist2)
3239  return 1;
3240  for (k = 0; k < nlist1; k++)
3241  if (QString::compare(list1[k],list2[k]) != 0)
3242  return 1;
3243  return 0;
3244 }
3245 
3246 void mne_transpose_dsquare(double **mat, int n)
3247 /*
3248  * In-place transpose of a square matrix
3249  */
3250 {
3251  int j,k;
3252  double val;
3253 
3254  for (j = 1; j < n; j++)
3255  for (k = 0; k < j; k++) {
3256  val = mat[j][k];
3257  mat[j][k] = mat[k][j];
3258  mat[k][j] = val;
3259  }
3260  return;
3261 }
3262 
3263 int mne_proj_op_apply_cov(MneProjOp* op, MneCovMatrix*& c)
3264 /*
3265  * Apply the projection operator to a covariance matrix
3266  */
3267 {
3268  double **dcov = NULL;
3269  int j,k,p;
3270  int do_complement = TRUE;
3271 
3272  if (op == NULL || op->nitems == 0)
3273  return OK;
3274 
3275  if (mne_name_list_match(op->names,op->nch,c->names,c->ncov) != OK) {
3276  qCritical("Incompatible data in mne_proj_op_apply_cov");
3277  return FAIL;
3278  }
3279 
3280  dcov = ALLOC_DCMATRIX_3(c->ncov,c->ncov);
3281  /*
3282  * Return the appropriate result
3283  */
3284  if (c->cov_diag) { /* Pick the diagonals */
3285  for (j = 0, p = 0; j < c->ncov; j++)
3286  for (k = 0; k < c->ncov; k++)
3287  dcov[j][k] = (j == k) ? c->cov_diag[j] : 0;
3288  }
3289  else { /* Return full matrix */
3290  for (j = 0, p = 0; j < c->ncov; j++)
3291  for (k = 0; k <= j; k++)
3292  dcov[j][k] = c->cov[p++];
3293  for (j = 0; j < c->ncov; j++)
3294  for (k = j+1; k < c->ncov; k++)
3295  dcov[j][k] = dcov[k][j];
3296  }
3297 
3298  /*
3299  * Project from front and behind
3300  */
3301  for (k = 0; k < c->ncov; k++) {
3302  if (mne_proj_op_proj_dvector(op,dcov[k],c->ncov,do_complement) != OK)
3303  return FAIL;
3304  }
3305 
3306  mne_transpose_dsquare(dcov,c->ncov);
3307 
3308  for (k = 0; k < c->ncov; k++)
3309  if (mne_proj_op_proj_dvector(op,dcov[k],c->ncov,do_complement) != OK)
3310  return FAIL;
3311 
3312  /*
3313  * Return the result
3314  */
3315  if (c->cov_diag) { /* Pick the diagonal elements */
3316  for (j = 0; j < c->ncov; j++) {
3317  c->cov_diag[j] = dcov[j][j];
3318  }
3319  FREE_3(c->cov); c->cov = NULL;
3320  }
3321  else { /* Put everything back */
3322  for (j = 0, p = 0; j < c->ncov; j++)
3323  for (k = 0; k <= j; k++)
3324  c->cov[p++] = dcov[j][k];
3325  }
3326 
3327  FREE_DCMATRIX_3(dcov);
3328 
3329  c->nproj = MneProjOp::mne_proj_op_affect(op,c->names,c->ncov);
3330  return OK;
3331 }
3332 
3333 //=============================================================================================================
3334 // DEFINE MEMBER METHODS
3335 //=============================================================================================================
3336 
3337 DipoleFitData::DipoleFitData()
3338 : mri_head_t (NULL)
3339 , meg_head_t (NULL)
3340 , coord_frame(FIFFV_COORD_UNKNOWN)
3341 , nmeg (0)
3342 , neeg (0)
3343 , ch_names (NULL)
3344 , pick (NULL)
3345 , meg_coils (NULL)
3346 , eeg_els (NULL)
3347 , eeg_model (NULL)
3348 , bem_model (NULL)
3349 , sphere_funcs (NULL)
3350 , bem_funcs (NULL)
3351 , funcs (NULL)
3352 , mag_dipole_funcs (NULL)
3353 , fixed_noise (FALSE)
3354 , noise_orig (NULL)
3355 , noise (NULL)
3356 , nave (1)
3357 , proj (NULL)
3358 , column_norm (COLUMN_NORM_NONE)
3359 , fit_mag_dipoles (FALSE)
3360 , user (NULL)
3361 , user_free (NULL)
3362 {
3363  r0[0] = 0.0f;
3364  r0[1] = 0.0f;
3365  r0[2] = 0.0f;
3366 }
3367 
3368 //=============================================================================================================
3369 
3371 {
3372  if(mri_head_t)
3373  delete mri_head_t;
3374  if(meg_head_t)
3375  delete meg_head_t;
3376 
3377  if(meg_coils)
3378  delete meg_coils;
3379  if(eeg_els)
3380  delete eeg_els;
3381 
3382  bemname.clear();
3383 
3384  mne_free_cov_3(noise);
3385  mne_free_cov_3(noise_orig);
3386  ch_names.clear();
3387 
3388  if(pick)
3389  delete pick;
3390  if(bem_model)
3391  delete bem_model;
3392 
3393  if(eeg_model)
3394  delete eeg_model;
3395  if (user_free)
3396  user_free(user);
3397 
3398  if(proj)
3399  delete proj;
3400 
3401  free_dipole_fit_funcs(sphere_funcs);
3402  free_dipole_fit_funcs(bem_funcs);
3403  free_dipole_fit_funcs(mag_dipole_funcs);
3404 }
3405 
3406 //=============================================================================================================
3407 
3408 int DipoleFitData::setup_forward_model(DipoleFitData *d, MneCTFCompDataSet* comp_data, FwdCoilSet *comp_coils)
3409 /*
3410  * Take care of some hairy details
3411  */
3412 {
3413  FwdCompData* comp;
3414  dipoleFitFuncs f;
3415  int fit_sphere_to_bem = TRUE;
3416 
3417  if (!d->bemname.isEmpty()) {
3418  /*
3419  * Set up the boundary-element model
3420  */
3421  QString bemsolname = FwdBemModel::fwd_bem_make_bem_sol_name(d->bemname);
3422  d->bemname = bemsolname;
3423 
3424  printf("\nSetting up the BEM model using %s...\n",d->bemname.toUtf8().constData());
3425  printf("\nLoading surfaces...\n");
3426  d->bem_model = FwdBemModel::fwd_bem_load_three_layer_surfaces(d->bemname);
3427  if (d->bem_model) {
3428  printf("Three-layer model surfaces loaded.\n");
3429  }
3430  else {
3431  d->bem_model = FwdBemModel::fwd_bem_load_homog_surface(d->bemname);
3432  if (!d->bem_model)
3433  goto out;
3434  printf("Homogeneous model surface loaded.\n");
3435  }
3436  if (d->neeg > 0 && d->bem_model->nsurf == 1) {
3437  qCritical("Cannot use a homogeneous model in EEG calculations.");
3438  goto out;
3439  }
3440  printf("\nLoading the solution matrix...\n");
3441  if (FwdBemModel::fwd_bem_load_recompute_solution(d->bemname,FWD_BEM_UNKNOWN,FALSE,d->bem_model) == FAIL)
3442  goto out;
3443  printf("Employing the head->MRI coordinate transform with the BEM model.\n");
3444  if (FwdBemModel::fwd_bem_set_head_mri_t(d->bem_model,d->mri_head_t) == FAIL)
3445  goto out;
3446  printf("BEM model %s is now set up\n",d->bem_model->sol_name.toUtf8().constData());
3447  /*
3448  * Find the best-fitting sphere
3449  */
3450  if (fit_sphere_to_bem) {
3451  MneSurfaceOld* inner_skull;
3452  float simplex_size = 2e-2;
3453  float R;
3454 
3455  if ((inner_skull = d->bem_model->fwd_bem_find_surface(FIFFV_BEM_SURF_ID_BRAIN)) == NULL)
3456  goto out;
3457 
3458  if (fit_sphere_to_points(inner_skull->rr,inner_skull->np,simplex_size,d->r0,&R) == FAIL)
3459  goto out;
3460 
3461  FiffCoordTransOld::fiff_coord_trans(d->r0,d->mri_head_t,TRUE);
3462  printf("Fitted sphere model origin : %6.1f %6.1f %6.1f mm rad = %6.1f mm.\n",
3463  1000*d->r0[X_3],1000*d->r0[Y_3],1000*d->r0[Z_3],1000*R);
3464  }
3465  d->bem_funcs = f = new_dipole_fit_funcs();
3466  if (d->nmeg > 0) {
3467  /*
3468  * Use the new compensated field computation
3469  * It works the same way independent of whether or not the compensation is in effect
3470  */
3471  comp = FwdCompData::fwd_make_comp_data(comp_data,d->meg_coils,comp_coils,
3472  FwdBemModel::fwd_bem_field,NULL,NULL,d->bem_model,NULL);
3473  if (!comp)
3474  goto out;
3475  printf("Compensation setup done.\n");
3476 
3477  printf("MEG solution matrix...");
3478  if (FwdBemModel::fwd_bem_specify_coils(d->bem_model,d->meg_coils) == FAIL)
3479  goto out;
3480  if (FwdBemModel::fwd_bem_specify_coils(d->bem_model,comp->comp_coils) == FAIL)
3481  goto out;
3482  printf("[done]\n");
3483 
3484  f->meg_field = FwdCompData::fwd_comp_field;
3485  f->meg_vec_field = NULL;
3486  f->meg_client = comp;
3487  f->meg_client_free = FwdCompData::fwd_free_comp_data;
3488  }
3489  if (d->neeg > 0) {
3490  printf("\tEEG solution matrix...");
3491  if (FwdBemModel::fwd_bem_specify_els(d->bem_model,d->eeg_els) == FAIL)
3492  goto out;
3493  printf("[done]\n");
3494  f->eeg_pot = FwdBemModel::fwd_bem_pot_els;
3495  f->eeg_vec_pot = NULL;
3496  f->eeg_client = d->bem_model;
3497  }
3498  }
3499  if (d->neeg > 0 && !d->eeg_model) {
3500  qCritical("EEG sphere model not defined.");
3501  goto out;
3502  }
3503  d->sphere_funcs = f = new_dipole_fit_funcs();
3504  if (d->neeg > 0) {
3505  VEC_COPY_3(d->eeg_model->r0,d->r0);
3506  f->eeg_pot = FwdEegSphereModel::fwd_eeg_spherepot_coil;
3507  f->eeg_vec_pot = FwdEegSphereModel::fwd_eeg_spherepot_coil_vec;
3508  f->eeg_client = d->eeg_model;
3509  }
3510  if (d->nmeg > 0) {
3511  /*
3512  * Use the new compensated field computation
3513  * It works the same way independent of whether or not the compensation is in effect
3514  */
3515  comp = FwdCompData::fwd_make_comp_data(comp_data,d->meg_coils,comp_coils,
3516  FwdBemModel::fwd_sphere_field,
3517  FwdBemModel::fwd_sphere_field_vec,
3518  NULL,
3519  d->r0,NULL);
3520  if (!comp)
3521  goto out;
3522  f->meg_field = FwdCompData::fwd_comp_field;
3523  f->meg_vec_field = FwdCompData::fwd_comp_field_vec;
3524  f->meg_client = comp;
3525  f->meg_client_free = FwdCompData::fwd_free_comp_data;
3526  }
3527  printf("Sphere model origin : %6.1f %6.1f %6.1f mm.\n",
3528  1000*d->r0[X_3],1000*d->r0[Y_3],1000*d->r0[Z_3]);
3529  /*
3530  * Finally add the magnetic dipole fitting functions (for special purposes)
3531  */
3532  d->mag_dipole_funcs = f = new_dipole_fit_funcs();
3533  if (d->nmeg > 0) {
3534  /*
3535  * Use the new compensated field computation
3536  * It works the same way independent of whether or not the compensation is in effect
3537  */
3538  comp = FwdCompData::fwd_make_comp_data(comp_data,d->meg_coils,comp_coils,
3539  FwdBemModel::fwd_mag_dipole_field,
3540  FwdBemModel::fwd_mag_dipole_field_vec,
3541  NULL,
3542  NULL,NULL);
3543  if (!comp)
3544  goto out;
3545  f->meg_field = FwdCompData::fwd_comp_field;
3546  f->meg_vec_field = FwdCompData::fwd_comp_field_vec;
3547  f->meg_client = comp;
3548  f->meg_client_free = FwdCompData::fwd_free_comp_data;
3549  }
3550  f->eeg_pot = FwdBemModel::fwd_mag_dipole_field;
3551  f->eeg_vec_pot = FwdBemModel::fwd_mag_dipole_field_vec;
3552  /*
3553  * Select the appropriate fitting function
3554  */
3555  d->funcs = !d->bemname.isEmpty() ? d->bem_funcs : d->sphere_funcs;
3556 
3557  fprintf (stderr,"\n");
3558  return OK;
3559 
3560 out :
3561  return FAIL;
3562 }
3563 
3564 //=============================================================================================================
3565 
3566 MneCovMatrix* DipoleFitData::ad_hoc_noise(FwdCoilSet *meg, FwdCoilSet *eeg, float grad_std, float mag_std, float eeg_std)
3567 /*
3568  * Specify constant noise values
3569  */
3570 {
3571  int nchan;
3572  double *stds;
3573  QStringList names, ch_names;
3574  int k,n;
3575 
3576  printf("Using standard noise values "
3577  "(MEG grad : %6.1f fT/cm MEG mag : %6.1f fT EEG : %6.1f uV)\n",
3578  1e13*grad_std,1e15*mag_std,1e6*eeg_std);
3579 
3580  nchan = 0;
3581  if (meg)
3582  nchan = nchan + meg->ncoil;
3583  if (eeg)
3584  nchan = nchan + eeg->ncoil;
3585 
3586  stds = MALLOC_3(nchan,double);
3587 
3588  n = 0;
3589  if (meg) {
3590  for (k = 0; k < meg->ncoil; k++, n++) {
3591  if (meg->coils[k]->is_axial_coil()) {
3592  stds[n] = mag_std*mag_std;
3593 #ifdef TEST_REF
3594  if (meg->coils[k]->type == FIFFV_COIL_CTF_REF_MAG ||
3595  meg->coils[k]->type == FIFFV_COIL_CTF_REF_GRAD ||
3596  meg->coils[k]->type == FIFFV_COIL_CTF_OFFDIAG_REF_GRAD)
3597  stds[n] = 1e6*stds[n];
3598 #endif
3599  }
3600  else
3601  stds[n] = grad_std*grad_std;
3602  ch_names.append(meg->coils[k]->chname);
3603  }
3604  }
3605  if (eeg) {
3606  for (k = 0; k < eeg->ncoil; k++, n++) {
3607  stds[n] = eeg_std*eeg_std;
3608  ch_names.append(eeg->coils[k]->chname);
3609  }
3610  }
3611  names = ch_names;
3612  return MneCovMatrix::mne_new_cov(FIFFV_MNE_NOISE_COV,nchan,names,NULL,stds);
3613 }
3614 
3615 //=============================================================================================================
3616 
3617 int DipoleFitData::make_projection(const QList<QString> &projnames,
3618  const QList<FiffChInfo>& chs,
3619  int nch,
3620  MneProjOp* *res)
3621 /*
3622  * Process the projection data
3623  */
3624 {
3625  MneProjOp* all = NULL;
3626  MneProjOp* one = NULL;
3627  int k,found;
3628  int neeg;
3629 
3630  for (k = 0, neeg = 0; k < nch; k++)
3631  if (chs[k].kind == FIFFV_EEG_CH)
3632  neeg++;
3633 
3634  if (projnames.size() == 0 && neeg == 0)
3635  return OK;
3636 
3637  for (k = 0; k < projnames.size(); k++) {
3638  if ((one = mne_read_proj_op_3(projnames[k])) == NULL)
3639  goto bad;
3640  if (one->nitems == 0) {
3641  printf("No linear projection information in %s.\n",projnames[k].toUtf8().data());
3642  if(one)
3643  delete one;
3644  one = NULL;
3645  }
3646  else {
3647  printf("Loaded projection from %s:\n",projnames[k].toUtf8().data());
3648  mne_proj_op_report_3(stderr,"\t",one);
3649  all = MneProjOp::mne_proj_op_combine(all,one);
3650  if(one)
3651  delete one;
3652  one = NULL;
3653  }
3654  }
3655 
3656  if (neeg > 0) {
3657  found = FALSE;
3658  if (all) {
3659  for (k = 0; k < all->nitems; k++)
3660  if (all->items[k]->kind == FIFFV_MNE_PROJ_ITEM_EEG_AVREF) {
3661  found = TRUE;
3662  break;
3663  }
3664  }
3665  if (!found) {
3666  if ((one = MneProjOp::mne_proj_op_average_eeg_ref(chs,nch)) != NULL) {
3667  printf("Average EEG reference projection added:\n");
3668  mne_proj_op_report_3(stderr,"\t",one);
3669  all = MneProjOp::mne_proj_op_combine(all,one);
3670  if(one)
3671  delete one;
3672  one = NULL;
3673  }
3674  }
3675  }
3676  if (all && MneProjOp::mne_proj_op_affect_chs(all,chs,nch) == 0) {
3677  printf("Projection will not have any effect on selected channels. Projection omitted.\n");
3678  if(all)
3679  delete all;
3680  all = NULL;
3681  }
3682  *res = all;
3683  return OK;
3684 
3685 bad :
3686  if(all)
3687  delete all;
3688  all = NULL;
3689  return FAIL;
3690 }
3691 
3692 //=============================================================================================================
3693 
3694 int DipoleFitData::scale_noise_cov(DipoleFitData *f, int nave)
3695 {
3696  float nave_ratio = ((float)f->nave)/(float)nave;
3697  int k;
3698 
3699  if (!f->noise)
3700  return OK;
3701 
3702  if (f->noise->cov != NULL) {
3703  printf("Decomposing the sensor noise covariance matrix...\n");
3704  if (mne_decompose_eigen_cov_3(f->noise) == FAIL)
3705  goto bad;
3706 
3707  for (k = 0; k < f->noise->ncov*(f->noise->ncov+1)/2; k++)
3708  f->noise->cov[k] = nave_ratio*f->noise->cov[k];
3709  for (k = 0; k < f->noise->ncov; k++) {
3710  f->noise->lambda[k] = nave_ratio*f->noise->lambda[k];
3711  if (f->noise->lambda[k] < 0.0)
3712  f->noise->lambda[k] = 0.0;
3713  }
3714  if (mne_add_inv_cov_3(f->noise) == FAIL)
3715  goto bad;
3716  }
3717  else {
3718  for (k = 0; k < f->noise->ncov; k++)
3719  f->noise->cov_diag[k] = nave_ratio*f->noise->cov_diag[k];
3720  printf("Decomposition not needed for a diagonal noise covariance matrix.\n");
3721  if (mne_add_inv_cov_3(f->noise) == FAIL)
3722  goto bad;
3723  }
3724  printf("Effective nave is now %d\n",nave);
3725  f->nave = nave;
3726  return OK;
3727 
3728 bad :
3729  return FAIL;
3730 }
3731 
3732 //=============================================================================================================
3733 
3734 int DipoleFitData::scale_dipole_fit_noise_cov(DipoleFitData *f, int nave)
3735 {
3736  float nave_ratio = ((float)f->nave)/(float)nave;
3737  int k;
3738 
3739  if (!f->noise)
3740  return OK;
3741  if (f->fixed_noise)
3742  return OK;
3743 
3744  if (f->noise->cov) {
3745  /*
3746  * Do the decomposition and check that the matrix is positive definite
3747  */
3748  printf("Decomposing the noise covariance...");
3749  if (f->noise->cov) {
3750  if (mne_decompose_eigen_cov_3(f->noise) == FAIL)
3751  goto bad;
3752  for (k = 0; k < f->noise->ncov; k++) {
3753  if (f->noise->lambda[k] < 0.0)
3754  f->noise->lambda[k] = 0.0;
3755  }
3756  }
3757  for (k = 0; k < f->noise->ncov*(f->noise->ncov+1)/2; k++)
3758  f->noise->cov[k] = nave_ratio*f->noise->cov[k];
3759  for (k = 0; k < f->noise->ncov; k++) {
3760  f->noise->lambda[k] = nave_ratio*f->noise->lambda[k];
3761  if (f->noise->lambda[k] < 0.0)
3762  f->noise->lambda[k] = 0.0;
3763  }
3764  if (mne_add_inv_cov_3(f->noise) == FAIL)
3765  goto bad;
3766  }
3767  else {
3768  for (k = 0; k < f->noise->ncov; k++)
3769  f->noise->cov_diag[k] = nave_ratio*f->noise->cov_diag[k];
3770  printf("Decomposition not needed for a diagonal noise covariance matrix.\n");
3771  if (mne_add_inv_cov_3(f->noise) == FAIL)
3772  goto bad;
3773  }
3774  printf("Effective nave is now %d\n",nave);
3775  f->nave = nave;
3776  return OK;
3777 
3778 bad :
3779  return FAIL;
3780 }
3781 
3782 //=============================================================================================================
3783 
3784 int DipoleFitData::select_dipole_fit_noise_cov(DipoleFitData *f, mshMegEegData d)
3785 /*
3786  * Do the channel selection and scale with nave
3787  */
3788 {
3789  int nave,j,k;
3790  float nonsel_w = 30;
3791  int min_nchan = 20;
3792 
3793  if (!f || !f->noise_orig)
3794  return OK;
3795  if (!d)
3796  nave = 1;
3797  else {
3798  if (d->nave < 0)
3799  nave = d->meas->current->nave;
3800  else
3801  nave = d->nave;
3802  }
3803  /*
3804  * Channel selection
3805  */
3806  if (d) {
3807  float *w = MALLOC_3(f->noise_orig->ncov,float);
3808  int nomit_meg,nomit_eeg,nmeg,neeg;
3809  double *val;
3810 
3811  nmeg = neeg = 0;
3812  nomit_meg = nomit_eeg = 0;
3813  for (k = 0; k < f->noise_orig->ncov; k++) {
3814  if (f->noise_orig->ch_class[k] == MNE_COV_CH_EEG)
3815  neeg++;
3816  else
3817  nmeg++;
3818  if (is_selected_in_data(d,f->noise_orig->names[k]))
3819  w[k] = 1.0;
3820  else {
3821  w[k] = nonsel_w;
3822  if (f->noise_orig->ch_class[k] == MNE_COV_CH_EEG)
3823  nomit_eeg++;
3824  else
3825  nomit_meg++;
3826  }
3827  }
3828  mne_free_cov_3(f->noise); f->noise = NULL;
3829  if (nmeg > 0 && nmeg-nomit_meg > 0 && nmeg-nomit_meg < min_nchan) {
3830  qCritical("Too few MEG channels remaining");
3831  FREE_3(w);
3832  return FAIL;
3833  }
3834  if (neeg > 0 && neeg-nomit_eeg > 0 && neeg-nomit_eeg < min_nchan) {
3835  qCritical("Too few EEG channels remaining");
3836  FREE_3(w);
3837  return FAIL;
3838  }
3839  f->noise = MneCovMatrix::mne_dup_cov(f->noise_orig);
3840  if (nomit_meg+nomit_eeg > 0) {
3841  if (f->noise->cov) {
3842  for (j = 0; j < f->noise->ncov; j++)
3843  for (k = 0; k <= j; k++) {
3844  val = f->noise->cov+mne_lt_packed_index_3(j,k);
3845  *val = w[j]*w[k]*(*val);
3846  }
3847  }
3848  else {
3849  for (j = 0; j < f->noise->ncov; j++) {
3850  val = f->noise->cov_diag+j;
3851  *val = w[j]*w[j]*(*val);
3852  }
3853  }
3854  }
3855  FREE_3(w);
3856  }
3857  else {
3858  if (f->noise && f->nave == nave)
3859  return OK;
3860  f->noise = MneCovMatrix::mne_dup_cov(f->noise_orig);
3861  }
3862 
3863  return scale_dipole_fit_noise_cov(f,nave);
3864 }
3865 
3866 //=============================================================================================================
3867 
3869  const QString &measname,
3870  const QString& bemname,
3871  Vector3f *r0,
3873  int accurate_coils,
3874  const QString &badname,
3875  const QString &noisename,
3876  float grad_std,
3877  float mag_std,
3878  float eeg_std,
3879  float mag_reg,
3880  float grad_reg,
3881  float eeg_reg,
3882  int diagnoise,
3883  const QList<QString> &projnames,
3884  int include_meg,
3885  int include_eeg)
3886 /*
3887  * Background work for modelling
3888  */
3889 {
3890  DipoleFitData* res = new DipoleFitData;
3891  int k;
3892  QStringList badlist;
3893  int nbad = 0;
3894  QStringList file_bads;
3895  int file_nbad;
3896  int coord_frame = FIFFV_COORD_HEAD;
3897  MneCovMatrix* cov;
3898  FwdCoilSet* templates = NULL;
3899  MneCTFCompDataSet* comp_data = NULL;
3900  FwdCoilSet* comp_coils = NULL;
3901 
3902  /*
3903  * Read the coordinate transformations
3904  */
3905  if (!mriname.isEmpty()) {
3906  if ((res->mri_head_t = FiffCoordTransOld::mne_read_mri_transform(mriname)) == NULL)
3907  goto bad;
3908  }
3909  else if (!bemname.isEmpty()) {
3910  qWarning("Source of MRI / head transform required for the BEM model is missing");
3911  goto bad;
3912  }
3913  else {
3914  float move[] = { 0.0, 0.0, 0.0 };
3915  float rot[3][3] = { { 1.0, 0.0, 0.0 },
3916  { 0.0, 1.0, 0.0 },
3917  { 0.0, 0.0, 1.0 } };
3918  res->mri_head_t = new FiffCoordTransOld(FIFFV_COORD_MRI,FIFFV_COORD_HEAD,rot,move);
3919  }
3920 
3921  FiffCoordTransOld::mne_print_coord_transform(stdout,res->mri_head_t);
3922  if ((res->meg_head_t = FiffCoordTransOld::mne_read_meas_transform(measname)) == NULL)
3923  goto bad;
3924  FiffCoordTransOld::mne_print_coord_transform(stdout,res->meg_head_t);
3925  /*
3926  * Read the bad channel lists
3927  */
3928  if (!badname.isEmpty()) {
3929  if (mne_read_bad_channels_3(badname,badlist,nbad) != OK)
3930  goto bad;
3931  printf("%d bad channels read from %s.\n",nbad,badname.toUtf8().data());
3932  }
3933  if (mne_read_bad_channel_list_3(measname,file_bads,file_nbad) == OK && file_nbad > 0) {
3934  if (badlist.isEmpty())
3935  nbad = 0;
3936  for (k = 0; k < file_nbad; k++) {
3937  badlist.append(file_bads[k]);
3938  nbad++;
3939  }
3940  file_bads.clear();
3941  printf("%d bad channels read from the data file.\n",file_nbad);
3942  }
3943  printf("%d bad channels total.\n",nbad);
3944  /*
3945  * Read the channel information
3946  */
3947  if (read_meg_eeg_ch_info(measname,
3948  include_meg,
3949  include_eeg,
3950  badlist,
3951  nbad,
3952  res->chs,
3953  &res->nmeg,
3954  &res->neeg) != OK)
3955  goto bad;
3956 
3957  if (res->nmeg > 0)
3958  printf("Will use %3d MEG channels from %s\n",res->nmeg,measname.toUtf8().data());
3959  if (res->neeg > 0)
3960  printf("Will use %3d EEG channels from %s\n",res->neeg,measname.toUtf8().data());
3961  {
3962  QString s = mne_channel_names_to_string_3(res->chs,
3963  res->nmeg+res->neeg);
3964  int n;
3965  mne_string_to_name_list_3(s,res->ch_names,n);
3966  }
3967  /*
3968  * Make coil definitions
3969  */
3970  res->coord_frame = coord_frame;
3971  if (coord_frame == FIFFV_COORD_HEAD) {
3972  //#ifdef USE_SHARE_PATH
3973  // char *coilfile = mne_compose_mne_name("share/mne","coil_def.dat");
3974  //#else
3975  // const char *path = "setup/mne";
3976  // const char *filename = "coil_def.dat";
3977  // const char *coilfile = mne_compose_mne_name(path,filename);
3978 
3979  // QString qPath("/usr/pubsw/packages/mne/stable/share/mne/coil_def.dat");
3980 
3981  QString qPath = QString(QCoreApplication::applicationDirPath() + "/resources/general/coilDefinitions/coil_def.dat");
3982  QFile file(qPath);
3983  if ( !QCoreApplication::startingUp() )
3984  qPath = QCoreApplication::applicationDirPath() + QString("/resources/general/coilDefinitions/coil_def.dat");
3985  else if (!file.exists())
3986  qPath = "./resources/general/coilDefinitions/coil_def.dat";
3987 
3988  char *coilfile = MALLOC_3(strlen(qPath.toUtf8().data())+1,char);
3989  strcpy(coilfile,qPath.toUtf8().data());
3990  //#endif
3991 
3992  if (!coilfile)
3993  goto bad;
3994  if ((templates = FwdCoilSet::read_coil_defs(coilfile)) == NULL) {
3995  FREE_3(coilfile);
3996  goto bad;
3997  }
3998 
3999  if ((res->meg_coils = templates->create_meg_coils(res->chs,
4000  res->nmeg,
4001  accurate_coils ? FWD_COIL_ACCURACY_ACCURATE : FWD_COIL_ACCURACY_NORMAL,
4002  res->meg_head_t)) == NULL)
4003  goto bad;
4004  if ((res->eeg_els = FwdCoilSet::create_eeg_els(res->chs.mid(res->nmeg),
4005  res->neeg,
4006  NULL)) == NULL)
4007  goto bad;
4008  printf("Head coordinate coil definitions created.\n");
4009  }
4010  else {
4011  qWarning("Cannot handle computations in %s coordinates",mne_coord_frame_name_3(coord_frame));
4012  goto bad;
4013  }
4014  /*
4015  * Forward model setup
4016  */
4017  res->bemname = bemname;
4018  if (r0) {
4019  res->r0[0] = (*r0)[0];
4020  res->r0[1] = (*r0)[1];
4021  res->r0[2] = (*r0)[2];
4022  }
4023  res->eeg_model = eeg_model;
4024  /*
4025  * Compensation data
4026  */
4027  if ((comp_data = MneCTFCompDataSet::mne_read_ctf_comp_data(measname)) == NULL)
4028  goto bad;
4029  if (comp_data->ncomp > 0) { /* Compensation channel information may be needed */
4030  QList<FiffChInfo> comp_chs;
4031  int ncomp = 0;
4032 
4033  printf("%d compensation data sets in %s\n",comp_data->ncomp,measname.toUtf8().data());
4034  QList<FiffChInfo> temp;
4035  if (mne_read_meg_comp_eeg_ch_info_3(measname,
4036  temp,
4037  NULL,
4038  comp_chs,
4039  &ncomp,
4040  temp,
4041  NULL,
4042  NULL,
4043  NULL) == FAIL)
4044  goto bad;
4045  if (ncomp > 0) {
4046  if ((comp_coils = templates->create_meg_coils(comp_chs,
4047  ncomp,
4048  FWD_COIL_ACCURACY_NORMAL,
4049  res->meg_head_t)) == NULL) {
4050  goto bad;
4051  }
4052  printf("%d compensation channels in %s\n",comp_coils->ncoil,measname.toUtf8().data());
4053  }
4054  }
4055  else { /* Get rid of the empty data set */
4056  if(comp_data)
4057  delete comp_data;
4058  comp_data = NULL;
4059  }
4060  /*
4061  * Ready to set up the forward model
4062  */
4063  if (setup_forward_model(res,comp_data,comp_coils) == FAIL)
4064  goto bad;
4065  res->column_norm = COLUMN_NORM_LOC;
4066  /*
4067  * Projection data should go here
4068  */
4069  if (make_projection(projnames,
4070  res->chs,
4071  res->nmeg+res->neeg,
4072  &res->proj) == FAIL)
4073  goto bad;
4074  if (res->proj && res->proj->nitems > 0) {
4075  printf("Final projection operator is:\n");
4076  mne_proj_op_report_3(stderr,"\t",res->proj);
4077 
4078  if (mne_proj_op_chs_3(res->proj,res->ch_names,res->nmeg+res->neeg) == FAIL)
4079  goto bad;
4080  if (mne_proj_op_make_proj(res->proj) == FAIL)
4081  goto bad;
4082  }
4083  else
4084  printf("No projection will be applied to the data.\n");
4085 
4086  /*
4087  * Noise covariance
4088  */
4089  if (!noisename.isEmpty()) {
4090  if ((cov = mne_read_cov(noisename,FIFFV_MNE_SENSOR_COV)) == NULL)
4091  goto bad;
4092  printf("Read a %s noise-covariance matrix from %s\n",
4093  cov->cov_diag ? "diagonal" : "full", noisename.toUtf8().data());
4094  }
4095  else {
4096  if ((cov = ad_hoc_noise(res->meg_coils,res->eeg_els,grad_std,mag_std,eeg_std)) == NULL)
4097  goto bad;
4098  }
4099  res->noise = mne_pick_chs_cov_omit(cov,
4100  res->ch_names,
4101  res->nmeg+res->neeg,
4102  TRUE,
4103  res->chs);
4104  if (res->noise == NULL) {
4105  mne_free_cov_3(cov);
4106  goto bad;
4107  }
4108 
4109  printf("Picked appropriate channels from the noise-covariance matrix.\n");
4110  mne_free_cov_3(cov);
4111 
4112  /*
4113  * Apply the projection operator to the noise-covariance matrix
4114  */
4115  if (res->proj && res->proj->nitems > 0 && res->proj->nvec > 0) {
4116  if (mne_proj_op_apply_cov(res->proj,res->noise) == FAIL)
4117  goto bad;
4118  printf("Projection applied to the covariance matrix.\n");
4119  }
4120 
4121  /*
4122  * Force diagonal noise covariance?
4123  */
4124  if (diagnoise) {
4125  mne_revert_to_diag_cov(res->noise);
4126  printf("Using only the main diagonal of the noise-covariance matrix.\n");
4127  }
4128 
4129  /*
4130  * Regularize the possibly deficient noise-covariance matrix
4131  */
4132  if (res->noise->cov) {
4133  float regs[3];
4134  int do_it;
4135 
4136  regs[MNE_COV_CH_MEG_MAG] = mag_reg;
4137  regs[MNE_COV_CH_MEG_GRAD] = grad_reg;
4138  regs[MNE_COV_CH_EEG] = eeg_reg;
4139  /*
4140  * Classify the channels
4141  */
4142  if (mne_classify_channels_cov(res->noise,
4143  res->chs,
4144  res->nmeg+res->neeg) == FAIL)
4145  goto bad;
4146  /*
4147  * Do we need to do anything?
4148  */
4149  for (k = 0, do_it = 0; k < res->noise->ncov; k++) {
4150  if (res->noise->ch_class[k] != MNE_COV_CH_UNKNOWN &&
4151  regs[res->noise->ch_class[k]] > 0.0)
4152  do_it++;
4153  }
4154  /*
4155  * Apply regularization if necessary
4156  */
4157  if (do_it > 0)
4158  mne_regularize_cov(res->noise,regs);
4159  else
4160  printf("No regularization applied to the noise-covariance matrix\n");
4161  }
4162 
4163  /*
4164  * Do the decomposition and check that the matrix is positive definite
4165  */
4166  printf("Decomposing the noise covariance...\n");
4167  if (res->noise->cov) {
4168  if (mne_decompose_eigen_cov_3(res->noise) == FAIL)
4169  goto bad;
4170  printf("Eigenvalue decomposition done.\n");
4171  for (k = 0; k < res->noise->ncov; k++) {
4172  if (res->noise->lambda[k] < 0.0)
4173  res->noise->lambda[k] = 0.0;
4174  }
4175  }
4176  else {
4177  printf("Decomposition not needed for a diagonal covariance matrix.\n");
4178  if (mne_add_inv_cov_3(res->noise) == FAIL)
4179  goto bad;
4180  }
4181 
4182  badlist.clear();
4183  delete templates;
4184  delete comp_coils;
4185  if(comp_data)
4186  delete comp_data;
4187  return res;
4188 
4189 bad : {
4190  badlist.clear();
4191  delete templates;
4192  delete comp_coils;
4193  if(comp_data)
4194  delete comp_data;
4195  if(res)
4196  delete res;
4197  return NULL;
4198  }
4199 }
4200 
4201 //=============================================================================================================
4202 
4203 //============================= dipole_forward.c =============================
4204 
4205 void print_fields(float *rd,
4206  float *Q,
4207  float time,
4208  float integ,
4209  DipoleFitData* fit,
4210  MneMeasData* data)
4211 
4212 {
4213  float *one = MALLOC_3(data->nchan,float);
4214  int k;
4215  float **fwd = NULL;
4216 
4217  if (mne_get_values_from_data_3(time,integ,data->current->data,data->current->np,data->nchan,data->current->tmin,
4218  1.0/data->current->tstep,FALSE,one) == FAIL) {
4219  printf("Cannot pick time: %7.1f ms\n",1000*time);
4220  return;
4221  }
4222  for (k = 0; k < data->nchan; k++)
4223  if (data->chs[k].chpos.coil_type == FIFFV_COIL_CTF_REF_GRAD ||
4224  data->chs[k].chpos.coil_type == FIFFV_COIL_CTF_OFFDIAG_REF_GRAD) {
4225  printf("%g ",1e15*one[k]);
4226  }
4227  printf("\n");
4228 
4229  fwd = ALLOC_CMATRIX_3(3,fit->nmeg+fit->neeg);
4230  if (DipoleFitData::compute_dipole_field(fit,rd,FALSE,fwd) == FAIL)
4231  goto out;
4232 
4233  for (k = 0; k < data->nchan; k++)
4234  if (data->chs[k].chpos.coil_type == FIFFV_COIL_CTF_REF_GRAD ||
4235  data->chs[k].chpos.coil_type == FIFFV_COIL_CTF_OFFDIAG_REF_GRAD) {
4236  printf("%g ",1e15*(Q[X_3]*fwd[X_3][k]+Q[Y_3]*fwd[Y_3][k]+Q[Z_3]*fwd[Z_3][k]));
4237  }
4238  printf("\n");
4239 
4240 out : {
4241  FREE_3(one);
4242  FREE_CMATRIX_3(fwd);
4243  }
4244  return;
4245 }
4246 
4247 //=============================================================================================================
4248 
4249 DipoleForward* dipole_forward(DipoleFitData* d,
4250  float **rd,
4251  int ndip,
4252  DipoleForward* old)
4253 /*
4254  * Compute the forward solution and do other nice stuff
4255  */
4256 {
4257  DipoleForward* res;
4258  float **this_fwd;
4259  float S[3];
4260  int k,p;
4261  /*
4262  * Allocate data if necessary
4263  */
4264  if (old && old->ndip == ndip && old->nch == d->nmeg+d->neeg) {
4265  res = old;
4266  }
4267  else {
4268  delete old; old = NULL;
4269  res = new DipoleForward;
4270  res->fwd = ALLOC_CMATRIX_3(3*ndip,d->nmeg+d->neeg);
4271  res->uu = ALLOC_CMATRIX_3(3*ndip,d->nmeg+d->neeg);
4272  res->vv = ALLOC_CMATRIX_3(3*ndip,3);
4273  res->sing = MALLOC_3(3*ndip,float);
4274  res->nch = d->nmeg+d->neeg;
4275  res->rd = ALLOC_CMATRIX_3(ndip,3);
4276  res->scales = MALLOC_3(3*ndip,float);
4277  res->ndip = ndip;
4278  }
4279 
4280  for (k = 0; k < ndip; k++) {
4281  VEC_COPY_3(res->rd[k],rd[k]);
4282  this_fwd = res->fwd + 3*k;
4283  /*
4284  * Calculate the field of three orthogonal dipoles
4285  */
4286  if ((DipoleFitData::compute_dipole_field(d,rd[k],TRUE,this_fwd)) == FAIL)
4287  goto bad;
4288  /*
4289  * Choice of column normalization
4290  * (componentwise normalization is not recommended)
4291  */
4292  if (d->column_norm == COLUMN_NORM_LOC || d->column_norm == COLUMN_NORM_COMP) {
4293  for (p = 0; p < 3; p++)
4294  S[p] = mne_dot_vectors_3(res->fwd[3*k+p],res->fwd[3*k+p],res->nch);
4295  if (d->column_norm == COLUMN_NORM_COMP) {
4296  for (p = 0; p < 3; p++)
4297  res->scales[3*k+p] = sqrt(S[p]);
4298  }
4299  else {
4300  /*
4301  * Divide by three or not?
4302  */
4303  res->scales[3*k+0] = res->scales[3*k+1] = res->scales[3*k+2] = sqrt(S[0]+S[1]+S[2])/3.0;
4304  }
4305  for (p = 0; p < 3; p++) {
4306  if (res->scales[3*k+p] > 0.0) {
4307  res->scales[3*k+p] = 1.0/res->scales[3*k+p];
4308  mne_scale_vector_3(res->scales[3*k+p],res->fwd[3*k+p],res->nch);
4309  }
4310  else
4311  res->scales[3*k+p] = 1.0;
4312  }
4313  }
4314  else {
4315  res->scales[3*k] = 1.0;
4316  res->scales[3*k+1] = 1.0;
4317  res->scales[3*k+2] = 1.0;
4318  }
4319  }
4320 
4321  /*
4322  * SVD
4323  */
4324  if (mne_svd_3(res->fwd,3*ndip,d->nmeg+d->neeg,res->sing,res->vv,res->uu) != 0)
4325  goto bad;
4326 
4327  return res;
4328 
4329 bad : {
4330  if (!old)
4331  delete res;
4332  return NULL;
4333  }
4334 }
4335 
4336 //=============================================================================================================
4337 
4338 DipoleForward* DipoleFitData::dipole_forward_one(DipoleFitData* d,
4339  float *rd,
4340  DipoleForward* old)
4341 /*
4342  * Convenience function to compute the field of one dipole
4343  */
4344 {
4345  float *rds[1];
4346  rds[0] = rd;
4347  return dipole_forward(d,rds,1,old);
4348 }
4349 
4350 //=============================================================================================================
4351 // fit_dipoles.c
4352 static float fit_eval(float *rd,int npar,void *user)
4353 /*
4354  * Calculate the residual sum of squares
4355  */
4356 {
4357  DipoleFitData* fit = (DipoleFitData*)user;
4358  DipoleForward* fwd;
4359  fitDipUser fuser = (fitDipUser)fit->user;
4360  double Bm2,one;
4361  int ncomp,c;
4362  (void)npar; // avoid unused variable warning.
4363 
4364  fwd = fuser->fwd = DipoleFitData::dipole_forward_one(fit,rd,fuser->fwd);
4365  ncomp = fwd->sing[2]/fwd->sing[0] > fuser->limit ? 3 : 2;
4366  if (fuser->report_dim)
4367  printf("ncomp = %d\n",ncomp);
4368 
4369  for (c = 0, Bm2 = 0.0; c < ncomp; c++) {
4370  one = mne_dot_vectors_3(fwd->uu[c],fuser->B,fwd->nch);
4371  Bm2 = Bm2 + one*one;
4372  }
4373  return fuser->B2-Bm2;
4374 }
4375 
4376 static int find_best_guess(float *B, /* The whitened data */
4377  int nch,
4378  GuessData* guess, /* Guesses */
4379  float limit, /* Pseudoradial component omission limit */
4380  int *bestp, /* Which is the best */
4381  float *goodp) /* Best goodness of fit */
4382 /*
4383  * Thanks to the precomputed SVD everything is really simple
4384  */
4385 {
4386  int k,c;
4387  double B2,Bm2,this_good,one;
4388  int best = -1;
4389  float good = 0.0;
4390  DipoleForward* fwd;
4391  int ncomp;
4392 
4393  B2 = mne_dot_vectors_3(B,B,nch);
4394  for (k = 0; k < guess->nguess; k++) {
4395  fwd = guess->guess_fwd[k];
4396  if (fwd->nch == nch) {
4397  ncomp = fwd->sing[2]/fwd->sing[0] > limit ? 3 : 2;
4398  for (c = 0, Bm2 = 0.0; c < ncomp; c++) {
4399  one = mne_dot_vectors_3(fwd->uu[c],B,nch);
4400  Bm2 = Bm2 + one*one;
4401  }
4402  this_good = 1.0 - (B2 - Bm2)/B2;
4403  if (this_good > good) {
4404  best = k;
4405  good = this_good;
4406  }
4407  }
4408  }
4409  if (best < 0) {
4410  printf("No reasonable initial guess found.");
4411  return FAIL;
4412  }
4413  *bestp = best;
4414  *goodp = good;
4415  return OK;
4416 }
4417 
4418 static float **make_initial_dipole_simplex(float *r0,
4419  float size)
4420 /*
4421  * Make the initial tetrahedron
4422  */
4423 {
4424  /*
4425  * For this definition of a regular tetrahedron, see
4426  *
4427  * http://mathworld.wolfram.com/Tetrahedron.html
4428  *
4429  */
4430  float x = sqrt(3.0)/3.0;
4431  float r = sqrt(6.0)/12.0;
4432  float R = 3*r;
4433  float d = x/2.0;
4434  float rr[][3] = { { x , 0.0, -r },
4435  { -d, 0.5, -r },
4436  { -d, -0.5, -r },
4437  { 0.0, 0.0, R } };
4438 
4439  float **simplex = ALLOC_CMATRIX_3(4,3);
4440  int j,k;
4441 
4442  for (j = 0; j < 4; j++) {
4443  VEC_COPY_3(simplex[j],rr[j]);
4444  for (k = 0; k < 3; k++)
4445  simplex[j][k] = size*simplex[j][k] + r0[k];
4446  }
4447  return simplex;
4448 }
4449 
4450 static int report_func(int loop,
4451  float *fitpar,
4452  int npar,
4453  double fval_lo,
4454  double fval_hi,
4455  double par_diff)
4456 /*
4457  * Report periodically
4458  */
4459 {
4460  float *r0 = fitpar;
4461  (void) npar; // avoid unused variable warning.
4462  fprintf(stdout,"loop %d rd %7.2f %7.2f %7.2f fval %g %g par diff %g\n",
4463  loop,1000*r0[0],1000*r0[1],1000*r0[2],fval_lo,fval_hi,1000*par_diff);
4464 
4465  return OK;
4466 }
4467 
4468 static int fit_Q(DipoleFitData* fit, /* The fit data */
4469  float *B, /* Measurement */
4470  float *rd, /* Dipole position */
4471  float limit, /* Radial component omission limit */
4472  float *Q, /* The result */
4473  int *ncomp,
4474  float *res) /* Residual sum of squares */
4475 /*
4476  * fit the dipole moment once the location is known
4477  */
4478 {
4479  int c;
4480  DipoleForward* fwd = DipoleFitData::dipole_forward_one(fit,rd,NULL);
4481  float Bm2,one;
4482 
4483  if (!fwd)
4484  return FAIL;
4485 
4486  *ncomp = fwd->sing[2]/fwd->sing[0] > limit ? 3 : 2;
4487 
4488  Q[0] = Q[1] = Q[2] = 0.0;
4489  for (c = 0, Bm2 = 0.0; c < *ncomp; c++) {
4490  one = mne_dot_vectors_3(fwd->uu[c],B,fwd->nch);
4491  mne_add_scaled_vector_to_3(fwd->vv[c],one/fwd->sing[c],Q,3);
4492  Bm2 = Bm2 + one*one;
4493  }
4494  /*
4495  * Counteract the effect of column normalization
4496  */
4497  for (c = 0; c < 3; c++)
4498  Q[c] = fwd->scales[c]*Q[c];
4499  *res = mne_dot_vectors_3(B,B,fwd->nch) - Bm2;
4500 
4501  delete fwd;
4502 
4503  return OK;
4504 }
4505 
4506 static float rtol(float *vals,int nval)
4507 
4508 {
4509  float minv,maxv;
4510  int k;
4511 
4512  minv = maxv = vals[0];
4513  for (k = 1; k < nval; k++) {
4514  if (vals[k] < minv)
4515  minv = vals[k];
4516  if (vals[k] > maxv)
4517  maxv = vals[k];
4518  }
4519  return 2.0*(maxv-minv)/(maxv+minv);
4520 }
4521 
4522 /*
4523  * This routine comes from Numerical recipes
4524  */
4525 
4526 #define ALPHA 1.0
4527 #define BETA 0.5
4528 #define GAMMA 2.0
4529 #define MIN_STOL_LOOP 5
4530 
4531 static float tryf (float **p,
4532  float *y,
4533  float *psum,
4534  int ndim,
4535  float (*func)(float *x,int npar,void *user_data), /* The function to be evaluated */
4536  void *user_data, /* Data to be passed to the above function in each evaluation */
4537  int ihi,
4538  int *neval,
4539  float fac)
4540 
4541 {
4542  int j;
4543  float fac1,fac2,ytry,*ptry;
4544 
4545  ptry = ALLOC_FLOAT_3(ndim);
4546  fac1 = (1.0-fac)/ndim;
4547  fac2 = fac1-fac;
4548  for (j = 0; j < ndim; j++)
4549  ptry[j] = psum[j]*fac1-p[ihi][j]*fac2;
4550  ytry = (*func)(ptry,ndim,user_data);
4551  ++(*neval);
4552  if (ytry < y[ihi]) {
4553  y[ihi] = ytry;
4554  for (j = 0; j < ndim; j++) {
4555  psum[j] += ptry[j]-p[ihi][j];
4556  p[ihi][j] = ptry[j];
4557  }
4558  }
4559  FREE_3(ptry);
4560  return ytry;
4561 }
4562 
4563 int simplex_minimize(float **p, /* The initial simplex */
4564  float *y, /* Function values at the vertices */
4565  int ndim, /* Number of variables */
4566  float ftol, /* Relative convergence tolerance */
4567  float stol,
4568  float (*func)(float *x,int npar,void *user_data),/* The function to be evaluated */
4569  void *user_data, /* Data to be passed to the above function in each evaluation */
4570  int max_eval, /* Maximum number of function evaluations */
4571  int *neval, /* Number of function evaluations */
4572  int report, /* How often to report (-1 = no_reporting) */
4573  int (*report_func)(int loop,
4574  float *fitpar, int npar,
4575  double fval_lo,
4576  double fval_hi,
4577  double par_diff)) /* The function to be called when reporting */
4578 /*
4579  * Minimization with the simplex algorithm
4580  * Modified from Numerical recipes
4581  */
4582 
4583 {
4584  int i,j,ilo,ihi,inhi;
4585  int mpts = ndim+1;
4586  float ytry,ysave,sum,rtol,*psum;
4587  double dsum,diff;
4588  int result = 0;
4589  int count = 0;
4590  int loop = 1;
4591 
4592  psum = ALLOC_FLOAT_3(ndim);
4593  *neval = 0;
4594  for (j = 0; j < ndim; j++) {
4595  for (i = 0,sum = 0.0; i<mpts; i++)
4596  sum += p[i][j];
4597  psum[j] = sum;
4598  }
4599  if (report_func != NULL && report > 0)
4600  (void)report_func (0,p[0],ndim,-1.0,-1.0,0.0);
4601 
4602  dsum = 0.0;
4603  for (;;count++,loop++) {
4604  ilo = 1;
4605  ihi = y[1]>y[2] ? (inhi = 2,1) : (inhi = 1,2);
4606  for (i = 0; i < mpts; i++) {
4607  if (y[i] < y[ilo]) ilo = i;
4608  if (y[i] > y[ihi]) {
4609  inhi = ihi;
4610  ihi = i;
4611  } else if (y[i] > y[inhi])
4612  if (i != ihi) inhi = i;
4613  }
4614  rtol = 2.0*std::fabs(y[ihi]-y[ilo])/(std::fabs(y[ihi])+std::fabs(y[ilo]));
4615  /*
4616  * Report that we are proceeding...
4617  */
4618  if (count == report && report_func != NULL) {
4619  if (report_func (loop,p[ilo],ndim,y[ilo],y[ihi],sqrt(dsum))) {
4620  printf("Interation interrupted.");
4621  result = -1;
4622  break;
4623  }
4624  count = 0;
4625  }
4626  if (rtol < ftol) break;
4627  if (*neval >= max_eval) {
4628  printf("Maximum number of evaluations exceeded.");
4629  result = -1;
4630  break;
4631  }
4632  if (stol > 0) { /* Has the simplex collapsed? */
4633  for (dsum = 0.0, j = 0; j < ndim; j++) {
4634  diff = p[ilo][j] - p[ihi][j];
4635  dsum += diff*diff;
4636  }
4637  if (loop > MIN_STOL_LOOP && sqrt(dsum) < stol)
4638  break;
4639  }
4640  ytry = tryf(p,y,psum,ndim,func,user_data,ihi,neval,-ALPHA);
4641  if (ytry <= y[ilo])
4642  ytry = tryf(p,y,psum,ndim,func,user_data,ihi,neval,GAMMA);
4643  else if (ytry >= y[inhi]) {
4644  ysave = y[ihi];
4645  ytry = tryf(p,y,psum,ndim,func,user_data,ihi,neval,BETA);
4646  if (ytry >= ysave) {
4647  for (i = 0; i < mpts; i++) {
4648  if (i != ilo) {
4649  for (j = 0; j < ndim; j++) {
4650  psum[j] = 0.5*(p[i][j]+p[ilo][j]);
4651  p[i][j] = psum[j];
4652  }
4653  y[i] = (*func)(psum,ndim,user_data);
4654  }
4655  }
4656  *neval += ndim;
4657  for (j = 0; j < ndim; j++) {
4658  for (i = 0,sum = 0.0; i < mpts; i++)
4659  sum += p[i][j];
4660  psum[j] = sum;
4661  }
4662  }
4663  }
4664  }
4665  FREE_3 (psum);
4666  return (result);
4667 }
4668 
4669 //=============================================================================================================
4670 // fit_dipoles.c
4671 bool DipoleFitData::fit_one(DipoleFitData* fit, /* Precomputed fitting data */
4672  GuessData* guess, /* The initial guesses */
4673  float time, /* Which time is it? */
4674  float *B, /* The field to fit */
4675  int verbose,
4676  ECD& res /* The fitted dipole */
4677  )
4678 {
4679  float **simplex = NULL; /* The simplex */
4680  float vals[4]; /* Values at the vertices */
4681  float limit = 0.2; /* (pseudo) radial component omission limit */
4682  float size = 1e-2; /* Size of the initial simplex */
4683  float ftol[] = { 1e-2, 1e-2 }; /* Tolerances on the the two passes */
4684  float atol[] = { 0.2e-3, 0.2e-3 }; /* If dipole movement between two iterations is less than this,
4685  we consider to have converged */
4686  int ntol = 2;
4687  int max_eval = 1000; /* Limit for fit function evaluations */
4688  int report_interval = verbose ? 1 : -1; /* How often to report the intermediate result */
4689 
4690  int best;
4691  float good,rd_guess[3],rd_final[3],Q[3],final_val;
4693  int k,p,neval,neval_tot,nchan,ncomp;
4694  int fit_fail;
4695 
4696  nchan = fit->nmeg+fit->neeg;
4697  user.fwd = NULL;
4698 
4699  if (MneProjOp::mne_proj_op_proj_vector(fit->proj,B,nchan,TRUE) == FAIL)
4700  goto bad;
4701 
4702  if (mne_whiten_one_data(B,B,nchan,fit->noise) == FAIL)
4703  goto bad;
4704  /*
4705  * Get the initial guess
4706  */
4707  if (find_best_guess(B,nchan,guess,limit,&best,&good) < 0)
4708  goto bad;
4709 
4710  user.limit = limit;
4711  user.B = B;
4712  user.B2 = mne_dot_vectors_3(B,B,nchan);
4713  user.fwd = NULL;
4714  user.report_dim = FALSE;
4715  fit->user = &user;
4716 
4717  VEC_COPY_3(rd_guess,guess->rr[best]);
4718  VEC_COPY_3(rd_final,guess->rr[best]);
4719 
4720  neval_tot = 0;
4721  fit_fail = FALSE;
4722  for (k = 0; k < ntol; k++) {
4723  /*
4724  * Do first pass with the sphere model
4725  */
4726  if (k == 0)
4727  fit->funcs = fit->sphere_funcs;
4728  else
4729  fit->funcs = !fit->bemname.isEmpty() ? fit->bem_funcs : fit->sphere_funcs;
4730 
4731  simplex = make_initial_dipole_simplex(rd_guess,size);
4732  for (p = 0; p < 4; p++)
4733  vals[p] = fit_eval(simplex[p],3,fit);
4734  if (simplex_minimize(simplex, /* The initial simplex */
4735  vals, /* Function values at the vertices */
4736  3, /* Number of variables */
4737  ftol[k], /* Relative convergence tolerance for the target function */
4738  atol[k], /* Absolute tolerance for the change in the parameters */
4739  fit_eval, /* The function to be evaluated */
4740  fit, /* Data to be passed to the above function in each evaluation */
4741  max_eval, /* Maximum number of function evaluations */
4742  &neval, /* Number of function evaluations */
4743  report_interval, /* How often to report (-1 = no_reporting) */
4744  report_func) != OK) {
4745  if (k == 0)
4746  goto bad;
4747  else {
4748  printf("\nWarning (t = %8.1f ms) : g = %6.1f %% final val = %7.3f rtol = %f\n",
4749  1000*time,100*(1 - vals[0]/user.B2),vals[0],rtol(vals,4));
4750  fit_fail = TRUE;
4751  }
4752  }
4753  VEC_COPY_3(rd_final,simplex[0]);
4754  VEC_COPY_3(rd_guess,simplex[0]);
4755  FREE_CMATRIX_3(simplex); simplex = NULL;
4756 
4757  neval_tot += neval;
4758  final_val = vals[0];
4759  }
4760  /*
4761  * Confidence limits should be computed here
4762  */
4763  /*
4764  * Compute the dipole moment at the final point
4765  */
4766  if (fit_Q(fit,user.B,rd_final,user.limit,Q,&ncomp,&final_val) == OK) {
4767  res.time = time;
4768  res.valid = true;
4769  for(int i = 0; i < 3; ++i)
4770  res.rd[i] = rd_final[i];
4771  for(int i = 0; i < 3; ++i)
4772  res.Q[i] = Q[i];
4773  res.good = 1.0 - final_val/user.B2;
4774  if (fit_fail)
4775  res.good = -res.good;
4776  res.khi2 = final_val;
4777  if (fit->proj)
4778  res.nfree = nchan-3-ncomp-fit->proj->nvec;
4779  else
4780  res.nfree = nchan-3-ncomp;
4781  res.neval = neval_tot;
4782  }
4783  else
4784  goto bad;
4785  delete user.fwd;
4786  FREE_CMATRIX_3(simplex);
4787 
4788  return true;
4789 
4790 bad : {
4791  delete user.fwd;
4792  FREE_CMATRIX_3(simplex);
4793  return false;
4794  }
4795 }
4796 
4797 //=============================================================================================================
4798 
4799 int DipoleFitData::compute_dipole_field(DipoleFitData* d, float *rd, int whiten, float **fwd)
4800 /*
4801  * Compute the field and take whitening and projection into account
4802  */
4803 {
4804  float *eeg_fwd[3];
4805  static float Qx[] = {1.0,0.0,0.0};
4806  static float Qy[] = {0.0,1.0,0.0};
4807  static float Qz[] = {0.0,0.0,1.0};
4808  int k;
4809  /*
4810  * Compute the fields
4811  */
4812  if (d->nmeg > 0) {
4813  if (d->funcs->meg_vec_field) {
4814  if (d->funcs->meg_vec_field(rd,d->meg_coils,fwd,d->funcs->meg_client) != OK)
4815  goto bad;
4816  }
4817  else {
4818  if (d->funcs->meg_field(rd,Qx,d->meg_coils,fwd[0],d->funcs->meg_client) != OK)
4819  goto bad;
4820  if (d->funcs->meg_field(rd,Qy,d->meg_coils,fwd[1],d->funcs->meg_client) != OK)
4821  goto bad;
4822  if (d->funcs->meg_field(rd,Qz,d->meg_coils,fwd[2],d->funcs->meg_client) != OK)
4823  goto bad;
4824  }
4825  }
4826 
4827  if (d->neeg > 0) {
4828  if (d->funcs->eeg_vec_pot) {
4829  eeg_fwd[0] = fwd[0]+d->nmeg;
4830  eeg_fwd[1] = fwd[1]+d->nmeg;
4831  eeg_fwd[2] = fwd[2]+d->nmeg;
4832  if (d->funcs->eeg_vec_pot(rd,d->eeg_els,eeg_fwd,d->funcs->eeg_client) != OK)
4833  goto bad;
4834  }
4835  else {
4836  if (d->funcs->eeg_pot(rd,Qx,d->eeg_els,fwd[0]+d->nmeg,d->funcs->eeg_client) != OK)
4837  goto bad;
4838  if (d->funcs->eeg_pot(rd,Qy,d->eeg_els,fwd[1]+d->nmeg,d->funcs->eeg_client) != OK)
4839  goto bad;
4840  if (d->funcs->eeg_pot(rd,Qz,d->eeg_els,fwd[2]+d->nmeg,d->funcs->eeg_client) != OK)
4841  goto bad;
4842  }
4843  }
4844 
4845  /*
4846  * Apply projection
4847  */
4848 #ifdef DEBUG
4849  fprintf(stdout,"orig : ");
4850  for (k = 0; k < 3; k++)
4851  fprintf(stdout,"%g ",sqrt(mne_dot_vectors_3(fwd[k],fwd[k],d->nmeg+d->neeg)));
4852  fprintf(stdout,"\n");
4853 #endif
4854 
4855  for (k = 0; k < 3; k++)
4856  if (MneProjOp::mne_proj_op_proj_vector(d->proj,fwd[k],d->nmeg+d->neeg,TRUE) == FAIL)
4857  goto bad;
4858 
4859 #ifdef DEBUG
4860  fprintf(stdout,"proj : ");
4861  for (k = 0; k < 3; k++)
4862  fprintf(stdout,"%g ",sqrt(mne_dot_vectors_3(fwd[k],fwd[k],d->nmeg+d->neeg)));
4863  fprintf(stdout,"\n");
4864 #endif
4865 
4866  /*
4867  * Whiten
4868  */
4869  if (d->noise && whiten) {
4870  if (mne_whiten_data(fwd,fwd,3,d->nmeg+d->neeg,d->noise) == FAIL)
4871  goto bad;
4872  }
4873 
4874 #ifdef DEBUG
4875  fprintf(stdout,"white : ");
4876  for (k = 0; k < 3; k++)
4877  fprintf(stdout,"%g ",sqrt(mne_dot_vectors_3(fwd[k],fwd[k],d->nmeg+d->neeg)));
4878  fprintf(stdout,"\n");
4879 #endif
4880 
4881  return OK;
4882 
4883 bad :
4884  return FAIL;
4885 }
QSharedPointer< FiffStream > SPtr
Definition: fiff_stream.h:107
#define FIFFV_MNE_NOISE_COV
#define FIFFV_MNE_COORD_MRI_VOXEL
GuessData class declaration.
static bool fit_one(DipoleFitData *fit, GuessData *guess, float time, float *B, int verbose, ECD &res)
FIFFLIB::FiffSparseMatrix * pick
FWDLIB::FwdCoilSet * meg_coils
FiffStream class declaration.
#define FIFFV_MNE_COORD_CTF_DEVICE
This structure is used in the compensated field calculations.
Definition: fwd_comp_data.h:87
QSharedPointer< FiffDirNode > SPtr
Definition: fiff_dir_node.h:77
int nfree
Definition: ecd.h:113
static DipoleFitData * setup_dipole_fit_data(const QString &mriname, const QString &measname, const QString &bemname, Eigen::Vector3f *r0, FWDLIB::FwdEegSphereModel *eeg_model, int accurate_coils, const QString &badname, const QString &noisename, float grad_std, float mag_std, float eeg_std, float mag_reg, float grad_reg, float eeg_reg, int diagnoise, const QList< QString > &projnames, int include_meg, int include_eeg)
Eigen::Vector3f r0
Definition: fiff_ch_pos.h:112
Eigen::Vector3f rd
Definition: ecd.h:109
Covariance matrix storage.
GuessData description.
Definition: guess_data.h:78
#define FIFFV_COIL_CTF_REF_GRAD
This defines a surface.
Electric Current Dipole description.
Definition: ecd.h:72
One MNE CTF Compensation Data Set description.
Dipole Fit Data implementation.
#define FIFF_MNE_COV_DIAG
MNELIB::MneCovMatrix * noise
MneCovMatrix class declaration.
QString chname
Definition: fwd_coil.h:167
MNELIB::MneCovMatrix * noise_orig
#define FIFF_MNE_COV_NFREE
#define FIFF_CH_INFO
Definition: fiff_file.h:456
Matrix specification with a channel list.
float time
Definition: ecd.h:108
#define FIFF_NCHAN
Definition: fiff_file.h:453
FIFFLIB::fiff_float_t * data
FIFFLIB::fiff_int_t coding
#define FIFF_MNE_COV_DIM
FWDLIB::FwdBemModel * bem_model
FWDLIB::FwdEegSphereModel * eeg_model
#define FIFFV_MNE_COORD_FS_TAL_LTZ
#define FIFF_MNE_COV
bool valid
Definition: ecd.h:107
MneSurfaceOld class declaration.
MNEProjItem class declaration.
MNE SSS Data description.
Definition: mne_sss_data.h:82
#define FIFF_MNE_PROJ_ITEM_ACTIVE
FIFFLIB::FiffCoordTransOld * mri_head_t
One linear projection item.
Definition: mne_proj_op.h:83
#define FIFF_NAME
Definition: fiff_file.h:485
float good
Definition: ecd.h:111
FwdBemModel class declaration.
float khi2
Definition: ecd.h:112
FwdCoilSet description.
Definition: fwd_coil_set.h:75
FwdCompData class declaration.
FwdCoilSet * create_meg_coils(const QList< FIFFLIB::FiffChInfo > &chs, int nch, int acc, const FIFFLIB::FiffCoordTransOld *t)
dipoleFitFuncs mag_dipole_funcs
Data associated with MNE computations for each mneMeasDataSet.
#define FIFFV_COIL_CTF_REF_MAG
Electric Current Dipole description.
FIFF File I/O routines.
Definition: fiff_stream.h:104
Coordinate transformation descriptor.
One linear projection item.
Definition: mne_proj_item.h:77
int neval
Definition: ecd.h:114
#define FIFF_MNE_COV_EIGENVALUES
Electric Current Dipole (ECD) class declaration.
FIFFLIB::fiff_int_t * inds
#define FIFF_MNE_COV_KIND
#define FIFF_COORD_TRANS
Definition: fiff_file.h:475
ToDo Old implementation use new fiff_id.h instead.
Definition: fiff_types.h:218
Channel info descriptor.
Definition: fiff_ch_info.h:74
#define FIFFV_MNE_COORD_RAS
Eigen::Vector3f Q
Definition: ecd.h:110
#define FIFFV_SSS_JOB_NOTHING
Definition: fiff_file.h:540
bool is_axial_coil() const
Definition: fwd_coil.cpp:268
static QStringList split_name_list(QString p_sNameList)
FIFFLIB::FiffCoordTransOld * meg_head_t
FWDLIB::FwdCoilSet * eeg_els
FIFFLIB::fiff_int_t * ptrs
#define FIFFV_MNE_COORD_CTF_HEAD
QList< FIFFLIB::FiffChInfo > chs
QSharedPointer< FiffTag > SPtr
Definition: fiff_tag.h:152
#define FIFFV_MNE_COORD_MNI_TAL
MneNamedMatrix * vecs
#define FIFFV_REF_MEG_CH
MNELIB::MneProjOp * proj
DipoleForward ** guess_fwd
Definition: guess_data.h:139
#define FIFFV_COIL_NONE
#define FIFFV_MNE_COORD_FS_TAL_GTZ
DipoleForward description.
easurement data representation in MNE calculations
Definition: mne_meas_data.h:92
fiff_int_t coil_type
Definition: fiff_ch_pos.h:111
Dipole Fit Data class declaration.
#define FIFF_DESCRIPTION
Definition: fiff_file.h:486