MNE-CPP  0.1.9
A Framework for Electrophysiology
compute_fwd.cpp
1 
2 
3 #include "compute_fwd.h"
4 
5 #include <fiff/fiff.h>
6 #include <fiff/c/fiff_coord_trans_old.h>
7 #include "../fwd_coil_set.h"
8 #include "../fwd_comp_data.h"
9 #include <mne/c/mne_ctf_comp_data_set.h>
10 #include "../fwd_eeg_sphere_model_set.h"
11 #include "../fwd_bem_model.h"
12 
13 #include <mne/c/mne_named_matrix.h>
14 #include <mne/c/mne_nearest.h>
15 #include <mne/c/mne_source_space_old.h>
16 #include <mne/mne_forwardsolution.h>
17 
18 #include <fiff/c/fiff_sparse_matrix.h>
19 
20 #include <fiff/fiff_types.h>
21 
22 #include <time.h>
23 
24 #include <Eigen/Dense>
25 
26 #include <QCoreApplication>
27 #include <QFile>
28 #include <QDir>
29 
30 using namespace Eigen;
31 using namespace FWDLIB;
32 using namespace FIFFLIB;
33 using namespace MNELIB;
34 
35 #ifndef TRUE
36 #define TRUE 1
37 #endif
38 
39 #ifndef FALSE
40 #define FALSE 0
41 #endif
42 
43 #ifndef FAIL
44 #define FAIL -1
45 #endif
46 
47 #ifndef OK
48 #define OK 0
49 #endif
50 
51 #define X_41 0
52 #define Y_41 1
53 #define Z_41 2
54 
55 #define ALLOC_ICMATRIX_41(x,y) mne_imatrix_41((x),(y))
56 
57 #define MALLOC_41(x,t) (t *)malloc((x)*sizeof(t))
58 #define REALLOC_41(x,y,t) (t *)((x == Q_NULLPTR) ? malloc((y)*sizeof(t)) : realloc((x),(y)*sizeof(t)))
59 
60 #define FREE_41(x) if ((char *)(x) != Q_NULLPTR) free((char *)(x))
61 
62 #define VEC_COPY_41(to,from) {\
63  (to)[X_41] = (from)[X_41];\
64  (to)[Y_41] = (from)[Y_41];\
65  (to)[Z_41] = (from)[Z_41];\
66  }
67 
68 #define VEC_DOT_41(x,y) ((x)[X_41]*(y)[X_41] + (x)[Y_41]*(y)[Y_41] + (x)[Z_41]*(y)[Z_41])
69 
70 #define VEC_LEN_41(x) sqrt(VEC_DOT_41(x,x))
71 
72 #define ALLOC_CMATRIX_41(x,y) mne_cmatrix_41((x),(y))
73 #define FREE_CMATRIX_41(m) mne_free_cmatrix_41((m))
74 #define FREE_ICMATRIX_41(m) mne_free_icmatrix_41((m))
75 
76 static void matrix_error_41(int kind, int nr, int nc)
77 
78 {
79  if (kind == 1)
80  printf("Failed to allocate memory pointers for a %d x %d matrix\n",nr,nc);
81  else if (kind == 2)
82  printf("Failed to allocate memory for a %d x %d matrix\n",nr,nc);
83  else
84  printf("Allocation error for a %d x %d matrix\n",nr,nc);
85  if (sizeof(void *) == 4) {
86  printf("This is probably because you seem to be using a computer with 32-bit architecture.\n");
87  printf("Please consider moving to a 64-bit platform.");
88  }
89  printf("Cannot continue. Sorry.\n");
90  exit(1);
91 }
92 
93 float **mne_cmatrix_41(int nr,int nc)
94 
95 {
96  int i;
97  float **m;
98  float *whole;
99 
100  m = MALLOC_41(nr,float *);
101  if (!m) matrix_error_41(1,nr,nc);
102  whole = MALLOC_41(nr*nc,float);
103  if (!whole) matrix_error_41(2,nr,nc);
104 
105  for(i=0;i<nr;i++)
106  m[i] = whole + i*nc;
107  return m;
108 }
109 
110 int **mne_imatrix_41(int nr,int nc)
111 
112 {
113  int i,**m;
114  int *whole;
115 
116  m = MALLOC_41(nr,int *);
117  if (!m) matrix_error_41(1,nr,nc);
118  whole = MALLOC_41(nr*nc,int);
119  if (!whole) matrix_error_41(2,nr,nc);
120  for(i=0;i<nr;i++)
121  m[i] = whole + i*nc;
122  return m;
123 }
124 
125 void mne_free_cmatrix_41 (float **m)
126 {
127  if (m) {
128  FREE_41(*m);
129  FREE_41(m);
130  }
131 }
132 
133 void mne_free_icmatrix_41 (int **m)
134 
135 {
136  if (m) {
137  FREE_41(*m);
138  FREE_41(m);
139  }
140 }
141 
142 fiffId get_file_id(const QString& name)
143 {
144  QFile file(name);
145  FiffStream::SPtr stream(new FiffStream(&file));
146  fiffId id;
147  if(!stream->open()) {
148  stream->close();
149  return Q_NULLPTR;
150  }
151  else {
152  id = MALLOC_41(1,fiffIdRec);
153  id->version = stream->id().version;
154  id->machid[0] = stream->id().machid[0];
155  id->machid[1] = stream->id().machid[1];
156  id->time = stream->id().time;
158  stream->close();
159  return id;
160  }
161 }
162 
163 //============================= mne_read_forward_solution.c =============================
164 
165 //int mne_read_meg_comp_eeg_ch_info_41(const QString& name,
166 // QList<FiffChInfo>& megp, /* MEG channels */
167 // int* nmegp,
168 // QList<FiffChInfo>& meg_compp,
169 // int* nmeg_compp,
170 // QList<FiffChInfo>& eegp, /* EEG channels */
171 // int* neegp,
172 // FiffCoordTransOld** meg_head_t,
173 // fiffId* idp) /* The measurement ID */
175 // * Read the channel information and split it into three arrays,
176 // * one for MEG, one for MEG compensation channels, and one for EEG
177 // */
178 //{
179 // QFile file(name);
180 // FiffStream::SPtr stream(new FiffStream(&file));
181 
182 // QList<FiffChInfo> chs;
183 // int nchan = 0;
184 // QList<FiffChInfo> meg;
185 // int nmeg = 0;
186 // QList<FiffChInfo> meg_comp;
187 // int nmeg_comp = 0;
188 // QList<FiffChInfo> eeg;
189 // int neeg = 0;
190 // fiffId id = Q_NULLPTR;
191 // QList<FiffDirNode::SPtr> nodes;
192 // FiffDirNode::SPtr info;
193 // FiffTag::SPtr t_pTag;
194 // FiffChInfo this_ch;
195 // FiffCoordTransOld* t = Q_NULLPTR;
196 // fiff_int_t kind, pos;
197 // int j,k,to_find;
198 
199 // if(!stream->open())
200 // goto bad;
201 
202 // nodes = stream->dirtree()->dir_tree_find(FIFFB_MNE_PARENT_MEAS_FILE);
203 
204 // if (nodes.size() == 0) {
205 // nodes = stream->dirtree()->dir_tree_find(FIFFB_MEAS_INFO);
206 // if (nodes.size() == 0) {
207 // qCritical ("Could not find the channel information.");
208 // goto bad;
209 // }
210 // }
211 // info = nodes[0];
212 // to_find = 0;
213 // for (k = 0; k < info->nent(); k++) {
214 // kind = info->dir[k]->kind;
215 // pos = info->dir[k]->pos;
216 // switch (kind) {
217 // case FIFF_NCHAN :
218 // if (!stream->read_tag(t_pTag,pos))
219 // goto bad;
220 // nchan = *t_pTag->toInt();
221 
222 // for (j = 0; j < nchan; j++) {
223 // chs.append(FiffChInfo());
224 // chs[j].scanNo = -1;
225 // }
226 // to_find = nchan;
227 // break;
228 
229 // case FIFF_PARENT_BLOCK_ID :
230 // if(!stream->read_tag(t_pTag, pos))
231 // goto bad;
232 // // id = t_pTag->toFiffID();
233 // *id = *(fiffId)t_pTag->data();
234 // break;
235 
236 // case FIFF_COORD_TRANS :
237 // if(!stream->read_tag(t_pTag, pos))
238 // goto bad;
239 // // t = t_pTag->toCoordTrans();
240 // t = FiffCoordTransOld::read_helper( t_pTag );
241 // if (t->from != FIFFV_COORD_DEVICE || t->to != FIFFV_COORD_HEAD)
242 // t = Q_NULLPTR;
243 // break;
244 
245 // case FIFF_CH_INFO : /* Information about one channel */
246 // if(!stream->read_tag(t_pTag, pos))
247 // goto bad;
248 // this_ch = t_pTag->toChInfo();
249 // if (this_ch.scanNo <= 0 || this_ch.scanNo > nchan) {
250 // printf ("FIFF_CH_INFO : scan # out of range %d (%d)!",this_ch.scanNo,nchan);
251 // goto bad;
252 // }
253 // else
254 // chs[this_ch.scanNo-1] = this_ch;
255 // to_find--;
256 // break;
257 // }
258 // }
259 // if (to_find != 0) {
260 // qCritical("Some of the channel information was missing.");
261 // goto bad;
262 // }
263 // if (t == Q_NULLPTR && meg_head_t != Q_NULLPTR) {
264 // /*
265 // * Try again in a more general fashion
266 // */
267 // if ((t = FiffCoordTransOld::mne_read_meas_transform(name)) == Q_NULLPTR) {
268 // qCritical("MEG -> head coordinate transformation not found.");
269 // goto bad;
270 // }
271 // }
272 // /*
273 // * Sort out the channels
274 // */
275 // for (k = 0; k < nchan; k++) {
276 // if (chs[k].kind == FIFFV_MEG_CH) {
277 // meg.append(chs[k]);
278 // nmeg++;
279 // } else if (chs[k].kind == FIFFV_REF_MEG_CH) {
280 // meg_comp.append(chs[k]);
281 // nmeg_comp++;
282 // } else if (chs[k].kind == FIFFV_EEG_CH) {
283 // eeg.append(chs[k]);
284 // neeg++;
285 // }
286 // }
287 // // fiff_close(in);
288 
289 // stream->close();
290 
291 // megp = meg;
292 // if(nmegp) {
293 // *nmegp = nmeg;
294 // }
295 
296 // meg_compp = meg_comp;
297 // if(nmeg_compp) {
298 // *nmeg_compp = nmeg_comp;
299 // }
300 
301 // eegp = eeg;
302 // if(neegp) {
303 // *neegp = neeg;
304 // }
305 
306 // if (idp == Q_NULLPTR) {
307 // FREE_41(id);
308 // } else {
309 // *idp = id;
310 // }
311 
312 // if (meg_head_t == Q_NULLPTR) {
313 // FREE_41(t);
314 // } else {
315 // *meg_head_t = t;
316 // }
317 
318 // return FIFF_OK;
319 
320 //bad : {
321 // // fiff_close(in);
322 // stream->close();
323 // FREE_41(id);
324 // // FREE_41(tag.data);
325 // FREE_41(t);
326 // return FIFF_FAIL;
327 // }
328 //}
329 
330 int ComputeFwd::mne_read_meg_comp_eeg_ch_info_41(FIFFLIB::FiffInfoBase::SPtr pFiffInfoBase,
331  QList<FiffChInfo>& listMegCh,
332  int& iNMeg,
333  QList<FiffChInfo>& listMegComp,
334  int& iNMegCmp,
335  QList<FiffChInfo>& listEegCh,
336  int &iNEeg,
337  FiffCoordTransOld** transDevHeadOld,
338  FiffId& id)
339 {
340  int iNumCh = pFiffInfoBase->nchan;
341  for (int k = 0; k < iNumCh; k++) {
342  if (pFiffInfoBase->chs[k].kind == FIFFV_MEG_CH) {
343  listMegCh.append(pFiffInfoBase->chs[k]);
344  iNMeg++;
345  } else if (pFiffInfoBase->chs[k].kind == FIFFV_REF_MEG_CH) {
346  listMegComp.append(pFiffInfoBase->chs[k]);
347  iNMegCmp++;
348  } else if (pFiffInfoBase->chs[k].kind == FIFFV_EEG_CH) {
349  listEegCh.append(pFiffInfoBase->chs[k]);
350  iNEeg++;
351  }
352  }
353 
354  if(!m_pSettings->meg_head_t) {
355  *transDevHeadOld = new FiffCoordTransOld(pFiffInfoBase->dev_head_t.toOld());
356  } else {
357  *transDevHeadOld = m_pSettings->meg_head_t;
358  }
359  if(!m_meg_head_t) {
360  qCritical("MEG -> head coordinate transformation not found.");
361  return FIFF_FAIL;
362  }
363  id = pFiffInfoBase->meas_id;
364  return OK;
365 }
366 
367 int mne_check_chinfo(const QList<FiffChInfo>& chs,
368  int nch)
369 /*
370  * Check that all EEG channels have reasonable locations
371  */
372 {
373  int k;
374  FiffChInfo ch;
375  float close = 0.02f;
376 
377  for (k = 0; k < nch; k++) {
378  if (chs.at(k).kind == FIFFV_EEG_CH) {
379  if (chs.at(k).chpos.r0.norm() < close) {
380  qCritical("Some EEG channels do not have locations assigned.");
381  return FAIL;
382  }
383  }
384  }
385 
386  return OK;
387 }
388 
389 //=============================================================================================================
390 // Temporary Helpers
391 //=============================================================================================================
392 
393 void write_id_old(FiffStream::SPtr& t_pStream, fiff_int_t kind, fiffId id)
394 {
395  fiffId t_id = id;
396  if(t_id->version == -1)
397  {
398  /* initialize random seed: */
399  srand ( time(Q_NULLPTR) );
400  double rand_1 = (double)(rand() % 100);rand_1 /= 100;
401  double rand_2 = (double)(rand() % 100);rand_2 /= 100;
402 
403  time_t seconds;
404  seconds = time (Q_NULLPTR);
405 
406  //fiff_int_t timezone = 5; // Matlab does not know the timezone
407  t_id->version = (1 << 16) | 2; // Version (1 << 16) | 2
408  t_id->machid[0] = 65536*rand_1; // Machine id is random for now
409  t_id->machid[1] = 65536*rand_2; // Machine id is random for now
410  t_id->time.secs = (int)seconds; //seconds since January 1, 1970 //3600*(24*(now-datenum(1970,1,1,0,0,0))+timezone);
411  t_id->time.usecs = 0; // Do not know how we could get this
412  }
413 
414  //
415  //
416  fiff_int_t datasize = 5*4; // The id comprises five integers
417 
418  *t_pStream << (qint32)kind;
419  *t_pStream << (qint32)FIFFT_ID_STRUCT;
420  *t_pStream << (qint32)datasize;
421  *t_pStream << (qint32)FIFFV_NEXT_SEQ;
422  //
423  // Collect the bits together for one write
424  //
425  qint32 data[5];
426  data[0] = t_id->version;
427  data[1] = t_id->machid[0];
428  data[2] = t_id->machid[1];
429  data[3] = t_id->time.secs;
430  data[4] = t_id->time.usecs;
431 
432  for(qint32 i = 0; i < 5; ++i)
433  *t_pStream << data[i];
434 }
435 
436 //=============================================================================================================
437 
438 void write_coord_trans_old(FiffStream::SPtr& t_pStream, const FiffCoordTransOld* trans)
439 {
440  //?typedef struct _fiffCoordTransRec {
441  // fiff_int_t from; /*!< Source coordinate system. */
442  // fiff_int_t to; /*!< Destination coordinate system. */
443  // fiff_float_t rot[3][3]; /*!< The forward transform (rotation part) */
444  // fiff_float_t move[3]; /*!< The forward transform (translation part) */
445  // fiff_float_t invrot[3][3]; /*!< The inverse transform (rotation part) */
446  // fiff_float_t invmove[3]; /*!< The inverse transform (translation part) */
447  //} *fiffCoordTrans, fiffCoordTransRec; /*!< Coordinate transformation descriptor */
448  fiff_int_t datasize = 4*2*12 + 4*2;
449 
450  *t_pStream << (qint32)FIFF_COORD_TRANS;
451  *t_pStream << (qint32)FIFFT_COORD_TRANS_STRUCT;
452  *t_pStream << (qint32)datasize;
453  *t_pStream << (qint32)FIFFV_NEXT_SEQ;
454 
455  //
456  // Start writing fiffCoordTransRec
457  //
458  *t_pStream << (qint32)trans->from;
459  *t_pStream << (qint32)trans->to;
460 
461  //
462  // The transform...
463  //
464  qint32 r, c;
465  for (r = 0; r < 3; ++r)
466  for (c = 0; c < 3; ++c)
467  *t_pStream << (float)trans->rot(r,c);
468  for (r = 0; r < 3; ++r)
469  *t_pStream << (float)trans->move[r];
470 
471  //
472  // ...and its inverse
473  //
474  for (r = 0; r < 3; ++r)
475  for (c = 0; c < 3; ++c)
476  *t_pStream << (float)trans->invrot(r,c);
477  for (r = 0; r < 3; ++r)
478  *t_pStream << (float)trans->invmove[r];
479 }
480 
481 static int **make_file_triangle_list_41(int **tris, int ntri)
482 /*
483  * In the file the numbering starts from one
484  */
485 {
486  int **res = ALLOC_ICMATRIX_41(ntri,3);
487  int j,k;
488 
489  for (j = 0; j < ntri; j++)
490  for (k = 0; k < 3; k++)
491  res[j][k] = tris[j][k]+1;
492  return res;
493 }
494 
495 void mne_write_bad_channel_list_new(FiffStream::SPtr& t_pStream, const QStringList& t_badList)//FILE *out, char **list, int nlist)
496 {
497 
498  t_pStream->start_block(FIFFB_MNE_BAD_CHANNELS);
499  t_pStream->write_name_list(FIFF_MNE_CH_NAME_LIST,t_badList);
500  t_pStream->end_block(FIFFB_MNE_BAD_CHANNELS);
501 
503 
504  // fiff_int_t bad_channel_block = FIFFB_MNE_BAD_CHANNELS;
505  // fiffTagRec bad_channel_block_tags[] = {
506  // { FIFF_BLOCK_START, FIFFT_INT, 0, FIFFV_NEXT_SEQ, Q_NULLPTR },
507  // { FIFF_MNE_CH_NAME_LIST, FIFFT_STRING, 0, FIFFV_NEXT_SEQ, Q_NULLPTR },
508  // { FIFF_BLOCK_END, FIFFT_INT, 0, FIFFV_NEXT_SEQ, Q_NULLPTR }};
509  // int nbad_channel_block_tags = 3;
510  // char *names = Q_NULLPTR;
511  // int k;
512 
513  // if (nlist <= 0)
514  // return OK;
515 
516  // names = mne_name_list_to_string(list,nlist);
517  // bad_channel_block_tags[0].size = sizeof(fiff_int_t);
518  // bad_channel_block_tags[0].data = &bad_channel_block;
519 
520  // bad_channel_block_tags[1].size = strlen(names);
521  // bad_channel_block_tags[1].data = names;
522 
523  // bad_channel_block_tags[2].size = sizeof(fiff_int_t);
524  // bad_channel_block_tags[2].data = &bad_channel_block;
525 
526  // for (k = 0; k < nbad_channel_block_tags; k++)
527  // if (fiff_write_tag(out,bad_channel_block_tags+k) == FIFF_FAIL) {
528  // FREE(names);
529  // return FAIL;
530  // }
531  // FREE(names);
532  // return OK;
533 }
534 
535 void fiff_write_float_matrix_old ( FiffStream::SPtr& t_pStream, /* Destination file name */
536  int kind, /* What kind of tag */
537  fiff_float_t **data, /* The data to write */
538  int rows,
539  int cols) /* Number of rows and columns */
540 /*
541  * Write out a 2D floating-point matrix
542  */
543 {
544  MatrixXf mat(rows,cols);
545 
546  for (int i = 0; i < rows; ++i) {
547  for(int j = 0; j < cols; ++j) {
548  mat(i,j) = data[i][j];
549  }
550  }
551 
552  t_pStream->write_float_matrix(kind, mat);
553 
554  // int res,*dims;
555  // fiffTagRec tag;
556  //#ifdef INTEL_X86_ARCH
557  // int c;
558  //#endif
559  // int k;
560  // int rowsize;
561 
562  // tag.kind = kind;
563  // tag.type = FIFFT_FLOAT | FIFFT_MATRIX;
564  // tag.size = rows*cols*sizeof(fiff_float_t) + 3*sizeof(fiff_int_t);
565  // tag.data = Q_NULLPTR;
566  // tag.next = FIFFV_NEXT_SEQ;
567 
568  // if ((res = fiff_write_tag_info(out,&tag)) == -1)
569  // return FIFF_FAIL;
570 
571  // rowsize = cols*sizeof(fiff_float_t);
572  // for (k = 0; k < rows; k++) {
573  //#ifdef INTEL_X86_ARCH
574  // for (c = 0; c < cols; c++)
575  // swap_float(data[k]+c);
576  //#endif
577  // if (fwrite (data[k],rowsize,1,out) != 1) {
578  // if (ferror(out))
579  // err_set_sys_error("fwrite");
580  // else
581  // err_set_error("fwrite failed");
582  //#ifdef INTEL_X86_ARCH
583  // for (c = 0; c < cols; c++)
584  // swap_float(data[k]+c);
585  //#endif
586  // return FIFF_FAIL;
587  // }
588  //#ifdef INTEL_X86_ARCH
589  // for (c = 0; c < cols; c++)
590  // swap_float(data[k]+c);
591  //#endif
592  // }
593  // dims = MALLOC_41(3,fiff_int_t);
594  // dims[0] = swap_int(cols);
595  // dims[1] = swap_int(rows);
596  // dims[2] = swap_int(2);
597  // if (fwrite (dims,3*sizeof(fiff_int_t),1,out) != 1) {
598  // if (ferror(out))
599  // err_set_sys_error("fwrite");
600  // else
601  // err_set_error("fwrite failed");
602  // FREE(dims);
603  // return FIFF_FAIL;
604  // }
605  // FREE(dims);
606  // return res;
607 }
608 
609 void fiff_write_int_matrix_old ( FiffStream::SPtr& t_pStream,
610  int kind, /* What kind of tag */
611  fiff_int_t **data, /* The data to write */
612  int rows,
613  int cols) /* Number of rows and columns */
614 /*
615  * Write out a 2D integer matrix
616  */
617 {
618  MatrixXi mat(rows,cols);
619 
620  for (int i = 0; i < rows; ++i) {
621  for(int j = 0; j < cols; ++j) {
622  mat(i,j) = data[i][j];
623  }
624  }
625 
626  t_pStream->write_int_matrix(kind, mat);
627 
628  // int res,*dims;
629  // fiffTagRec tag;
630  //#ifdef INTEL_X86_ARCH
631  // int c;
632  //#endif
633  // int k;
634  // int rowsize;
635 
636  // tag.kind = kind;
637  // tag.type = FIFFT_INT | FIFFT_MATRIX;
638  // tag.size = rows*cols*sizeof(fiff_int_t) + 3*sizeof(fiff_int_t);
639  // tag.data = Q_NULLPTR;
640  // tag.next = FIFFV_NEXT_SEQ;
641 
642  // if ((res = fiff_write_tag_info(out,&tag)) == -1)
643  // return -1;
644 
645  // rowsize = cols*sizeof(fiff_int_t);
646  // for (k = 0; k < rows; k++) {
647  //#ifdef INTEL_X86_ARCH
648  // for (c = 0; c < cols; c++)
649  // data[k][c] = swap_int(data[k][c]);
650  //#endif
651  // if (fwrite (data[k],rowsize,1,out) != 1) {
652  // if (ferror(out))
653  // err_set_sys_error("fwrite");
654  // else
655  // err_set_error("fwrite failed");
656  // return -1;
657  //#ifdef INTEL_X86_ARCH
658  // for (c = 0; c < cols; c++)
659  // data[k][c] = swap_int(data[k][c]);
660  //#endif
661  // }
662  //#ifdef INTEL_X86_ARCH
663  // for (c = 0; c < cols; c++)
664  // data[k][c] = swap_int(data[k][c]);
665  //#endif
666  // }
667  // dims = MALLOC_41(3,fiff_int_t);
668  // dims[0] = swap_int(cols);
669  // dims[1] = swap_int(rows);
670  // dims[2] = swap_int(2);
671  // if (fwrite (dims,3*sizeof(fiff_int_t),1,out) != 1) {
672  // if (ferror(out))
673  // err_set_sys_error("fwrite");
674  // else
675  // err_set_error("fwrite failed");
676  // FREE(dims);
677  // return -1;
678  // }
679  // FREE(dims);
680  // return res;
681 }
682 
683 int fiff_write_float_sparse_matrix_old(FiffStream::SPtr& t_pStream, int kind, FiffSparseMatrix* mat)
684 /*
685  * Write a sparse matrix
686  */
687 {
688  FiffTag::SPtr tag;
689  int k;
690  int type;
691  int datasize,idxsize,ptrsize;
692  int two = 2;
693 
694  datasize = mat->nz * sizeof(fiff_float_t);
695  idxsize = mat->nz * sizeof(fiff_int_t);
696  if (mat->coding == FIFFTS_MC_CCS)
697  ptrsize = (mat->n+1) * sizeof(fiff_int_t);
698  else if (mat->coding == FIFFTS_MC_RCS)
699  ptrsize = (mat->m+1) * sizeof(fiff_int_t);
700  else {
701  qCritical("Incomprehensible sparse matrix coding");
702  return FIFF_FAIL;
703  }
704 
705  if (datasize <= 0 || idxsize <= 0 || ptrsize <= 0) {
706  qCritical("fiff_write_float_ccs_matrix: negative vector size(s) in sparse matrix!\n");
707  return FIFF_FAIL;
708  }
709 
710 // tag.kind = kind;
711  if (mat->coding == FIFFTS_MC_CCS)
712  type = FIFFT_FLOAT | FIFFT_CCS_MATRIX;//tag.type = FIFFT_FLOAT | FIFFT_CCS_MATRIX;
713  else if (mat->coding == FIFFTS_MC_RCS)
714  type = FIFFT_FLOAT | FIFFT_RCS_MATRIX;//tag.type = FIFFT_FLOAT | FIFFT_RCS_MATRIX;
715  else {
716  qCritical("Incomprehensible sparse matrix coding");
717  return FIFF_FAIL;
718  }
719 
720 // tag.size = datasize+idxsize+ptrsize+4*sizeof(fiff_int_t);
721 // tag.data = Q_NULLPTR;
722 // tag.next = FIFFV_NEXT_SEQ;
723 
724  //Write Tag Info
725  *t_pStream << (qint32)kind;
726  *t_pStream << (qint32)type;
727  *t_pStream << (qint32)(datasize+idxsize+ptrsize+4*sizeof(fiff_int_t));
728  *t_pStream << (qint32)FIFFV_NEXT_SEQ;
729 // if (fiff_write_tag_info(out,&tag) == FIFF_FAIL)
730 // return FIFF_FAIL;
731 
732  /*
733  * Write data
734  */
735  for(k = 0; k < mat->nz; ++k)
736  *t_pStream << mat->data[k];
737 
738 // /*
739 // * Write data with swapping
740 // */
741 //#ifdef INTEL_X86_ARCH
742 // for (k = 0; k < mat->nz; k++)
743 // swap_floatp(mat->data+k);
744 //#endif
745 // res = fwrite (mat->data,datasize,1,out);
746 //#ifdef INTEL_X86_ARCH
747 // for (k = 0; k < mat->nz; k++)
748 // swap_floatp(mat->data+k);
749 //#endif
750 // if (res != 1)
751 // goto fwrite_error;
752 
753  /*
754  * Write indexes
755  */
756  for(k = 0; k < mat->nz; ++k)
757  *t_pStream << mat->inds[k];
758 
759 // /*
760 // * Write indexes with swapping
761 // */
762 //#ifdef INTEL_X86_ARCH
763 // for (k = 0; k < mat->nz; k++)
764 // swap_intp(mat->inds+k);
765 //#endif
766 // res = fwrite (mat->inds,idxsize,1,out);
767 //#ifdef INTEL_X86_ARCH
768 // for (k = 0; k < mat->nz; k++)
769 // swap_intp(mat->inds+k);
770 //#endif
771 // if (res != 1)
772 // goto fwrite_error;
773 
774  if (mat->coding == FIFFTS_MC_CCS) {
775 
776  for(k = 0; k < mat->n+1; ++k)
777  *t_pStream << mat->ptrs[k];
778 
779 //#ifdef INTEL_X86_ARCH
780 // for (k = 0; k < mat->n+1; k++)
781 // swap_intp(mat->ptrs+k);
782 //#endif
783 // res = fwrite (mat->ptrs,ptrsize,1,out);
784 //#ifdef INTEL_X86_ARCH
785 // for (k = 0; k < mat->n+1; k++)
786 // swap_intp(mat->ptrs+k);
787 //#endif
788 // if (res != 1)
789 // goto fwrite_error;
790  }
791  else { /* Row-compressed format */
792 
793  for(k = 0; k < mat->m+1; ++k)
794  *t_pStream << mat->ptrs[k];
795 
796 //#ifdef INTEL_X86_ARCH
797 // for (k = 0; k < mat->m+1; k++)
798 // swap_intp(mat->ptrs+k);
799 //#endif
800 // res = fwrite (mat->ptrs,ptrsize,1,out);
801 //#ifdef INTEL_X86_ARCH
802 // for (k = 0; k < mat->m+1; k++)
803 // swap_intp(mat->ptrs+k);
804 //#endif
805 // if (res != 1)
806 // goto fwrite_error;
807  }
808  /*
809  * Write the dimensions
810  */
811  *t_pStream << (qint32)mat->nz;
812 // val = swap_int(mat->nz);
813 // if (fwrite (&val,sizeof(fiff_int_t),1,out) != 1)
814 // goto fwrite_error;
815 
816  *t_pStream << (qint32)mat->m;
817 // val = swap_int(mat->m);
818 // if (fwrite (&val,sizeof(fiff_int_t),1,out) != 1)
819 // goto fwrite_error;
820 
821  *t_pStream << (qint32)mat->n;
822 // val = swap_int(mat->n);
823 // if (fwrite (&val,sizeof(fiff_int_t),1,out) != 1)
824 // goto fwrite_error;
825 
826  *t_pStream << (qint32)two;
827 // val = swap_int(two);
828 // if (fwrite (&val,sizeof(fiff_int_t),1,out) != 1)
829 // goto fwrite_error;
830  return FIFF_OK;
831 
832 //fwrite_error : {
837 // return FIFF_FAIL;
838 // }
839 }
840 
841 static int comp_points2(const void *vp1,const void *vp2)
842 
843 {
844  MneNearest* v1 = (MneNearest*)vp1;
845  MneNearest* v2 = (MneNearest*)vp2;
846 
847  if (v1->vert > v2->vert)
848  return 1;
849  else if (v1->vert == v2->vert)
850  return 0;
851  else
852  return -1;
853 }
854 
855 void mne_sort_nearest_by_vertex(MneNearest* points, int npoint)
856 
857 {
858  if (npoint > 1 && points != Q_NULLPTR)
859  qsort(points,npoint,sizeof(MneNearest),comp_points2);
860  return;
861 }
862 
863 FiffSparseMatrix* mne_create_sparse_rcs(int nrow, /* Number of rows */
864  int ncol, /* Number of columns */
865  int *nnz, /* Number of non-zero elements on each row */
866  int **colindex, /* Column indices of non-zero elements on each row */
867  float **vals) /* The nonzero elements on each row
868  * If Q_NULLPTR, the matrix will be all zeroes */
869 
870 {
871  FiffSparseMatrix* sparse = Q_NULLPTR;
872  int j,k,nz,ptr,size,ind;
873  int stor_type = FIFFTS_MC_RCS;
874 
875  for (j = 0, nz = 0; j < nrow; j++)
876  nz = nz + nnz[j];
877 
878  if (nz <= 0) {
879  qCritical("No nonzero elements specified.");
880  return Q_NULLPTR;
881  }
882  size = nz*(sizeof(fiff_float_t) + sizeof(fiff_int_t)) +
883  (nrow+1)*(sizeof(fiff_int_t));
884 
885  sparse = new FiffSparseMatrix;
886  sparse->coding = stor_type;
887  sparse->m = nrow;
888  sparse->n = ncol;
889  sparse->nz = nz;
890  sparse->data = (float *)malloc(size);
891  sparse->inds = (int *)(sparse->data+nz);
892  sparse->ptrs = sparse->inds+nz;
893 
894  for (j = 0, nz = 0; j < nrow; j++) {
895  ptr = -1;
896  for (k = 0; k < nnz[j]; k++) {
897  if (ptr < 0)
898  ptr = nz;
899  ind = sparse->inds[nz] = colindex[j][k];
900  if (ind < 0 || ind >= ncol) {
901  qCritical("Column index out of range in mne_create_sparse_rcs");
902  goto bad;
903  }
904  if (vals)
905  sparse->data[nz] = vals[j][k];
906  else
907  sparse->data[nz] = 0.0;
908  nz++;
909  }
910  sparse->ptrs[j] = ptr;
911  }
912  sparse->ptrs[nrow] = nz;
913  for (j = nrow-1; j >= 0; j--) /* Take care of the empty rows */
914  if (sparse->ptrs[j] < 0)
915  sparse->ptrs[j] = sparse->ptrs[j+1];
916  return sparse;
917 
918 bad : {
919  if(sparse)
920  delete sparse;
921  return Q_NULLPTR;
922  }
923 }
924 
925 FiffSparseMatrix* mne_pick_lower_triangle_rcs(FiffSparseMatrix* mat)
926 /*
927  * Fill in upper triangle with the lower triangle values
928  */
929 {
930  int *nnz = Q_NULLPTR;
931  int **colindex = Q_NULLPTR;
932  float **vals = Q_NULLPTR;
933  FiffSparseMatrix* res = Q_NULLPTR;
934  int i,j,k;
935 
936  if (mat->coding != FIFFTS_MC_RCS) {
937  qCritical("The input matrix to mne_add_upper_triangle_rcs must be in RCS format");
938  goto out;
939  }
940  if (mat->m != mat->n) {
941  qCritical("The input matrix to mne_pick_lower_triangle_rcs must be square");
942  goto out;
943  }
944  /*
945  * Pick the lower triangle elements
946  */
947  nnz = MALLOC_41(mat->m,int);
948  colindex = MALLOC_41(mat->m,int *);
949  vals = MALLOC_41(mat->m,float *);
950  for (i = 0; i < mat->m; i++) {
951  nnz[i] = mat->ptrs[i+1] - mat->ptrs[i];
952  if (nnz[i] > 0) {
953  colindex[i] = MALLOC_41(nnz[i],int);
954  vals[i] = MALLOC_41(nnz[i],float);
955  for (j = mat->ptrs[i], k = 0; j < mat->ptrs[i+1]; j++) {
956  if (mat->inds[j] <= i) {
957  vals[i][k] = mat->data[j];
958  colindex[i][k] = mat->inds[j];
959  k++;
960  }
961  }
962  nnz[i] = k;
963  }
964  else {
965  colindex[i] = Q_NULLPTR;
966  vals[i] = Q_NULLPTR;
967  }
968  }
969  /*
970  * Assemble the matrix
971  */
972  res = mne_create_sparse_rcs(mat->m,mat->n,nnz,colindex,vals);
973 
974 out : {
975  for (i = 0; i < mat->m; i++) {
976  FREE_41(colindex[i]);
977  FREE_41(vals[i]);
978  }
979  FREE_41(nnz);
980  FREE_41(vals);
981  FREE_41(colindex);
982  return res;
983  }
984 }
985 
986 static int write_volume_space_info(FiffStream::SPtr& t_pStream, MneSourceSpaceOld* ss, int selected_only)
987 /*
988  * Write the vertex neighbors and other information for a volume source space
989  */
990 {
991  int ntot,nvert;
992  int *nneighbors = Q_NULLPTR;
993  int *neighbors = Q_NULLPTR;
994  int *inuse_map = Q_NULLPTR;
995  int nneigh,*neigh;
996  int k,p;
997  int res = FAIL;
998 
999  if (ss->type != FIFFV_MNE_SPACE_VOLUME)
1000  return OK;
1001  if (!ss->neighbor_vert || !ss->nneighbor_vert)
1002  return OK;
1003  if (selected_only) {
1004  inuse_map = MALLOC_41(ss->np,int);
1005  for (k = 0,p = 0, ntot = 0; k < ss->np; k++) {
1006  if (ss->inuse[k]) {
1007  ntot += ss->nneighbor_vert[k];
1008  inuse_map[k] = p++;
1009  }
1010  else
1011  inuse_map[k] = -1;
1012  }
1013  nneighbors = MALLOC_41(ss->nuse,int);
1014  neighbors = MALLOC_41(ntot,int);
1015  /*
1016  * Pick the neighbors and fix the vertex numbering to refer
1017  * to the vertices in use only
1018  */
1019  for (k = 0, nvert = 0, ntot = 0; k < ss->np; k++) {
1020  if (ss->inuse[k]) {
1021  neigh = ss->neighbor_vert[k];
1022  nneigh = ss->nneighbor_vert[k];
1023  nneighbors[nvert++] = nneigh;
1024  for (p = 0; p < nneigh; p++)
1025  neighbors[ntot++] = neigh[p] < 0 ? -1 : inuse_map[neigh[p]];
1026  }
1027  }
1028  }
1029  else {
1030  for (k = 0, ntot = 0; k < ss->np; k++)
1031  ntot += ss->nneighbor_vert[k];
1032  nneighbors = MALLOC_41(ss->np,int);
1033  neighbors = MALLOC_41(ntot,int);
1034  nvert = ss->np;
1035  for (k = 0, ntot = 0; k < ss->np; k++) {
1036  neigh = ss->neighbor_vert[k];
1037  nneigh = ss->nneighbor_vert[k];
1038  nneighbors[k] = nneigh;
1039  for (p = 0; p < nneigh; p++)
1040  neighbors[ntot++] = neigh[p];
1041  }
1042  }
1043 
1044  t_pStream->write_int(FIFF_MNE_SOURCE_SPACE_NNEIGHBORS,nneighbors,nvert);
1045  // tag.next = FIFFV_NEXT_SEQ;
1046  // tag.kind = FIFF_MNE_SOURCE_SPACE_NNEIGHBORS;
1047  // tag.type = FIFFT_INT;
1048  // tag.size = nvert*sizeof(fiff_int_t);
1049  // tag.data = (fiff_byte_t *)nneighbors;
1050  // if (fiff_write_tag(out,&tag) == FIFF_FAIL)
1051  // goto out;
1052 
1053  t_pStream->write_int(FIFF_MNE_SOURCE_SPACE_NEIGHBORS,neighbors,ntot);
1054  // tag.next = FIFFV_NEXT_SEQ;
1055  // tag.kind = FIFF_MNE_SOURCE_SPACE_NEIGHBORS;
1056  // tag.type = FIFFT_INT;
1057  // tag.size = ntot*sizeof(fiff_int_t);
1058  // tag.data = (fiff_byte_t *)neighbors;
1059  // if (fiff_write_tag(out,&tag) == FIFF_FAIL)
1060  // goto out;
1061 
1062  /*
1063  * Write some additional stuff
1064  */
1065  if (!selected_only) {
1066  if (ss->voxel_surf_RAS_t) {
1067  write_coord_trans_old(t_pStream, ss->voxel_surf_RAS_t);//t_pStream->write_coord_trans(ss->voxel_surf_RAS_t);
1068 
1069  t_pStream->write_int(FIFF_MNE_SOURCE_SPACE_VOXEL_DIMS,ss->vol_dims,3);
1070  // tag.next = FIFFV_NEXT_SEQ;
1071  // tag.kind = FIFF_MNE_SOURCE_SPACE_VOXEL_DIMS;
1072  // tag.type = FIFFT_INT;
1073  // tag.size = 3*sizeof(fiff_int_t);
1074  // tag.data = (fiff_byte_t *)ss->vol_dims;
1075  // if (fiff_write_tag(out,&tag) == FIFF_FAIL)
1076  // goto out;
1077  }
1078  if (ss->interpolator && !ss->MRI_volume.isEmpty()) {
1079  t_pStream->start_block(FIFFB_MNE_PARENT_MRI_FILE);
1080  if (ss->MRI_surf_RAS_RAS_t)
1081  write_coord_trans_old(t_pStream, ss->MRI_surf_RAS_RAS_t);//t_pStream->write_coord_trans(ss->MRI_surf_RAS_RAS_t);
1082  if (ss->MRI_voxel_surf_RAS_t)
1083  write_coord_trans_old(t_pStream, ss->MRI_voxel_surf_RAS_t);//t_pStream->write_coord_trans(ss->MRI_voxel_surf_RAS_t);
1084  t_pStream->write_string(FIFF_MNE_FILE_NAME,ss->MRI_volume);
1085  if (ss->interpolator)
1086  fiff_write_float_sparse_matrix_old(t_pStream,FIFF_MNE_SOURCE_SPACE_INTERPOLATOR,ss->interpolator);
1087  if (ss->MRI_vol_dims[0] > 0 && ss->MRI_vol_dims[1] > 0 && ss->MRI_vol_dims[2] > 0) {
1088  t_pStream->write_int(FIFF_MRI_WIDTH,&ss->MRI_vol_dims[0]);
1089  t_pStream->write_int(FIFF_MRI_HEIGHT,&ss->MRI_vol_dims[1]);
1090  t_pStream->write_int(FIFF_MRI_DEPTH,&ss->MRI_vol_dims[2]);
1091  }
1092  t_pStream->end_block(FIFFB_MNE_PARENT_MRI_FILE);
1093  }
1094  }
1095  else {
1096  if (ss->interpolator && !ss->MRI_volume.isEmpty()) {
1097  t_pStream->write_string(FIFF_MNE_SOURCE_SPACE_MRI_FILE,ss->MRI_volume);
1098  qCritical("Cannot write the interpolator for selection yet");
1099  goto out;
1100  }
1101  }
1102  res = OK;
1103  goto out;
1104 
1105 out : {
1106  FREE_41(inuse_map);
1107  FREE_41(nneighbors);
1108  FREE_41(neighbors);
1109  return res;
1110  }
1111 }
1112 
1113 int mne_write_one_source_space(FiffStream::SPtr& t_pStream, MneSourceSpaceOld* ss,bool selected_only)
1114 {
1115  float **sel = Q_NULLPTR;
1116  int **tris = Q_NULLPTR;
1117  int *nearest = Q_NULLPTR;
1118  float *nearest_dist = Q_NULLPTR;
1119  int p,pp;
1120 
1121  if (ss->np <= 0) {
1122  qCritical("No points in the source space being saved");
1123  goto bad;
1124  }
1125 
1126  t_pStream->start_block(FIFFB_MNE_SOURCE_SPACE);
1127 
1128  /*
1129  * General information
1130  */
1131  if (ss->type != FIFFV_MNE_SPACE_UNKNOWN)
1132  t_pStream->write_int(FIFF_MNE_SOURCE_SPACE_TYPE,&ss->type);
1133  if (ss->id != FIFFV_MNE_SURF_UNKNOWN)
1134  t_pStream->write_int(FIFF_MNE_SOURCE_SPACE_ID,&ss->id);
1135  if (!ss->subject.isEmpty() && ss->subject.size() > 0) {
1136  QString subj(ss->subject);
1137  t_pStream->write_string(FIFF_SUBJ_HIS_ID,subj);
1138  }
1139 
1140  t_pStream->write_int(FIFF_MNE_COORD_FRAME,&ss->coord_frame);
1141 
1142  if (selected_only) {
1143  if (ss->nuse == 0) {
1144  qCritical("No vertices in use. Cannot write active-only vertices from this source space");
1145  goto bad;
1146  }
1147  sel = ALLOC_CMATRIX_41(ss->nuse,3);
1148  t_pStream->write_int(FIFF_MNE_SOURCE_SPACE_NPOINTS,&ss->nuse);
1149  for (p = 0, pp = 0; p < ss->np; p++) {
1150  if (ss->inuse[p]) {
1151  sel[pp][X_41] = ss->rr[p][X_41];
1152  sel[pp][Y_41] = ss->rr[p][Y_41];
1153  sel[pp][Z_41] = ss->rr[p][Z_41];
1154  pp++;
1155  }
1156  }
1157  fiff_write_float_matrix_old(t_pStream, FIFF_MNE_SOURCE_SPACE_POINTS, sel,ss->nuse,3);
1158 
1159  for (p = 0, pp = 0; p < ss->np; p++) {
1160  if (ss->inuse[p]) {
1161  sel[pp][X_41] = ss->nn[p][X_41];
1162  sel[pp][Y_41] = ss->nn[p][Y_41];
1163  sel[pp][Z_41] = ss->nn[p][Z_41];
1164  pp++;
1165  }
1166  }
1167  fiff_write_float_matrix_old(t_pStream, FIFF_MNE_SOURCE_SPACE_NORMALS, sel,ss->nuse,3);
1168  FREE_CMATRIX_41(sel); sel = Q_NULLPTR;
1169 #ifdef WRONG
1170  /*
1171  * This code is incorrect because the numbering in the nuse triangulation refers to the complete source space
1172  */
1173  if (ss->nuse_tri > 0) { /* Write the triangulation information */
1174  /*
1175  * The 'use' triangulation is identical to the complete one
1176  */
1177  if (fiff_write_int_tag(out,FIFF_MNE_SOURCE_SPACE_NTRI,ss->nuse_tri) == FIFF_FAIL)
1178  goto bad;
1179  tris = make_file_triangle_list(ss->use_itris,ss->nuse_tri);
1180  if (fiff_write_int_matrix(out,FIFF_MNE_SOURCE_SPACE_TRIANGLES,tris,
1181  ss->nuse_tri,3) == FIFF_FAIL)
1182  goto bad;
1183 
1184  if (fiff_write_int_tag(out,FIFF_MNE_SOURCE_SPACE_NUSE_TRI,ss->nuse_tri) == FIFF_FAIL)
1185  goto bad;
1186  if (fiff_write_int_matrix(out,FIFF_MNE_SOURCE_SPACE_USE_TRIANGLES,tris,
1187  ss->nuse_tri,3) == FIFF_FAIL)
1188  goto bad;
1189  FREE_ICMATRIX(tris); tris = Q_NULLPTR;
1190  }
1191 #endif
1192  }
1193  else {
1194  // fiffTagRec tag;
1195  t_pStream->write_int(FIFF_MNE_SOURCE_SPACE_NPOINTS,&ss->np);
1196 
1197  fiff_write_float_matrix_old(t_pStream, FIFF_MNE_SOURCE_SPACE_POINTS, ss->rr, ss->np, 3);
1198 
1199  fiff_write_float_matrix_old(t_pStream, FIFF_MNE_SOURCE_SPACE_NORMALS, ss->nn, ss->np, 3);
1200 
1201  if (ss->nuse > 0 && ss->inuse) {
1202  t_pStream->write_int(FIFF_MNE_SOURCE_SPACE_SELECTION,ss->inuse,ss->np);
1203  // tag.next = 0;
1204  // tag.kind = FIFF_MNE_SOURCE_SPACE_SELECTION;
1205  // tag.type = FIFFT_INT;
1206  // tag.size = (ss->np)*sizeof(fiff_int_t);
1207  // tag.data = (fiff_byte_t *)(ss->inuse);
1208  // if (fiff_write_tag(out,&tag) == FIFF_FAIL)
1209  // goto bad;
1210  t_pStream->write_int(FIFF_MNE_SOURCE_SPACE_NUSE,&ss->nuse);
1211  // if (fiff_write_int_tag (out, FIFF_MNE_SOURCE_SPACE_NUSE,ss->nuse) == FIFF_FAIL)
1212  // goto bad;
1213  }
1214  if (ss->ntri > 0) { /* Write the triangulation information */
1215  t_pStream->write_int(FIFF_MNE_SOURCE_SPACE_NTRI,&ss->ntri);
1216  // if (fiff_write_int_tag(out,FIFF_MNE_SOURCE_SPACE_NTRI,ss->ntri) == FIFF_FAIL)
1217  // goto bad;
1218  tris = make_file_triangle_list_41(ss->itris,ss->ntri);
1219 
1220  fiff_write_int_matrix_old(t_pStream, FIFF_MNE_SOURCE_SPACE_TRIANGLES, tris, ss->ntri, 3);
1221  // if (fiff_write_int_matrix(out,FIFF_MNE_SOURCE_SPACE_TRIANGLES,tris,
1222  // ss->ntri,3) == FIFF_FAIL)
1223  // goto bad;
1224  FREE_ICMATRIX_41(tris); tris = Q_NULLPTR;
1225  }
1226  if (ss->nuse_tri > 0) { /* Write the triangulation information for the vertices in use */
1227  t_pStream->write_int(FIFF_MNE_SOURCE_SPACE_NUSE_TRI,&ss->nuse_tri);
1228  // if (fiff_write_int_tag(out,FIFF_MNE_SOURCE_SPACE_NUSE_TRI,ss->nuse_tri) == FIFF_FAIL)
1229  // goto bad;
1230  tris = make_file_triangle_list_41(ss->use_itris,ss->nuse_tri);
1231 
1232  fiff_write_int_matrix_old(t_pStream, FIFF_MNE_SOURCE_SPACE_USE_TRIANGLES, tris, ss->nuse_tri, 3);
1233  // if (fiff_write_int_matrix(out,FIFF_MNE_SOURCE_SPACE_USE_TRIANGLES,tris,
1234  // ss->nuse_tri,3) == FIFF_FAIL)
1235  // goto bad;
1236  FREE_ICMATRIX_41(tris); tris = Q_NULLPTR;
1237  }
1238  if (ss->nearest) { /* Write the patch information */
1239  nearest = MALLOC_41(ss->np,int);
1240  nearest_dist = MALLOC_41(ss->np,float);
1241 
1242  mne_sort_nearest_by_vertex(ss->nearest,ss->np);
1243  for (p = 0; p < ss->np; p++) {
1244  nearest[p] = ss->nearest[p].nearest;
1245  nearest_dist[p] = ss->nearest[p].dist;
1246  }
1247 
1248  t_pStream->write_int(FIFF_MNE_SOURCE_SPACE_NEAREST,nearest,ss->np);
1249  // tag.next = FIFFV_NEXT_SEQ;
1250  // tag.kind = FIFF_MNE_SOURCE_SPACE_NEAREST;
1251  // tag.type = FIFFT_INT;
1252  // tag.size = (ss->np)*sizeof(fiff_int_t);
1253  // tag.data = (fiff_byte_t *)(nearest);
1254  // if (fiff_write_tag(out,&tag) == FIFF_FAIL)
1255  // goto bad;
1256 
1257  t_pStream->write_float(FIFF_MNE_SOURCE_SPACE_NEAREST_DIST,nearest_dist,ss->np);
1258  // tag.next = FIFFV_NEXT_SEQ;
1259  // tag.kind = FIFF_MNE_SOURCE_SPACE_NEAREST_DIST;
1260  // tag.type = FIFFT_FLOAT;
1261  // tag.size = (ss->np)*sizeof(fiff_float_t);
1262  // tag.data = (fiff_byte_t *)(nearest_dist);
1263  // if (fiff_write_tag(out,&tag) == FIFF_FAIL)
1264  // goto bad;
1265 
1266  FREE_41(nearest); nearest = Q_NULLPTR;
1267  FREE_41(nearest_dist); nearest_dist = Q_NULLPTR;
1268  }
1269  if (ss->dist) { /* Distance information */
1270  FiffSparseMatrix* m = mne_pick_lower_triangle_rcs(ss->dist);
1271  if (!m)
1272  goto bad;
1273  if (fiff_write_float_sparse_matrix_old(t_pStream,FIFF_MNE_SOURCE_SPACE_DIST,m) == FIFF_FAIL) {
1274  if(m)
1275  delete m;
1276  goto bad;
1277  }
1278  if(m)
1279  delete m;
1280 
1281  t_pStream->write_float(FIFF_MNE_SOURCE_SPACE_DIST_LIMIT,&ss->dist_limit);
1282  }
1283  }
1284  /*
1285  * Volume source spaces have additional information
1286  */
1287  // if (write_volume_space_info(out,ss,selected_only) == FIFF_FAIL)
1288  // goto bad;
1289 
1290  t_pStream->end_block(FIFFB_MNE_SOURCE_SPACE);
1291  return FIFF_OK;
1292 
1293 bad : {
1294  FREE_ICMATRIX_41(tris);
1295  FREE_CMATRIX_41(sel);
1296  FREE_41(nearest);
1297  FREE_41(nearest_dist);
1298  return FIFF_FAIL;
1299  }
1300 }
1301 
1302 QString mne_name_list_to_string_41(const QStringList& list)
1303 /*
1304  * Convert a string array to a colon-separated string
1305  */
1306 {
1307  int nlist = list.size();
1308  QString res;
1309  if (nlist == 0 || list.isEmpty())
1310  return res;
1311 // res[0] = '\0';
1312  for (int k = 0; k < nlist-1; k++) {
1313  res += list[k];
1314  res += ":";
1315  }
1316  res += list[nlist-1];
1317  return res;
1318 }
1319 
1320 int mne_write_named_matrix( FiffStream::SPtr& t_pStream,
1321  int kind,
1322  MneNamedMatrix* mat)
1323 /*
1324  * Write a block which contains information about one named matrix
1325  */
1326 {
1327  QString names;
1328 
1329  t_pStream->start_block(FIFFB_MNE_NAMED_MATRIX);
1330 
1331  t_pStream->write_int(FIFF_MNE_NROW,&mat->nrow);
1332  t_pStream->write_int(FIFF_MNE_NCOL,&mat->ncol);
1333  if (!mat->rowlist.isEmpty()) {
1334  names = mne_name_list_to_string_41(mat->rowlist);
1335  t_pStream->write_string(FIFF_MNE_ROW_NAMES,names);
1336  }
1337  if (!mat->collist.isEmpty()) {
1338  names = mne_name_list_to_string_41(mat->collist);
1339  t_pStream->write_string(FIFF_MNE_COL_NAMES,names);
1340  }
1341  fiff_write_float_matrix_old (t_pStream,kind,mat->data,mat->nrow,mat->ncol);
1342 
1343  t_pStream->end_block(FIFFB_MNE_NAMED_MATRIX);
1344 
1345  return FIFF_OK;
1346 }
1347 
1348 bool fiff_put_dir (FiffStream::SPtr& t_pStream, const QList<FiffDirEntry::SPtr>& dir)
1349 /*
1350  * Put in new directory
1351  */
1352 {
1353  int nent = dir.size();
1354  int k;
1355  FiffTag::SPtr t_pTag;
1356  fiff_int_t dirpos;
1357 
1358  for (k = 0; k < nent; k++) {
1359  if (dir[k]->kind == FIFF_DIR_POINTER) {
1360  /*
1361  * Read current value of directory pointer
1362  */
1363  if (!t_pStream->read_tag(t_pTag,dir[k]->pos)) {
1364  fprintf (stderr,"Could not read FIFF_DIR_POINTER!\n");
1365  return false;
1366  }
1367  /*
1368  * If there is no directory, append the new one
1369  */
1370  dirpos = *t_pTag->toInt();//2GB restriction
1371  if (dirpos <= 0)
1372  dirpos = -1;
1373  /*
1374  * Put together the directory tag
1375  */
1376  dirpos = (fiff_int_t)t_pStream->write_dir_entries(dir);//2GB restriction
1377  if (dirpos < 0)
1378  printf ("Could not update directory!\n");
1379  else {
1380  t_pTag->setNum(dirpos);
1381 // t_pStream->write_tag(t_pTag,dir[k]->pos);
1382 
1383  t_pStream->write_dir_pointer(dirpos, dir[k]->pos);
1384 
1385 // t_pStream->device()->seek(dir[k]->pos);
1386 
1387 // fiff_int_t datasize = 1 * 4;
1388 
1389 // *t_pStream << (qint32)t_pTag->kind;
1390 // *t_pStream << (qint32)t_pTag->type;
1391 // *t_pStream << (qint32)datasize;
1392 // *t_pStream << (qint32)t_pTag->next;
1393 
1394 // *t_pStream << dirpos;
1395 
1396  }
1397  return true;
1398  }
1399  }
1400  printf ("Could not find place for directory!\n");
1401  return false;
1402 }
1403 
1404 //============================= write_solution.c =============================
1405 
1406 bool write_solution(const QString& name, /* Destination file */
1407  MneSourceSpaceOld* *spaces, /* The source spaces */
1408  int nspace,
1409  const QString& mri_file, /* MRI file and data obtained from there */
1410  fiffId mri_id,
1411  FiffCoordTransOld* mri_head_t,
1412  const QString& meas_file, /* MEG file and data obtained from there */
1413  FiffId meas_id,
1414  FiffCoordTransOld* meg_head_t,
1415  QList<FiffChInfo> meg_chs,
1416  int nmeg,
1417  QList<FiffChInfo> eeg_chs,
1418  int neeg,
1419  int fixed_ori, /* Fixed orientation dipoles? */
1420  int coord_frame, /* Coordinate frame employed in the forward calculations */
1421  FiffNamedMatrix& meg_solution,
1422  FiffNamedMatrix& eeg_solution,
1423  FiffNamedMatrix& meg_solution_grad,
1424  FiffNamedMatrix& eeg_solution_grad,
1425  bool bDoGrad)
1426 {
1427  // New Stuff
1428  QFile file(name);
1429 
1430  QFile fileIn(name);
1431  FiffStream::SPtr t_pStreamIn;
1432 
1433  int nvert;
1434  int k;
1435 
1436  //
1437  // Open the file, create directory
1438  //
1439 
1440  // Create the file and save the essentials
1441  FiffStream::SPtr t_pStream = FiffStream::start_file(file);
1442 
1443  t_pStream->start_block(FIFFB_MNE);
1444 
1445  /*
1446  * Information from the MRI file
1447  */
1448  {
1449  t_pStream->start_block(FIFFB_MNE_PARENT_MRI_FILE);
1450 
1451  t_pStream->write_string(FIFF_MNE_FILE_NAME, mri_file);
1452  if (mri_id != Q_NULLPTR)
1453  write_id_old(t_pStream, FIFF_PARENT_FILE_ID, mri_id);//t_pStream->write_id(FIFF_PARENT_FILE_ID, mri_id);
1454  write_coord_trans_old(t_pStream, mri_head_t);//t_pStream->write_coord_trans(mri_head_t);
1455 
1456  t_pStream->end_block(FIFFB_MNE_PARENT_MRI_FILE);
1457  }
1458 
1459  /*
1460  * Information from the MEG file
1461  */
1462  {
1463  QStringList file_bads;
1464 
1465  t_pStream->start_block(FIFFB_MNE_PARENT_MEAS_FILE);
1466 
1467  t_pStream->write_string(FIFF_MNE_FILE_NAME, meas_file);
1468  if (!meas_id.isEmpty()) {
1469  t_pStream->write_id(FIFF_PARENT_BLOCK_ID, meas_id);
1470  }
1471  write_coord_trans_old(t_pStream, meg_head_t);//t_pStream->write_coord_trans(meg_head_t);
1472 
1473  int nchan = nmeg+neeg;
1474  t_pStream->write_int(FIFF_NCHAN,&nchan);
1475 
1476  FiffChInfo chInfo;
1477  int k, p;
1478  for (k = 0, p = 0; k < nmeg; k++) {
1479 // meg_chs[k].scanNo = ++p;
1480 // chInfo.scanNo = meg_chs[k].scanNo;
1481 // chInfo.logNo = meg_chs[k].logNo;
1482 // chInfo.kind = meg_chs[k].kind;
1483 // chInfo.range = meg_chs[k].range;
1484 // chInfo.cal = meg_chs[k].cal;
1485 // chInfo.chpos.coil_type = meg_chs[k].chpos.coil_type;
1486 // chInfo.chpos.r0[0] = meg_chs[k].chpos.r0[0];
1487 // chInfo.chpos.r0[1] = meg_chs[k].chpos.r0[1];
1488 // chInfo.chpos.r0[2] = meg_chs[k].chpos.r0[2];
1489 // chInfo.chpos.ex[0] = meg_chs[k].chpos.ex[0];
1490 // chInfo.chpos.ex[1] = meg_chs[k].chpos.ex[1];
1491 // chInfo.chpos.ex[2] = meg_chs[k].chpos.ex[2];
1492 // chInfo.chpos.ey[0] = meg_chs[k].chpos.ey[0];
1493 // chInfo.chpos.ey[1] = meg_chs[k].chpos.ey[1];
1494 // chInfo.chpos.ey[2] = meg_chs[k].chpos.ey[2];
1495 // chInfo.chpos.ez[0] = meg_chs[k].chpos.ez[0];
1496 // chInfo.chpos.ez[1] = meg_chs[k].chpos.ez[1];
1497 // chInfo.chpos.ez[2] = meg_chs[k].chpos.ez[2];
1498 // chInfo.unit = meg_chs[k].unit;
1499 // chInfo.unit_mul = meg_chs[k].unit_mul;
1500 // chInfo.ch_name = QString(meg_chs[k].ch_name);
1501 
1502  t_pStream->write_ch_info(meg_chs[k]);
1503  }
1504 
1505  for (k = 0; k < neeg; k++) {
1506 // eeg_chs[k].scanNo = ++p;
1507 // chInfo.scanNo = eeg_chs[k].scanNo;
1508 // chInfo.logNo = eeg_chs[k].logNo;
1509 // chInfo.kind = eeg_chs[k].kind;
1510 // chInfo.range = eeg_chs[k].range;
1511 // chInfo.cal = eeg_chs[k].cal;
1512 // chInfo.chpos.coil_type = eeg_chs[k].chpos.coil_type;
1513 // chInfo.chpos.r0[0] = eeg_chs[k].chpos.r0[0];
1514 // chInfo.chpos.r0[1] = eeg_chs[k].chpos.r0[1];
1515 // chInfo.chpos.r0[2] = eeg_chs[k].chpos.r0[2];
1516 // chInfo.chpos.ex[0] = eeg_chs[k].chpos.ex[0];
1517 // chInfo.chpos.ex[1] = eeg_chs[k].chpos.ex[1];
1518 // chInfo.chpos.ex[2] = eeg_chs[k].chpos.ex[2];
1519 // chInfo.chpos.ey[0] = eeg_chs[k].chpos.ey[0];
1520 // chInfo.chpos.ey[1] = eeg_chs[k].chpos.ey[1];
1521 // chInfo.chpos.ey[2] = eeg_chs[k].chpos.ey[2];
1522 // chInfo.chpos.ez[0] = eeg_chs[k].chpos.ez[0];
1523 // chInfo.chpos.ez[1] = eeg_chs[k].chpos.ez[1];
1524 // chInfo.chpos.ez[2] = eeg_chs[k].chpos.ez[2];
1525 // chInfo.unit = eeg_chs[k].unit;
1526 // chInfo.unit_mul = eeg_chs[k].unit_mul;
1527 // chInfo.ch_name = QString(eeg_chs[k].ch_name);
1528 
1529  t_pStream->write_ch_info(eeg_chs[k]);
1530  }
1531  /*
1532  * Copy the bad channel list from the measurement file
1533  */
1534 
1535  //
1536  // mne_read_bad_channel_list replacement
1537  //
1538  QFile fileBad(meas_file);
1539  FiffStream::SPtr t_pStreamBads(new FiffStream(&fileBad));
1540  if(!t_pStreamBads->open())
1541  return false;
1542  file_bads = t_pStreamBads->read_bad_channels(t_pStreamBads->dirtree());
1543 
1544  //
1545  // mne_write_bad_channel_list replacement
1546  //
1547  mne_write_bad_channel_list_new(t_pStream,file_bads);
1548 
1549  t_pStream->end_block(FIFFB_MNE_PARENT_MEAS_FILE);
1550  }
1551 
1552  /*
1553  * Write the source spaces (again)
1554  */
1555  for (k = 0, nvert = 0; k < nspace; k++) {
1556  if (mne_write_one_source_space(t_pStream,spaces[k],FALSE) == FIFF_FAIL)
1557  goto bad;
1558  nvert += spaces[k]->nuse;
1559  }
1560 
1561  /*
1562  * MEG forward solution
1563  */
1564  if (nmeg > 0) {
1565  t_pStream->start_block(FIFFB_MNE_FORWARD_SOLUTION);
1566 
1567  int val = FIFFV_MNE_MEG;
1568  t_pStream->write_int(FIFF_MNE_INCLUDED_METHODS,&val);
1569  t_pStream->write_int(FIFF_MNE_COORD_FRAME,&coord_frame);
1570  val = fixed_ori ? FIFFV_MNE_FIXED_ORI : FIFFV_MNE_FREE_ORI;
1571  t_pStream->write_int(FIFF_MNE_SOURCE_ORIENTATION,&val);
1572  t_pStream->write_int(FIFF_MNE_SOURCE_SPACE_NPOINTS,&nvert);
1573  t_pStream->write_int(FIFF_NCHAN,&nmeg);
1574 
1575  meg_solution.transpose_named_matrix();
1576  t_pStream->write_named_matrix(FIFF_MNE_FORWARD_SOLUTION,meg_solution);
1577  meg_solution.transpose_named_matrix();
1578 // if (mne_write_named_matrix(t_pStream,FIFF_MNE_FORWARD_SOLUTION,meg_solution) == FIFF_FAIL)
1579 // goto bad;
1580  if (bDoGrad) {
1581  meg_solution_grad.transpose_named_matrix();
1582  t_pStream->write_named_matrix(FIFF_MNE_FORWARD_SOLUTION_GRAD,meg_solution_grad);
1583  meg_solution_grad.transpose_named_matrix();
1584  }
1585 
1586 // if (mne_write_named_matrix(t_pStream,FIFF_MNE_FORWARD_SOLUTION_GRAD,meg_solution_grad) == FIFF_FAIL)
1587 // goto bad;
1588  t_pStream->end_block(FIFFB_MNE_FORWARD_SOLUTION);
1589  }
1590  /*
1591  * EEG forward solution
1592  */
1593  if (neeg > 0) {
1594  t_pStream->start_block(FIFFB_MNE_FORWARD_SOLUTION);
1595 
1596  int val = FIFFV_MNE_EEG;
1597  t_pStream->write_int(FIFF_MNE_INCLUDED_METHODS,&val);
1598  t_pStream->write_int(FIFF_MNE_COORD_FRAME,&coord_frame);
1599  val = fixed_ori ? FIFFV_MNE_FIXED_ORI : FIFFV_MNE_FREE_ORI;
1600  t_pStream->write_int(FIFF_MNE_SOURCE_ORIENTATION,&val);
1601  t_pStream->write_int(FIFF_NCHAN,&neeg);
1602  t_pStream->write_int(FIFF_MNE_SOURCE_SPACE_NPOINTS,&nvert);
1603 
1604  eeg_solution.transpose_named_matrix();
1605  t_pStream->write_named_matrix(FIFF_MNE_FORWARD_SOLUTION,eeg_solution);
1606 // if (mne_write_named_matrix(t_pStream,FIFF_MNE_FORWARD_SOLUTION,eeg_solution) == FIFF_FAIL)
1607  eeg_solution.transpose_named_matrix();
1608 // goto bad;
1609  if (bDoGrad) {
1610  eeg_solution_grad.transpose_named_matrix();
1611  t_pStream->write_named_matrix(FIFF_MNE_FORWARD_SOLUTION_GRAD,eeg_solution_grad);
1612  eeg_solution_grad.transpose_named_matrix();
1613  }
1614 // if (mne_write_named_matrix(t_pStream,FIFF_MNE_FORWARD_SOLUTION_GRAD,eeg_solution_grad) == FIFF_FAIL)
1615 // goto bad;
1616 
1617  t_pStream->end_block(FIFFB_MNE_FORWARD_SOLUTION);
1618  }
1619 
1620  t_pStream->end_block(FIFFB_MNE);
1621  t_pStream->end_file();
1622  t_pStream->close();
1623  t_pStream.clear();
1624 
1625  /*
1626  * Add directory
1627  */
1628  t_pStreamIn = FiffStream::open_update(fileIn);
1629 
1630  if (!fiff_put_dir(t_pStreamIn,t_pStreamIn->dir()))
1631  goto bad;
1632  if(t_pStreamIn)
1633  t_pStreamIn->close();
1634 
1635  return true;
1636 
1637 bad : {
1638  if(t_pStream)
1639  t_pStream->close();
1640  if(t_pStreamIn)
1641  t_pStreamIn->close();
1642  return false;
1643  }
1644 }
1645 
1646 /*
1647  * Process the environment information
1648  */
1649 
1650 bool mne_attach_env(const QString& name, const QString& command)
1651 /*
1652  * Add the environment info for future reference
1653  */
1654 {
1655  int insert_blocks[] = { FIFFB_MNE , FIFFB_MEAS, FIFFB_MRI, FIFFB_BEM, -1 };
1656  QString cwd = QDir::currentPath();
1657 
1658  int b,k, insert;
1659  FiffTag::SPtr t_pTag;
1660 // QList<FiffTag::SPtr> tags;
1661  QFile fileInOut(name);
1662  FiffStream::SPtr t_pStreamInOut;
1663 
1664 // if (fiff_new_file_id(&id) == FIFF_FAIL)
1665 // return false;
1666  FiffId id(FiffId::new_file_id());
1667 
1668 //#ifdef DEBUG
1669 // printf("\n");
1670 // printf("cwd = %s\n",cwd);
1671 // printf("com = %s\n",command);
1672 // printf("envid = %s\n",mne_format_file_id(&id));
1673 //#endif
1674 
1675  if (!fileInOut.exists()) {
1676  qCritical("File %s does not exist. Cannot attach env info.",name.toUtf8().constData());
1677  return false;
1678  }
1679 // if (!fileInOut.isWritable()) {
1680 // qCritical("File %s is not writable. Cannot attach env info.",name.toUtf8().constData());
1681 // return false;
1682 // }
1683  /*
1684  * Open the file to modify
1685  */
1686  if (!(t_pStreamInOut = FiffStream::open_update(fileInOut)))
1687  return false;
1688  /*
1689  * Find an appropriate position to insert
1690  */
1691  for (insert = -1, b = 0; insert_blocks[b] >= 0; b++) {
1692  for (k = 0; k < t_pStreamInOut->nent(); k++) {
1693  if (t_pStreamInOut->dir()[k]->kind == FIFF_BLOCK_START) {
1694  if (!t_pStreamInOut->read_tag(t_pTag, t_pStreamInOut->dir()[k]->pos))
1695  return false;
1696  if (*(t_pTag->toInt()) == insert_blocks[b]) {
1697  insert = k;
1698  break;
1699  }
1700  }
1701  }
1702  if (insert >= 0)
1703  break;
1704  }
1705  if (insert < 0) {
1706  qCritical("Suitable place for environment insertion not found.");
1707  return false;
1708  }
1709 
1710  /*
1711  * Do not build the list of tags to insert -> Do insertion right away
1712  */
1713 
1714  // Modified of fiff_insert_after
1715  int where = insert;
1716  /*
1717  * Insert new tags into a file
1718  * The directory information in dest is updated
1719  */
1720  if (where < 0 || where >= t_pStreamInOut->nent()-1) {
1721  qCritical("illegal insertion point in fiff_insert_after!");
1722  return false;
1723  }
1724 
1725  FiffTag::SPtr t_pTagNext;
1726  QList<FiffDirEntry::SPtr> old_dir = t_pStreamInOut->dir();
1727  QList<FiffDirEntry::SPtr> this_ent = old_dir.mid(where);//this_ent = old_dir + where;
1728 
1729  if (!t_pStreamInOut->read_tag(t_pTagNext, this_ent[0]->pos))
1730  return false;
1731  /*
1732  * Update next info to be sequential
1733  */
1734  qint64 next_tmp = t_pStreamInOut->device()->pos();
1735  /*
1736  * Go to the end of the file
1737  */
1738  t_pStreamInOut->device()->seek(fileInOut.size());//SEEK_END
1739  /*
1740  * Allocate new directory
1741  * Copy the beginning of old directory
1742  */
1743  QList<FiffDirEntry::SPtr> new_dir = old_dir.mid(0,where+1);
1744 
1745  /*
1746  * Save the old size for future purposes
1747  */
1748  qint64 old_end = t_pStreamInOut->device()->pos();
1749  /*
1750  * Write tags, check for errors
1751  */
1752 
1753  //Don't use the for loop here instead do it explicitly for specific tags
1754  FiffDirEntry::SPtr new_this;
1755 
1756  new_this = FiffDirEntry::SPtr(new FiffDirEntry);
1757  new_this->kind = FIFF_BLOCK_START;
1758  new_this->type = FIFFT_INT;
1759  new_this->size = 1 * 4;
1760  new_this->pos = t_pStreamInOut->start_block(FIFFB_MNE_ENV);
1761  new_dir.append(new_this);
1762 
1763  new_this = FiffDirEntry::SPtr(new FiffDirEntry);
1764  new_this->kind = FIFF_BLOCK_ID;
1765  new_this->type = FIFFT_ID_STRUCT;
1766  new_this->size = 5 * 4;
1767  new_this->pos = t_pStreamInOut->write_id(FIFF_BLOCK_ID,id);
1768  new_dir.append(new_this);
1769 
1770  new_this = FiffDirEntry::SPtr(new FiffDirEntry);
1771  new_this->kind = FIFF_MNE_ENV_WORKING_DIR;
1772  new_this->type = FIFFT_STRING;
1773  new_this->size = cwd.size();
1774  new_this->pos = t_pStreamInOut->write_string(FIFF_MNE_ENV_WORKING_DIR,cwd);
1775  new_dir.append(new_this);
1776 
1777  new_this = FiffDirEntry::SPtr(new FiffDirEntry);
1778  new_this->kind = FIFF_MNE_ENV_COMMAND_LINE;
1779  new_this->type = FIFFT_STRING;
1780  new_this->size = command.size();
1781  new_this->pos = t_pStreamInOut->write_string(FIFF_MNE_ENV_COMMAND_LINE,command);
1782  new_dir.append(new_this);
1783 
1784  new_this = FiffDirEntry::SPtr(new FiffDirEntry);
1785  new_this->kind = FIFF_BLOCK_END;
1786  new_this->type = FIFFT_INT;
1787  new_this->size = 1 * 4;
1788 
1789  new_this->pos = t_pStreamInOut->end_block(FIFFB_MNE_ENV,next_tmp);
1790  new_dir.append(new_this);
1791 
1792  /*
1793  * Copy the rest of the old directory
1794  */
1795  new_dir.append(old_dir.mid(where+1));
1796  /*
1797  * Now, it is time to update the braching tag
1798  * If something goes wrong here, we cannot be sure that
1799  * the file is readable. Let's hope for the best...
1800  */
1801  t_pTagNext->next = (qint32)old_end;//2GB cut of
1802  t_pStreamInOut->write_tag(t_pTagNext,this_ent[0]->pos);
1803 
1804  /*
1805  * Update
1806  */
1807  t_pStreamInOut->dir() = new_dir;
1808 
1809  // Finished fiff_insert_after
1810 
1811  return true;
1812 }
1813 
1814 //=============================================================================================================
1815 // STATIC DEFINITIONS
1816 //=============================================================================================================
1817 
1818 //=============================================================================================================
1819 // DEFINE MEMBER METHODS
1820 //=============================================================================================================
1821 
1822 ComputeFwd::ComputeFwd(ComputeFwdSettings::SPtr pSettings)
1823  : sol(new FiffNamedMatrix)
1824  , sol_grad(new FiffNamedMatrix)
1825  , m_meg_forward(new FiffNamedMatrix)
1826  , m_meg_forward_grad(new FiffNamedMatrix)
1827  , m_eeg_forward(new FiffNamedMatrix)
1828  , m_eeg_forward_grad(new FiffNamedMatrix)
1829  , m_pSettings(pSettings)
1830 {
1831  initFwd();
1832 }
1833 
1834 //=============================================================================================================
1835 
1836 ComputeFwd::~ComputeFwd()
1837 {
1838  //ToDo Garbage collection
1839  for (int k = 0; k < m_iNSpace; k++)
1840  if(m_spaces[k])
1841  delete m_spaces[k];
1842  if(m_mri_head_t)
1843  delete m_mri_head_t;
1844  if(m_meg_head_t)
1845  delete m_meg_head_t;
1846  if(m_megcoils)
1847  delete m_megcoils;
1848  if(m_eegels)
1849  delete m_eegels;
1850  if(m_eegModel)
1851  delete m_eegModel;
1852 }
1853 
1854 //=============================================================================================================
1855 
1856 void ComputeFwd::initFwd()
1857 {
1858  // TODO: This only temporary until we have the fwd dlibrary refactored. This is only done in order to provide easy testing in test_forward_solution.
1859  m_spaces = Q_NULLPTR;
1860  m_iNSpace = 0;
1861  m_iNSource = 0;
1862 
1863  m_mri_head_t = Q_NULLPTR;
1864  m_meg_head_t = Q_NULLPTR;
1865 
1866  m_listMegChs = QList<FiffChInfo>();
1867  m_listEegChs = QList<FiffChInfo>();
1868  m_listCompChs = QList<FiffChInfo>();
1869 
1870  if(!m_pSettings->compute_grad) {
1871  m_meg_forward_grad = Q_NULLPTR;
1872  m_eeg_forward_grad = Q_NULLPTR;
1873  }
1874 
1875  int iNMeg = 0;
1876  int iNEeg = 0;
1877  int iNComp = 0;
1878 
1879  m_templates = Q_NULLPTR;
1880  m_megcoils = Q_NULLPTR;
1881  m_compcoils = Q_NULLPTR;
1882  m_compData = Q_NULLPTR;
1883  m_eegels = Q_NULLPTR;
1884  m_eegModels = Q_NULLPTR;
1885  m_iNChan = 0;
1886 
1887  int k;
1888  m_mri_id = Q_NULLPTR;
1889  m_meas_id.clear();
1890 
1891  FILE* out = Q_NULLPTR;
1893  m_eegModel = Q_NULLPTR;
1894  m_bemModel = Q_NULLPTR;
1895 
1896  // Report the setup
1897 
1898  // printf("\n");
1899  // mne_print_version_info(stderr,argv[0],PROGRAM_VERSION,__DATE__,__TIME__);
1900  printf("\n");
1901  printf("Source space : %s\n",m_pSettings->srcname.toUtf8().constData());
1902  if (!(m_pSettings->transname.isEmpty()) || !(m_pSettings->mriname.isEmpty())) {
1903  printf("MRI -> head transform source : %s\n",!(m_pSettings->mriname.isEmpty()) ? m_pSettings->mriname.toUtf8().constData() : m_pSettings->transname.toUtf8().constData());
1904  } else {
1905  printf("MRI and head coordinates are assumed to be identical.\n");
1906  }
1907  printf("Measurement data : %s\n",m_pSettings->measname.toUtf8().constData());
1908  if (!m_pSettings->bemname.isEmpty()) {
1909  printf("BEM model : %s\n",m_pSettings->bemname.toUtf8().constData());
1910  } else {
1911  printf("Sphere model : origin at (% 7.2f % 7.2f % 7.2f) mm\n",
1912  1000.0f*m_pSettings->r0[X_41],1000.0f*m_pSettings->r0[Y_41],1000.0f*m_pSettings->r0[Z_41]);
1913  if (m_pSettings->include_eeg) {
1914  printf("\n");
1915 
1916  if (m_pSettings->eeg_model_file.isEmpty()) {
1917  qCritical("!!!!!!!!!!TODO: default_eeg_model_file();");
1918  // m_pSettings->eeg_model_file = default_eeg_model_file();
1919  }
1920  m_eegModels = FwdEegSphereModelSet::fwd_load_eeg_sphere_models(m_pSettings->eeg_model_file,m_eegModels);
1921  m_eegModels->fwd_list_eeg_sphere_models(stderr);
1922 
1923  if (m_pSettings->eeg_model_name.isEmpty()) {
1924  m_pSettings->eeg_model_name = QString("Default");
1925  }
1926  if ((m_eegModel = m_eegModels->fwd_select_eeg_sphere_model(m_pSettings->eeg_model_name)) == Q_NULLPTR) {
1927  return;
1928  }
1929 
1930  if (!m_eegModel->fwd_setup_eeg_sphere_model(m_pSettings->eeg_sphere_rad,m_pSettings->use_equiv_eeg,3)) {
1931  return;
1932  }
1933 
1934  printf("Using EEG sphere model \"%s\" with scalp radius %7.1f mm\n",
1935  m_pSettings->eeg_model_name.toUtf8().constData(),1000*m_pSettings->eeg_sphere_rad);
1936  printf("%s the electrode locations to scalp\n",m_pSettings->scale_eeg_pos ? "Scale" : "Do not scale");
1937 
1938  m_eegModel->scale_pos = m_pSettings->scale_eeg_pos;
1939  VEC_COPY_41(m_eegModel->r0,m_pSettings->r0);
1940 
1941  printf("\n");
1942  }
1943  }
1944  printf("%s field computations\n",m_pSettings->accurate ? "Accurate" : "Standard");
1945  printf("Do computations in %s coordinates.\n",FiffCoordTransOld::mne_coord_frame_name(m_pSettings->coord_frame));
1946  printf("%s source orientations\n",m_pSettings->fixed_ori ? "Fixed" : "Free");
1947  if (m_pSettings->compute_grad) {
1948  printf("Compute derivatives with respect to source location coordinates\n");
1949  }
1950  printf("Destination for the solution : %s\n",m_pSettings->solname.toUtf8().constData());
1951  if (m_pSettings->do_all) {
1952  printf("Calculate solution for all source locations.\n");
1953  }
1954  if (m_pSettings->nlabel > 0) {
1955  printf("Source space will be restricted to sources in %d labels\n",m_pSettings->nlabel);
1956  }
1957 
1958  // Read the source locations
1959 
1960  printf("\n");
1961  printf("Reading %s...\n",m_pSettings->srcname.toUtf8().constData());
1962  if (MneSurfaceOrVolume::mne_read_source_spaces(m_pSettings->srcname,&m_spaces,&m_iNSpace) != OK) {
1963  return;
1964  }
1965  for (k = 0, m_iNSource = 0; k < m_iNSpace; k++) {
1966  if (m_pSettings->do_all) {
1967  MneSurfaceOrVolume::enable_all_sources(m_spaces[k]);
1968  }
1969  m_iNSource += m_spaces[k]->nuse;
1970  }
1971  if (m_iNSource == 0) {
1972  qCritical("No sources are active in these source spaces. --all option should be used.");
1973  return;
1974  }
1975  printf("Read %d source spaces a total of %d active source locations\n", m_iNSpace,m_iNSource);
1976  if (MneSurfaceOrVolume::restrict_sources_to_labels(m_spaces,m_iNSpace,m_pSettings->labels,m_pSettings->nlabel) == FAIL) {
1977  return;
1978  }
1979 
1980  // Read the MRI -> head coordinate transformation
1981  printf("\n");
1982  if (!m_pSettings->mriname.isEmpty()) {
1983  if ((m_mri_head_t = FiffCoordTransOld::mne_read_mri_transform(m_pSettings->mriname)) == Q_NULLPTR) {
1984  return;
1985  }
1986  if ((m_mri_id = get_file_id(m_pSettings->mriname)) == Q_NULLPTR) {
1987  qCritical("Couln't read MRI file id (How come?)");
1988  return;
1989  }
1990  }
1991  else if (!m_pSettings->transname.isEmpty()) {
1992  FiffCoordTransOld* t;
1993  if ((t = FiffCoordTransOld::mne_read_FShead2mri_transform(m_pSettings->transname.toUtf8().data())) == Q_NULLPTR) {
1994  return;
1995  }
1996  m_mri_head_t = t->fiff_invert_transform();
1997  if(t) {
1998  delete t;
1999  }
2000  } else {
2001  m_mri_head_t = FiffCoordTransOld::mne_identity_transform(FIFFV_COORD_MRI,FIFFV_COORD_HEAD);
2002  }
2003  FiffCoordTransOld::mne_print_coord_transform(stdout,m_mri_head_t);
2004 
2005  // Read the channel information and the MEG device -> head coordinate transformation
2006  // replace mne_read_meg_comp_eeg_ch_info_41()
2007 
2008  if(!m_pSettings->pFiffInfo) {
2009 
2010  QFile measname(m_pSettings->measname);
2011  FIFFLIB::FiffDirNode::SPtr DirNode;
2012  FiffStream::SPtr pStream(new FiffStream(&measname));
2013  FIFFLIB::FiffInfo fiffInfo;
2014  if(!pStream->open()) {
2015  qCritical() << "Could not open Stream.";
2016  return;
2017  }
2018 
2019  //Get Fiff info
2020  if(!pStream->read_meas_info(pStream->dirtree(), fiffInfo, DirNode)){
2021  qCritical() << "Could not find the channel information.";
2022  return;
2023  }
2024  pStream->close();
2025  m_pInfoBase = QSharedPointer<FIFFLIB::FiffInfo>(new FiffInfo(fiffInfo));
2026  } else {
2027  m_pInfoBase = m_pSettings->pFiffInfo;
2028  }
2029  if(!m_pInfoBase) {
2030  qCritical ("ComputeFwd::initFwd(): no FiffInfo");
2031  return;
2032  }
2033  if (mne_read_meg_comp_eeg_ch_info_41(m_pInfoBase,
2034  m_listMegChs,
2035  iNMeg,
2036  m_listCompChs,
2037  iNComp,
2038  m_listEegChs,
2039  iNEeg,
2040  &m_meg_head_t,
2041  m_meas_id) != OK) {
2042  return;
2043  }
2044 
2045  m_iNChan = iNMeg + iNEeg;
2046 
2047  printf("\n");
2048  if (iNMeg > 0) {
2049  printf("Read %3d MEG channels from %s\n",iNMeg,m_pSettings->measname.toUtf8().constData());
2050  }
2051  if (iNComp > 0) {
2052  printf("Read %3d MEG compensation channels from %s\n",iNComp,m_pSettings->measname.toUtf8().constData());
2053  }
2054  if (iNEeg > 0) {
2055  printf("Read %3d EEG channels from %s\n",iNEeg,m_pSettings->measname.toUtf8().constData());
2056  }
2057  if (!m_pSettings->include_meg) {
2058  printf("MEG not requested. MEG channels omitted.\n");
2059  m_listMegChs.clear();
2060  m_listCompChs.clear();
2061  iNMeg = 0;
2062  iNComp = 0;
2063  }
2064  else
2065  FiffCoordTransOld::mne_print_coord_transform(stdout,m_meg_head_t);
2066  if (!m_pSettings->include_eeg) {
2067  printf("EEG not requested. EEG channels omitted.\n");
2068  m_listEegChs.clear();
2069  iNEeg = 0;
2070  } else {
2071  if (mne_check_chinfo(m_listEegChs,iNEeg) != OK) {
2072  return;
2073  }
2074  }
2075 
2076  // Create coil descriptions with transformation to head or MRI frame
2077 
2078  if (m_pSettings->include_meg) {
2079  qPath = QString(QCoreApplication::applicationDirPath() + "/../resources/general/coilDefinitions/coil_def.dat");
2080  file.setFileName(qPath);
2081  if ( !QCoreApplication::startingUp() ) {
2082  qPath = QCoreApplication::applicationDirPath() + QString("/../resources/general/coilDefinitions/coil_def.dat");
2083  } else if (!file.exists()) {
2084  qPath = "../resources/general/coilDefinitions/coil_def.dat";
2085  }
2086 
2087  m_templates = FwdCoilSet::read_coil_defs(qPath);
2088  if (!m_templates) {
2089  return;
2090  }
2091 
2092  // Compensation data
2093 
2094  if ((m_compData = MneCTFCompDataSet::mne_read_ctf_comp_data(m_pSettings->measname)) == Q_NULLPTR) {
2095  return;
2096  }
2097  // Compensation channel information may be needed
2098  if (m_compData->ncomp > 0) {
2099  printf("%d compensation data sets in %s\n",m_compData->ncomp,m_pSettings->measname.toUtf8().constData());
2100  } else {
2101  m_listCompChs.clear();
2102  iNComp = 0;
2103 
2104  if(m_compData) {
2105  delete m_compData;
2106  }
2107  m_compData = Q_NULLPTR;
2108  }
2109  }
2110  if (m_pSettings->coord_frame == FIFFV_COORD_MRI) {
2111  FiffCoordTransOld* head_mri_t = m_mri_head_t->fiff_invert_transform();
2112  FiffCoordTransOld* meg_mri_t = FiffCoordTransOld::fiff_combine_transforms(FIFFV_COORD_DEVICE,FIFFV_COORD_MRI,m_meg_head_t,head_mri_t);
2113  if (meg_mri_t == Q_NULLPTR) {
2114  return;
2115  }
2116  if ((m_megcoils = m_templates->create_meg_coils(m_listMegChs,
2117  iNMeg,
2118  m_pSettings->accurate ? FWD_COIL_ACCURACY_ACCURATE : FWD_COIL_ACCURACY_NORMAL,
2119  meg_mri_t)) == Q_NULLPTR) {
2120  return;
2121  }
2122  if (iNComp > 0) {
2123  if ((m_compcoils = m_templates->create_meg_coils(m_listCompChs,
2124  iNComp,
2125  FWD_COIL_ACCURACY_NORMAL,
2126  meg_mri_t)) == Q_NULLPTR) {
2127  return;
2128  }
2129  }
2130  if ((m_eegels = FwdCoilSet::create_eeg_els(m_listEegChs,
2131  iNEeg,
2132  head_mri_t)) == Q_NULLPTR) {
2133  return;
2134  }
2135 
2136  FREE_41(head_mri_t);
2137  printf("MRI coordinate coil definitions created.\n");
2138  } else {
2139  if ((m_megcoils = m_templates->create_meg_coils(m_listMegChs,
2140  iNMeg,
2141  m_pSettings->accurate ? FWD_COIL_ACCURACY_ACCURATE : FWD_COIL_ACCURACY_NORMAL,
2142  m_meg_head_t)) == Q_NULLPTR) {
2143  return;
2144  }
2145 
2146  if (iNComp > 0) {
2147  if ((m_compcoils = m_templates->create_meg_coils(m_listCompChs,
2148  iNComp,
2149  FWD_COIL_ACCURACY_NORMAL,m_meg_head_t)) == Q_NULLPTR) {
2150  return;
2151  }
2152  }
2153  if ((m_eegels = FwdCoilSet::create_eeg_els(m_listEegChs,
2154  iNEeg,
2155  Q_NULLPTR)) == Q_NULLPTR) {
2156  return;
2157  }
2158  printf("Head coordinate coil definitions created.\n");
2159  }
2160 
2161  // Transform the source spaces into the appropriate coordinates
2162 
2163  if (MneSurfaceOrVolume::mne_transform_source_spaces_to(m_pSettings->coord_frame,m_mri_head_t,m_spaces,m_iNSpace) != OK) {
2164  return;
2165  }
2166  printf("Source spaces are now in %s coordinates.\n",FiffCoordTransOld::mne_coord_frame_name(m_pSettings->coord_frame));
2167 
2168  // Prepare the BEM model if necessary
2169 
2170  if (!m_pSettings->bemname.isEmpty()) {
2171  QString bemsolname = FwdBemModel::fwd_bem_make_bem_sol_name(m_pSettings->bemname);
2172  // FREE(bemname);
2173  m_pSettings->bemname = bemsolname;
2174 
2175  printf("\nSetting up the BEM model using %s...\n",m_pSettings->bemname.toUtf8().constData());
2176  printf("\nLoading surfaces...\n");
2177  m_bemModel = FwdBemModel::fwd_bem_load_three_layer_surfaces(m_pSettings->bemname);
2178 
2179  if (m_bemModel) {
2180  printf("Three-layer model surfaces loaded.\n");
2181  }
2182  else {
2183  m_bemModel = FwdBemModel::fwd_bem_load_homog_surface(m_pSettings->bemname);
2184  if (!m_bemModel) {
2185  return;
2186  }
2187  printf("Homogeneous model surface loaded.\n");
2188  }
2189  if (iNEeg > 0 && m_bemModel->nsurf == 1) {
2190  qCritical("Cannot use a homogeneous model in EEG calculations.");
2191  return;
2192  }
2193  printf("\nLoading the solution matrix...\n");
2194  if (FwdBemModel::fwd_bem_load_recompute_solution(m_pSettings->bemname.toUtf8().data(),FWD_BEM_UNKNOWN,FALSE,m_bemModel) == FAIL) {
2195  return;
2196  }
2197  if (m_pSettings->coord_frame == FIFFV_COORD_HEAD) {
2198  printf("Employing the head->MRI coordinate transform with the BEM model.\n");
2199  if (FwdBemModel::fwd_bem_set_head_mri_t(m_bemModel,m_mri_head_t) == FAIL) {
2200  return;
2201  }
2202  }
2203  printf("BEM model %s is now set up\n",m_bemModel->sol_name.toUtf8().constData());
2204  } else {
2205  printf("Using the sphere model.\n");
2206  }
2207  printf ("\n");
2208 
2209  // Try to circumvent numerical problems by excluding points too close our ouside the inner skull surface
2210 
2211  if (m_pSettings->filter_spaces) {
2212  if (!m_pSettings->mindistoutname.isEmpty()) {
2213  out = fopen(m_pSettings->mindistoutname.toUtf8().constData(),"w");
2214  if (out == Q_NULLPTR) {
2215  qCritical() << m_pSettings->mindistoutname;
2216  return;
2217  }
2218  printf("Omitted source space points will be output to : %s\n",m_pSettings->mindistoutname.toUtf8().constData());
2219  }
2220  MneSurfaceOrVolume::filter_source_spaces(m_pSettings->mindist,
2221  m_pSettings->bemname.toUtf8().data(),
2222  m_mri_head_t,
2223  m_spaces,
2224  m_iNSpace,out,m_pSettings->use_threads);
2225  if(out != Q_NULLPTR)
2226  fclose(out);
2227  }
2228 }
2229 
2230 //=========================================================================================================
2231 
2232 void ComputeFwd::calculateFwd()
2233 {
2234  int iNMeg = 0;
2235  int iNEeg = 0;
2236 
2237  if(m_megcoils) {
2238  iNMeg = m_megcoils->ncoil;
2239  }
2240  if(m_eegels) {
2241  iNEeg = m_eegels->ncoil;
2242  }
2243  if (!m_bemModel) {
2244  m_pSettings->use_threads = false;
2245  }
2246 
2247  // check if source spaces are still in head space
2248  if(m_spaces[0]->coord_frame != FIFFV_COORD_HEAD) {
2249  if (MneSurfaceOrVolume::mne_transform_source_spaces_to(m_pSettings->coord_frame,m_mri_head_t,m_spaces,m_iNSpace) != OK) {
2250  return;
2251  }
2252  }
2253 
2254  // Do the actual computation
2255  if (iNMeg > 0) {
2256  if ((FwdBemModel::compute_forward_meg(m_spaces,
2257  m_iNSpace,
2258  m_megcoils,
2259  m_compcoils,
2260  m_compData,
2261  m_pSettings->fixed_ori,
2262  m_bemModel,
2263  &m_pSettings->r0,
2264  m_pSettings->use_threads,
2265  *m_meg_forward.data(),
2266  *m_meg_forward_grad.data(),
2267  m_pSettings->compute_grad)) == FAIL) {
2268  return;
2269  }
2270  }
2271  if (iNEeg > 0) {
2272  if ((FwdBemModel::compute_forward_eeg(m_spaces,
2273  m_iNSpace,
2274  m_eegels,
2275  m_pSettings->fixed_ori,
2276  m_bemModel,
2277  m_eegModel,
2278  m_pSettings->use_threads,
2279  *m_eeg_forward.data(),
2280  *m_eeg_forward_grad.data(),
2281  m_pSettings->compute_grad))== FAIL) {
2282  return;
2283  }
2284  }
2285  if(iNMeg > 0 && iNEeg > 0) {
2286  if(m_meg_forward->data.cols() != m_eeg_forward->data.cols()) {
2287  qWarning() << "The MEG and EEG forward solutions do not match";
2288  return;
2289  }
2290  sol->clear();
2291  sol->data = MatrixXd(m_meg_forward->nrow + m_eeg_forward->nrow, m_meg_forward->ncol);
2292  sol->data.block(0,0,m_meg_forward->nrow,m_meg_forward->ncol) = m_meg_forward->data;
2293  sol->data.block(m_meg_forward->nrow,0,m_eeg_forward->nrow,m_eeg_forward->ncol) = m_eeg_forward->data;
2294  sol->nrow = m_meg_forward->nrow + m_eeg_forward->nrow;
2295  sol->row_names.append(m_meg_forward->row_names);
2296  } else if (iNMeg > 0) {
2297  sol = m_meg_forward;
2298  } else {
2299  sol = m_eeg_forward;
2300  }
2301 
2302  if(m_pSettings->compute_grad) {
2303  if(iNMeg > 0 && iNEeg > 0) {
2304  if(m_meg_forward_grad->data.cols() != m_eeg_forward_grad->data.cols()) {
2305  qWarning() << "The MEG and EEG forward solutions do not match";
2306  return;
2307  }
2308  sol_grad->clear();
2309  sol_grad->data = MatrixXd(m_meg_forward_grad->nrow + m_eeg_forward_grad->nrow, m_meg_forward_grad->ncol);
2310  sol_grad->data.block(0,0,m_meg_forward_grad->nrow,m_meg_forward_grad->ncol) = m_meg_forward_grad->data;
2311  sol_grad->data.block(m_meg_forward_grad->nrow,0,m_eeg_forward_grad->nrow,m_eeg_forward_grad->ncol) = m_eeg_forward_grad->data;
2312  sol_grad->nrow = m_meg_forward_grad->nrow + m_eeg_forward_grad->nrow;
2313  sol_grad->row_names.append(m_meg_forward_grad->row_names);
2314  } else if (iNMeg > 0) {
2315  sol_grad = m_meg_forward_grad;
2316  } else {
2317  sol_grad = m_eeg_forward_grad;
2318  }
2319  }
2320 }
2321 
2322 //=========================================================================================================
2323 
2324 void ComputeFwd::updateHeadPos(FiffCoordTransOld* transDevHeadOld)
2325 {
2326 
2327  int iNMeg = 0;
2328  if(m_megcoils) {
2329  iNMeg = m_megcoils->ncoil;
2330  }
2331 
2332  int iNComp = 0;
2333  if(m_compcoils) {
2334  iNComp = m_compcoils->ncoil;
2335  }
2336 // // transformation in head space
2337 // FiffCoordTransOld* transHeadHeadOld = new FiffCoordTransOld;
2338 //
2339 // // get transformation matrix in Headspace between two meg -> head transformation
2340 // // ToDo: the following part has to be tested
2341 // if(transDevHeadOld->from != m_megcoils->coord_frame) {
2342 // if(transDevHeadOld->to != m_megcoils->coord_frame) {
2343 // qWarning() << "ComputeFwd::updateHeadPos: Target Space not in Head Space - Space: " << transDevHeadOld->to;
2344 // return;
2345 // }
2346 // // calclulate rotation
2347 // Matrix3f matRot = m_meg_head_t->rot; // meg->head to origin
2348 // Matrix3f matRotDest = transDevHeadOld->rot; // meg->head after movement
2349 // Matrix3f matRotNew = Matrix3f::Zero(3,3); // rotation between origin and updated point
2350 //
2351 // // get Quaternions
2352 // Quaternionf quat(matRot);
2353 // Quaternionf quatDest(matRotDest);
2354 // Quaternionf quatNew;
2355 // // Rotation between two rotation matrices is calculated by multiplying with the inverse
2356 // // this is computitional easier with quaternions because for inverse we just have to flip signs
2357 // quatNew = quat*quatDest.inverse();
2358 // matRotNew = quatNew.toRotationMatrix(); // rotation between origin and updated point
2359 //
2360 // // calculate translation
2361 // // translation between to points = vector rsulting from substraction
2362 // VectorXf vecTrans = m_meg_head_t->move;
2363 // VectorXf vecTransDest = transDevHeadOld->move;
2364 // VectorXf vecTransNew = vecTrans - vecTransDest;
2365 //
2366 // // update new transformation matrix
2367 // transHeadHeadOld->from = m_megcoils->coord_frame;
2368 // transHeadHeadOld->to = m_megcoils->coord_frame;
2369 // transHeadHeadOld->rot = matRotNew;
2370 // transHeadHeadOld->move = vecTransNew;
2371 // transHeadHeadOld->add_inverse(transHeadHeadOld);
2372 // }
2373 //
2374 // FwdCoilSet* megcoilsNew = m_megcoils->dup_coil_set(transHeadHeadOld);
2375 
2376  // create new coilset with updated head position
2377  if (m_pSettings->coord_frame == FIFFV_COORD_MRI) {
2378  FiffCoordTransOld* head_mri_t = m_mri_head_t->fiff_invert_transform();
2379  FiffCoordTransOld* meg_mri_t = FiffCoordTransOld::fiff_combine_transforms(FIFFV_COORD_DEVICE,FIFFV_COORD_MRI,transDevHeadOld,head_mri_t);
2380  if (meg_mri_t == Q_NULLPTR) {
2381  return;
2382  }
2383  if ((m_megcoils = m_templates->create_meg_coils(m_listMegChs,
2384  iNMeg,
2385  m_pSettings->accurate ? FWD_COIL_ACCURACY_ACCURATE : FWD_COIL_ACCURACY_NORMAL,
2386  meg_mri_t)) == Q_NULLPTR) {
2387  return;
2388  }
2389  if (iNComp > 0) {
2390  if ((m_compcoils = m_templates->create_meg_coils(m_listCompChs,
2391  iNComp,
2392  FWD_COIL_ACCURACY_NORMAL,
2393  meg_mri_t)) == Q_NULLPTR) {
2394  return;
2395  }
2396  }
2397  } else {
2398  if ((m_megcoils = m_templates->create_meg_coils(m_listMegChs,
2399  iNMeg,
2400  m_pSettings->accurate ? FWD_COIL_ACCURACY_ACCURATE : FWD_COIL_ACCURACY_NORMAL,
2401  transDevHeadOld)) == Q_NULLPTR) {
2402  return;
2403  }
2404 
2405  if (iNComp > 0) {
2406  if ((m_compcoils = m_templates->create_meg_coils(m_listCompChs,
2407  iNComp,
2408  FWD_COIL_ACCURACY_NORMAL,transDevHeadOld)) == Q_NULLPTR) {
2409  return;
2410  }
2411  }
2412  }
2413 
2414  // check if source spaces are still in head space
2415  if(m_spaces[0]->coord_frame != FIFFV_COORD_HEAD) {
2416  if (MneSurfaceOrVolume::mne_transform_source_spaces_to(m_pSettings->coord_frame,m_mri_head_t,m_spaces,m_iNSpace) != OK) {
2417  return;
2418  }
2419  }
2420 
2421  // recompute meg forward
2422  if ((FwdBemModel::compute_forward_meg(m_spaces,
2423  m_iNSpace,
2424  m_megcoils,
2425  m_compcoils,
2426  m_compData, // we might have to update this too
2427  m_pSettings->fixed_ori,
2428  m_bemModel,
2429  &m_pSettings->r0,
2430  m_pSettings->use_threads,
2431  *m_meg_forward.data(),
2432  *m_meg_forward_grad.data(),
2433  m_pSettings->compute_grad)) == FAIL) {
2434  return;
2435  }
2436 
2437  // Update new Transformation Matrix
2438  m_meg_head_t = new FiffCoordTransOld(*transDevHeadOld);
2439  // update solution
2440  sol->data.block(0,0,m_meg_forward->nrow,m_meg_forward->ncol) = m_meg_forward->data;
2441  if(m_pSettings->compute_grad) {
2442  sol_grad->data.block(0,0,m_meg_forward_grad->nrow,m_meg_forward_grad->ncol) = m_meg_forward_grad->data;
2443  }
2444 }
2445 
2446 //=========================================================================================================
2447 
2448 void ComputeFwd::storeFwd(const QString& sSolName)
2449 {
2450  // We are ready to spill it out
2451  // Transform the source spaces back into MRI coordinates
2452  if (MneSourceSpaceOld::mne_transform_source_spaces_to(FIFFV_COORD_MRI,m_mri_head_t,m_spaces,m_iNSpace) != OK) {
2453  return;
2454  }
2455 
2456  int iNMeg = m_megcoils->ncoil;
2457  int iNEeg = m_eegels->ncoil;
2458 
2459  QString sName;
2460 
2461  if(sSolName == "default") {
2462  sName = m_pSettings->solname;
2463  } else {
2464  sName = sSolName;
2465  }
2466 
2467  printf("\nwriting %s...",sSolName.toUtf8().constData());
2468 
2469  if (!write_solution(sName, /* Destination file */
2470  m_spaces, /* The source spaces */
2471  m_iNSpace,
2472  m_pSettings->mriname,m_mri_id, /* MRI file and data obtained from there */
2473  m_mri_head_t,
2474  m_pSettings->measname,m_meas_id, /* MEG file and data obtained from there */
2475  m_meg_head_t,
2476  m_listMegChs,
2477  iNMeg,
2478  m_listEegChs,
2479  iNEeg,
2480  m_pSettings->fixed_ori, /* Fixed orientation dipoles? */
2481  m_pSettings->coord_frame, /* Coordinate frame */
2482  *m_meg_forward.data(),
2483  *m_eeg_forward.data(),
2484  *m_meg_forward_grad.data(),
2485  *m_eeg_forward_grad.data(),
2486  m_pSettings->compute_grad)) {
2487  return;
2488  }
2489 
2490  if (!mne_attach_env(m_pSettings->solname,m_pSettings->command)) {
2491  return;
2492  }
2493  printf("done\n");
2494  printf("\nFinished.\n");
2495 }
FIFFLIB::FiffInfoBase::SPtr
QSharedPointer< FiffInfoBase > SPtr
Definition: fiff_info_base.h:78
FWDLIB::FwdBemModel::fwd_bem_make_bem_sol_name
static QString fwd_bem_make_bem_sol_name(const QString &name)
Definition: fwd_bem_model.cpp:468
FIFF_MNE_SOURCE_SPACE_TYPE
#define FIFF_MNE_SOURCE_SPACE_TYPE
Definition: fiff_constants.h:349
FWDLIB::ComputeFwd::m_meg_forward_grad
QSharedDataPointer< FIFFLIB::FiffNamedMatrix > m_meg_forward_grad
Definition: compute_fwd.h:148
mne_source_space_old.h
MneSourceSpaceOld class declaration.
FWDLIB::ComputeFwd::storeFwd
void storeFwd(const QString &sSolName="default")
Definition: compute_fwd.cpp:2448
FIFF_MNE_SOURCE_SPACE_USE_TRIANGLES
#define FIFF_MNE_SOURCE_SPACE_USE_TRIANGLES
Definition: fiff_constants.h:354
FIFFLIB::_fiffIdRec
ToDo Old implementation use new fiff_id.h instead.
Definition: fiff_types_mne-c.h:214
FIFFLIB::FiffNamedMatrix::transpose_named_matrix
void transpose_named_matrix()
Definition: fiff_named_matrix.cpp:96
FIFFLIB::_fiffIdRec::time
fiffTimeRec time
Definition: fiff_types_mne-c.h:217
MNELIB::MneNearest
This is used in the patch definitions.
Definition: mne_nearest.h:77
FWDLIB::FwdBemModel::compute_forward_meg
static int compute_forward_meg(MNELIB::MneSourceSpaceOld **spaces, int nspace, FwdCoilSet *coils, FwdCoilSet *comp_coils, MNELIB::MneCTFCompDataSet *comp_data, bool fixed_ori, FwdBemModel *bem_model, Eigen::Vector3f *r0, bool use_threads, FIFFLIB::FiffNamedMatrix &resp, FIFFLIB::FiffNamedMatrix &resp_grad, bool bDoGRad)
Definition: fwd_bem_model.cpp:3100
FWDLIB::FwdCoilSet::create_meg_coils
FwdCoilSet * create_meg_coils(const QList< FIFFLIB::FiffChInfo > &chs, int nch, int acc, const FIFFLIB::FiffCoordTransOld *t)
Definition: fwd_coil_set.cpp:363
FIFFLIB::FiffCoordTransOld::to
FIFFLIB::fiff_int_t to
Definition: fiff_coord_trans_old.h:196
FIFFLIB::FiffStream::SPtr
QSharedPointer< FiffStream > SPtr
Definition: fiff_stream.h:107
FIFFLIB::FiffDirEntry::SPtr
QSharedPointer< FiffDirEntry > SPtr
Definition: fiff_dir_entry.h:70
FWDLIB::ComputeFwd::~ComputeFwd
virtual ~ComputeFwd()
Definition: compute_fwd.cpp:1836
fiff_sparse_matrix.h
FiffSparseMatrix class declaration.
FIFFLIB::FiffInfo
FIFF measurement file information.
Definition: fiff_info.h:84
FIFFLIB::FiffSparseMatrix
Data associated with MNE computations for each mneMeasDataSet.
Definition: fiff_sparse_matrix.h:74
FWDLIB::ComputeFwd::m_eeg_forward_grad
QSharedDataPointer< FIFFLIB::FiffNamedMatrix > m_eeg_forward_grad
Definition: compute_fwd.h:150
FIFFLIB::FiffDirNode::SPtr
QSharedPointer< FiffDirNode > SPtr
Definition: fiff_dir_node.h:76
FIFF_SUBJ_HIS_ID
#define FIFF_SUBJ_HIS_ID
Definition: fiff_file.h:573
FWDLIB::FwdEegSphereModelSet::fwd_list_eeg_sphere_models
void fwd_list_eeg_sphere_models(FILE *f)
Definition: fwd_eeg_sphere_model_set.cpp:260
FWDLIB::FwdBemModel::compute_forward_eeg
static int compute_forward_eeg(MNELIB::MneSourceSpaceOld **spaces, int nspace, FwdCoilSet *els, bool fixed_ori, FwdBemModel *bem_model, FwdEegSphereModel *m, bool use_threads, FIFFLIB::FiffNamedMatrix &resp, FIFFLIB::FiffNamedMatrix &resp_grad, bool bDoGrad)
Definition: fwd_bem_model.cpp:3364
FWDLIB::FwdEegSphereModel::fwd_setup_eeg_sphere_model
bool fwd_setup_eeg_sphere_model(float rad, bool fit_berg_scherg, int nfit)
Definition: fwd_eeg_sphere_model.cpp:1006
FIFFLIB::FiffSparseMatrix::coding
FIFFLIB::fiff_int_t coding
Definition: fiff_sparse_matrix.h:125
FIFF_MNE_SOURCE_SPACE_POINTS
#define FIFF_MNE_SOURCE_SPACE_POINTS
Definition: fiff_constants.h:341
FIFF_MNE_SOURCE_SPACE_NUSE
#define FIFF_MNE_SOURCE_SPACE_NUSE
Definition: fiff_constants.h:345
FIFFLIB::_fiffTimeRec::usecs
fiff_int_t usecs
Definition: fiff_types_mne-c.h:184
FIFF_MNE_SOURCE_SPACE_DIST
#define FIFF_MNE_SOURCE_SPACE_DIST
Definition: fiff_constants.h:360
FIFFLIB::FiffCoordTransOld::from
FIFFLIB::fiff_int_t from
Definition: fiff_coord_trans_old.h:195
FIFFLIB::FiffCoordTransOld::invrot
Eigen::Matrix3f invrot
Definition: fiff_coord_trans_old.h:199
FIFFLIB::FiffCoordTransOld::rot
Eigen::Matrix3f rot
Definition: fiff_coord_trans_old.h:197
fiff.h
FIFF class declaration, which provides static wrapper functions to stay consistent with mne matlab to...
FIFFLIB::FiffId
Universially unique identifier.
Definition: fiff_id.h:68
FIFF_MNE_SOURCE_SPACE_DIST_LIMIT
#define FIFF_MNE_SOURCE_SPACE_DIST_LIMIT
Definition: fiff_constants.h:361
FWDLIB::FwdEegSphereModel::r0
Eigen::Vector3f r0
Definition: fwd_eeg_sphere_model.h:328
k
int k
Definition: fiff_tag.cpp:322
FIFF_MNE_SOURCE_SPACE_NTRI
#define FIFF_MNE_SOURCE_SPACE_NTRI
Definition: fiff_constants.h:351
FIFFV_REF_MEG_CH
#define FIFFV_REF_MEG_CH
Definition: fiff_constants.h:270
FWDLIB::ComputeFwd::m_meg_forward
QSharedDataPointer< FIFFLIB::FiffNamedMatrix > m_meg_forward
Definition: compute_fwd.h:147
FIFFLIB::FiffCoordTransOld
Coordinate transformation descriptor.
Definition: fiff_coord_trans_old.h:80
FIFFB_BEM
#define FIFFB_BEM
Definition: fiff_file.h:403
FIFFLIB::FiffDirEntry
Directory entry description.
Definition: fiff_dir_entry.h:67
FIFFLIB::FiffSparseMatrix::ptrs
FIFFLIB::fiff_int_t * ptrs
Definition: fiff_sparse_matrix.h:131
FIFF_MNE_SOURCE_SPACE_SELECTION
#define FIFF_MNE_SOURCE_SPACE_SELECTION
Definition: fiff_constants.h:344
FIFFLIB::_fiffIdRec::machid
fiff_int_t machid[2]
Definition: fiff_types_mne-c.h:216
FWDLIB::FwdEegSphereModel::scale_pos
int scale_pos
Definition: fwd_eeg_sphere_model.h:336
FIFF_MNE_SOURCE_SPACE_MRI_FILE
#define FIFF_MNE_SOURCE_SPACE_MRI_FILE
Definition: fiff_constants.h:358
FIFF_MNE_COORD_FRAME
#define FIFF_MNE_COORD_FRAME
Definition: fiff_constants.h:331
FIFFLIB::FiffId::isEmpty
bool isEmpty() const
Definition: fiff_id.h:175
FWDLIB::ComputeFwdSettings::SPtr
QSharedPointer< ComputeFwdSettings > SPtr
Definition: compute_fwd_settings.h:82
FIFF_COORD_TRANS
#define FIFF_COORD_TRANS
Definition: fiff_file.h:475
FWDLIB::ComputeFwd::sol
QSharedDataPointer< FIFFLIB::FiffNamedMatrix > sol
Definition: compute_fwd.h:144
FIFF_MNE_SOURCE_SPACE_NORMALS
#define FIFF_MNE_SOURCE_SPACE_NORMALS
Definition: fiff_constants.h:342
FWDLIB::FwdEegSphereModelSet::fwd_load_eeg_sphere_models
static FwdEegSphereModelSet * fwd_load_eeg_sphere_models(const QString &p_sFileName, FwdEegSphereModelSet *now)
Definition: fwd_eeg_sphere_model_set.cpp:158
FIFFLIB::_fiffIdRec::version
fiff_int_t version
Definition: fiff_types_mne-c.h:215
FWDLIB::ComputeFwd::calculateFwd
void calculateFwd()
Definition: compute_fwd.cpp:2232
FWDLIB::ComputeFwd::m_eeg_forward
QSharedDataPointer< FIFFLIB::FiffNamedMatrix > m_eeg_forward
Definition: compute_fwd.h:149
FWDLIB::FwdCoilSet::read_coil_defs
static FwdCoilSet * read_coil_defs(const QString &name)
Definition: fwd_coil_set.cpp:416
FIFF_MNE_SOURCE_SPACE_NUSE_TRI
#define FIFF_MNE_SOURCE_SPACE_NUSE_TRI
Definition: fiff_constants.h:353
FIFFLIB::FiffTag::SPtr
QSharedPointer< FiffTag > SPtr
Definition: fiff_tag.h:152
mne_nearest.h
MneNearest class declaration.
FWDLIB::ComputeFwd::sol_grad
QSharedDataPointer< FIFFLIB::FiffNamedMatrix > sol_grad
Definition: compute_fwd.h:145
MNELIB::MneSourceSpaceOld
This defines a source space.
Definition: mne_source_space_old.h:76
FIFFLIB::FiffChInfo
Channel info descriptor.
Definition: fiff_ch_info.h:74
FIFFB_MRI
#define FIFFB_MRI
Definition: fiff_file.h:390
FWDLIB::FwdCoilSet::create_eeg_els
static FwdCoilSet * create_eeg_els(const QList< FIFFLIB::FiffChInfo > &chs, int nch, const FIFFLIB::FiffCoordTransOld *t)
Definition: fwd_coil_set.cpp:390
MNELIB::MneNamedMatrix
Matrix specification with a channel list.
Definition: mne_named_matrix.h:84
FIFF_MNE_ENV_WORKING_DIR
#define FIFF_MNE_ENV_WORKING_DIR
Definition: fiff_constants.h:401
FIFFLIB::FiffSparseMatrix::inds
FIFFLIB::fiff_int_t * inds
Definition: fiff_sparse_matrix.h:130
mne_forwardsolution.h
MNEForwardSolution class declaration, which provides the forward solution including the source space ...
FIFFLIB::_fiffTimeRec::secs
fiff_int_t secs
Definition: fiff_types_mne-c.h:183
FIFF_MNE_SOURCE_SPACE_ID
#define FIFF_MNE_SOURCE_SPACE_ID
Definition: fiff_constants.h:348
FIFF_NCHAN
#define FIFF_NCHAN
Definition: fiff_file.h:453
FIFFLIB::FiffStream
FIFF File I/O routines.
Definition: fiff_stream.h:104
mne_named_matrix.h
MNE Named Matrix (MneNamedMatrix) class declaration.
FIFF_MNE_SOURCE_SPACE_VOXEL_DIMS
#define FIFF_MNE_SOURCE_SPACE_VOXEL_DIMS
Definition: fiff_constants.h:356
FIFF_MNE_ENV_COMMAND_LINE
#define FIFF_MNE_ENV_COMMAND_LINE
Definition: fiff_constants.h:402
fiff_types.h
Definitions for describing the objects in a FIFF file.
FIFFLIB::FiffNamedMatrix
A named matrix.
Definition: fiff_named_matrix.h:76
FIFF_MNE_SOURCE_ORIENTATION
#define FIFF_MNE_SOURCE_ORIENTATION
Definition: fiff_constants.h:367
FIFF_MNE_SOURCE_SPACE_NEAREST
#define FIFF_MNE_SOURCE_SPACE_NEAREST
Definition: fiff_constants.h:346
FIFFLIB::FiffCoordTransOld::invmove
Eigen::Vector3f invmove
Definition: fiff_coord_trans_old.h:200
FWDLIB::ComputeFwd::updateHeadPos
void updateHeadPos(FIFFLIB::FiffCoordTransOld *transDevHeadOld)
Definition: compute_fwd.cpp:2324
FIFF_MNE_SOURCE_SPACE_TRIANGLES
#define FIFF_MNE_SOURCE_SPACE_TRIANGLES
Definition: fiff_constants.h:352
FIFFLIB::FiffSparseMatrix::m
FIFFLIB::fiff_int_t m
Definition: fiff_sparse_matrix.h:126
FIFFLIB::FiffId::clear
void clear()
Definition: fiff_id.cpp:115
FIFF_MNE_SOURCE_SPACE_NEAREST_DIST
#define FIFF_MNE_SOURCE_SPACE_NEAREST_DIST
Definition: fiff_constants.h:347
FIFFLIB::FiffCoordTransOld::move
Eigen::Vector3f move
Definition: fiff_coord_trans_old.h:198
FIFF_MNE_SOURCE_SPACE_NPOINTS
#define FIFF_MNE_SOURCE_SPACE_NPOINTS
Definition: fiff_constants.h:343
FWDLIB::FwdEegSphereModelSet::fwd_select_eeg_sphere_model
FwdEegSphereModel * fwd_select_eeg_sphere_model(const QString &p_sName)
Definition: fwd_eeg_sphere_model_set.cpp:234
fiff_coord_trans_old.h
FiffCoordTransOld class declaration.
FIFF_MNE_SOURCE_SPACE_INTERPOLATOR
#define FIFF_MNE_SOURCE_SPACE_INTERPOLATOR
Definition: fiff_constants.h:357
FIFFLIB::FiffSparseMatrix::n
FIFFLIB::fiff_int_t n
Definition: fiff_sparse_matrix.h:127
FIFF_MNE_FILE_NAME
#define FIFF_MNE_FILE_NAME
Definition: fiff_constants.h:336
mne_ctf_comp_data_set.h
MneCTFCompDataSet class declaration.
compute_fwd.h
Compute Forward Setting class declaration.
FIFFLIB::FiffSparseMatrix::data
FIFFLIB::fiff_float_t * data
Definition: fiff_sparse_matrix.h:129
FIFFLIB::FiffSparseMatrix::nz
FIFFLIB::fiff_int_t nz
Definition: fiff_sparse_matrix.h:128
Generated on Wed Aug 9 2023 17:55:41 for MNE-CPP by   doxygen 1.8.17