MNE-CPP 0.1.9
A Framework for Electrophysiology
Loading...
Searching...
No Matches
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
30using namespace Eigen;
31using namespace FWDLIB;
32using namespace FIFFLIB;
33using 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
76static 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
93float **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
110int **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
125void mne_free_cmatrix_41 (float **m)
126{
127 if (m) {
128 FREE_41(*m);
129 FREE_41(m);
130 }
131}
132
133void mne_free_icmatrix_41 (int **m)
134
135{
136 if (m) {
137 FREE_41(*m);
138 FREE_41(m);
139 }
140}
141
142fiffId 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
330int 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
367int 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
393void 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
438void 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
481static 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
495void 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
535void 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
609void 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
683int 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
841static 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
855void 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
863FiffSparseMatrix* 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
918bad : {
919 if(sparse)
920 delete sparse;
921 return Q_NULLPTR;
922 }
923}
924
925FiffSparseMatrix* 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
974out : {
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
986static 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
1105out : {
1106 FREE_41(inuse_map);
1107 FREE_41(nneighbors);
1108 FREE_41(neighbors);
1109 return res;
1110 }
1111}
1112
1113int 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
1293bad : {
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
1302QString 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
1320int 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
1348bool 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
1406bool 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
1637bad : {
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
1650bool 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
1822ComputeFwd::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
1836ComputeFwd::~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
1856void 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 int iNMeg = 0;
1871 int iNEeg = 0;
1872 int iNComp = 0;
1873
1874 m_templates = Q_NULLPTR;
1875 m_megcoils = Q_NULLPTR;
1876 m_compcoils = Q_NULLPTR;
1877 m_compData = Q_NULLPTR;
1878 m_eegels = Q_NULLPTR;
1879 m_eegModels = Q_NULLPTR;
1880 m_iNChan = 0;
1881
1882 int k;
1883 m_mri_id = Q_NULLPTR;
1884 m_meas_id.clear();
1885
1886 FILE* out = Q_NULLPTR;
1888 m_eegModel = Q_NULLPTR;
1889 m_bemModel = Q_NULLPTR;
1890
1891 // Report the setup
1892
1893 // printf("\n");
1894 // mne_print_version_info(stderr,argv[0],PROGRAM_VERSION,__DATE__,__TIME__);
1895 printf("\n");
1896 printf("Source space : %s\n",m_pSettings->srcname.toUtf8().constData());
1897 if (!(m_pSettings->transname.isEmpty()) || !(m_pSettings->mriname.isEmpty())) {
1898 printf("MRI -> head transform source : %s\n",!(m_pSettings->mriname.isEmpty()) ? m_pSettings->mriname.toUtf8().constData() : m_pSettings->transname.toUtf8().constData());
1899 } else {
1900 printf("MRI and head coordinates are assumed to be identical.\n");
1901 }
1902 printf("Measurement data : %s\n",m_pSettings->measname.toUtf8().constData());
1903 if (!m_pSettings->bemname.isEmpty()) {
1904 printf("BEM model : %s\n",m_pSettings->bemname.toUtf8().constData());
1905 } else {
1906 printf("Sphere model : origin at (% 7.2f % 7.2f % 7.2f) mm\n",
1907 1000.0f*m_pSettings->r0[X_41],1000.0f*m_pSettings->r0[Y_41],1000.0f*m_pSettings->r0[Z_41]);
1908 if (m_pSettings->include_eeg) {
1909 printf("\n");
1910
1911 if (m_pSettings->eeg_model_file.isEmpty()) {
1912 qCritical("!!!!!!!!!!TODO: default_eeg_model_file();");
1913 // m_pSettings->eeg_model_file = default_eeg_model_file();
1914 }
1915 m_eegModels = FwdEegSphereModelSet::fwd_load_eeg_sphere_models(m_pSettings->eeg_model_file,m_eegModels);
1916 m_eegModels->fwd_list_eeg_sphere_models(stderr);
1917
1918 if (m_pSettings->eeg_model_name.isEmpty()) {
1919 m_pSettings->eeg_model_name = QString("Default");
1920 }
1921 if ((m_eegModel = m_eegModels->fwd_select_eeg_sphere_model(m_pSettings->eeg_model_name)) == Q_NULLPTR) {
1922 return;
1923 }
1924
1925 if (!m_eegModel->fwd_setup_eeg_sphere_model(m_pSettings->eeg_sphere_rad,m_pSettings->use_equiv_eeg,3)) {
1926 return;
1927 }
1928
1929 printf("Using EEG sphere model \"%s\" with scalp radius %7.1f mm\n",
1930 m_pSettings->eeg_model_name.toUtf8().constData(),1000*m_pSettings->eeg_sphere_rad);
1931 printf("%s the electrode locations to scalp\n",m_pSettings->scale_eeg_pos ? "Scale" : "Do not scale");
1932
1933 m_eegModel->scale_pos = m_pSettings->scale_eeg_pos;
1934 VEC_COPY_41(m_eegModel->r0,m_pSettings->r0);
1935
1936 printf("\n");
1937 }
1938 }
1939 printf("%s field computations\n",m_pSettings->accurate ? "Accurate" : "Standard");
1940 printf("Do computations in %s coordinates.\n",FiffCoordTransOld::mne_coord_frame_name(m_pSettings->coord_frame));
1941 printf("%s source orientations\n",m_pSettings->fixed_ori ? "Fixed" : "Free");
1942 if (m_pSettings->compute_grad) {
1943 printf("Compute derivatives with respect to source location coordinates\n");
1944 }
1945 printf("Destination for the solution : %s\n",m_pSettings->solname.toUtf8().constData());
1946 if (m_pSettings->do_all) {
1947 printf("Calculate solution for all source locations.\n");
1948 }
1949 if (m_pSettings->nlabel > 0) {
1950 printf("Source space will be restricted to sources in %d labels\n",m_pSettings->nlabel);
1951 }
1952
1953 // Read the source locations
1954
1955 printf("\n");
1956 printf("Reading %s...\n",m_pSettings->srcname.toUtf8().constData());
1957 if (MneSurfaceOrVolume::mne_read_source_spaces(m_pSettings->srcname,&m_spaces,&m_iNSpace) != OK) {
1958 return;
1959 }
1960 for (k = 0, m_iNSource = 0; k < m_iNSpace; k++) {
1961 if (m_pSettings->do_all) {
1962 MneSurfaceOrVolume::enable_all_sources(m_spaces[k]);
1963 }
1964 m_iNSource += m_spaces[k]->nuse;
1965 }
1966 if (m_iNSource == 0) {
1967 qCritical("No sources are active in these source spaces. --all option should be used.");
1968 return;
1969 }
1970 printf("Read %d source spaces a total of %d active source locations\n", m_iNSpace,m_iNSource);
1971 if (MneSurfaceOrVolume::restrict_sources_to_labels(m_spaces,m_iNSpace,m_pSettings->labels,m_pSettings->nlabel) == FAIL) {
1972 return;
1973 }
1974
1975 // Read the MRI -> head coordinate transformation
1976 printf("\n");
1977 if (!m_pSettings->mriname.isEmpty()) {
1978 if ((m_mri_head_t = FiffCoordTransOld::mne_read_mri_transform(m_pSettings->mriname)) == Q_NULLPTR) {
1979 return;
1980 }
1981 if ((m_mri_id = get_file_id(m_pSettings->mriname)) == Q_NULLPTR) {
1982 qCritical("Couln't read MRI file id (How come?)");
1983 return;
1984 }
1985 }
1986 else if (!m_pSettings->transname.isEmpty()) {
1987 FiffCoordTransOld* t;
1988 if ((t = FiffCoordTransOld::mne_read_FShead2mri_transform(m_pSettings->transname.toUtf8().data())) == Q_NULLPTR) {
1989 return;
1990 }
1991 m_mri_head_t = t->fiff_invert_transform();
1992 if(t) {
1993 delete t;
1994 }
1995 } else {
1996 m_mri_head_t = FiffCoordTransOld::mne_identity_transform(FIFFV_COORD_MRI,FIFFV_COORD_HEAD);
1997 }
1998 FiffCoordTransOld::mne_print_coord_transform(stdout,m_mri_head_t);
1999
2000 // Read the channel information and the MEG device -> head coordinate transformation
2001 // replace mne_read_meg_comp_eeg_ch_info_41()
2002
2003 if(!m_pSettings->pFiffInfo) {
2004
2005 QFile measname(m_pSettings->measname);
2006 FIFFLIB::FiffDirNode::SPtr DirNode;
2007 FiffStream::SPtr pStream(new FiffStream(&measname));
2008 FIFFLIB::FiffInfo fiffInfo;
2009 if(!pStream->open()) {
2010 qCritical() << "Could not open Stream.";
2011 return;
2012 }
2013
2014 //Get Fiff info
2015 if(!pStream->read_meas_info(pStream->dirtree(), fiffInfo, DirNode)){
2016 qCritical() << "Could not find the channel information.";
2017 return;
2018 }
2019 pStream->close();
2020 m_pInfoBase = QSharedPointer<FIFFLIB::FiffInfo>(new FiffInfo(fiffInfo));
2021 } else {
2022 m_pInfoBase = m_pSettings->pFiffInfo;
2023 }
2024 if(!m_pInfoBase) {
2025 qCritical ("ComputeFwd::initFwd(): no FiffInfo");
2026 return;
2027 }
2028 if (mne_read_meg_comp_eeg_ch_info_41(m_pInfoBase,
2029 m_listMegChs,
2030 iNMeg,
2031 m_listCompChs,
2032 iNComp,
2033 m_listEegChs,
2034 iNEeg,
2035 &m_meg_head_t,
2036 m_meas_id) != OK) {
2037 return;
2038 }
2039
2040 m_iNChan = iNMeg + iNEeg;
2041
2042 printf("\n");
2043 if (iNMeg > 0) {
2044 printf("Read %3d MEG channels from %s\n",iNMeg,m_pSettings->measname.toUtf8().constData());
2045 }
2046 if (iNComp > 0) {
2047 printf("Read %3d MEG compensation channels from %s\n",iNComp,m_pSettings->measname.toUtf8().constData());
2048 }
2049 if (iNEeg > 0) {
2050 printf("Read %3d EEG channels from %s\n",iNEeg,m_pSettings->measname.toUtf8().constData());
2051 }
2052 if (!m_pSettings->include_meg) {
2053 printf("MEG not requested. MEG channels omitted.\n");
2054 m_listMegChs.clear();
2055 m_listCompChs.clear();
2056 iNMeg = 0;
2057 iNComp = 0;
2058 }
2059 else
2060 FiffCoordTransOld::mne_print_coord_transform(stdout,m_meg_head_t);
2061 if (!m_pSettings->include_eeg) {
2062 printf("EEG not requested. EEG channels omitted.\n");
2063 m_listEegChs.clear();
2064 iNEeg = 0;
2065 } else {
2066 if (mne_check_chinfo(m_listEegChs,iNEeg) != OK) {
2067 return;
2068 }
2069 }
2070
2071 // Create coil descriptions with transformation to head or MRI frame
2072
2073 if (m_pSettings->include_meg) {
2074 qPath = QString(QCoreApplication::applicationDirPath() + "/../resources/general/coilDefinitions/coil_def.dat");
2075 file.setFileName(qPath);
2076 if ( !QCoreApplication::startingUp() ) {
2077 qPath = QCoreApplication::applicationDirPath() + QString("/../resources/general/coilDefinitions/coil_def.dat");
2078 } else if (!file.exists()) {
2079 qPath = "../resources/general/coilDefinitions/coil_def.dat";
2080 }
2081
2082 m_templates = FwdCoilSet::read_coil_defs(qPath);
2083 if (!m_templates) {
2084 return;
2085 }
2086
2087 // Compensation data
2088
2089 if ((m_compData = MneCTFCompDataSet::mne_read_ctf_comp_data(m_pSettings->measname)) == Q_NULLPTR) {
2090 return;
2091 }
2092 // Compensation channel information may be needed
2093 if (m_compData->ncomp > 0) {
2094 printf("%d compensation data sets in %s\n",m_compData->ncomp,m_pSettings->measname.toUtf8().constData());
2095 } else {
2096 m_listCompChs.clear();
2097 iNComp = 0;
2098
2099 if(m_compData) {
2100 delete m_compData;
2101 }
2102 m_compData = Q_NULLPTR;
2103 }
2104 }
2105 if (m_pSettings->coord_frame == FIFFV_COORD_MRI) {
2106 FiffCoordTransOld* head_mri_t = m_mri_head_t->fiff_invert_transform();
2107 FiffCoordTransOld* meg_mri_t = FiffCoordTransOld::fiff_combine_transforms(FIFFV_COORD_DEVICE,FIFFV_COORD_MRI,m_meg_head_t,head_mri_t);
2108 if (meg_mri_t == Q_NULLPTR) {
2109 return;
2110 }
2111 if ((m_megcoils = m_templates->create_meg_coils(m_listMegChs,
2112 iNMeg,
2113 m_pSettings->accurate ? FWD_COIL_ACCURACY_ACCURATE : FWD_COIL_ACCURACY_NORMAL,
2114 meg_mri_t)) == Q_NULLPTR) {
2115 return;
2116 }
2117 if (iNComp > 0) {
2118 if ((m_compcoils = m_templates->create_meg_coils(m_listCompChs,
2119 iNComp,
2120 FWD_COIL_ACCURACY_NORMAL,
2121 meg_mri_t)) == Q_NULLPTR) {
2122 return;
2123 }
2124 }
2125 if ((m_eegels = FwdCoilSet::create_eeg_els(m_listEegChs,
2126 iNEeg,
2127 head_mri_t)) == Q_NULLPTR) {
2128 return;
2129 }
2130
2131 FREE_41(head_mri_t);
2132 printf("MRI coordinate coil definitions created.\n");
2133 } else {
2134 if ((m_megcoils = m_templates->create_meg_coils(m_listMegChs,
2135 iNMeg,
2136 m_pSettings->accurate ? FWD_COIL_ACCURACY_ACCURATE : FWD_COIL_ACCURACY_NORMAL,
2137 m_meg_head_t)) == Q_NULLPTR) {
2138 return;
2139 }
2140
2141 if (iNComp > 0) {
2142 if ((m_compcoils = m_templates->create_meg_coils(m_listCompChs,
2143 iNComp,
2144 FWD_COIL_ACCURACY_NORMAL,m_meg_head_t)) == Q_NULLPTR) {
2145 return;
2146 }
2147 }
2148 if ((m_eegels = FwdCoilSet::create_eeg_els(m_listEegChs,
2149 iNEeg,
2150 Q_NULLPTR)) == Q_NULLPTR) {
2151 return;
2152 }
2153 printf("Head coordinate coil definitions created.\n");
2154 }
2155
2156 // Transform the source spaces into the appropriate coordinates
2157
2158 if (MneSurfaceOrVolume::mne_transform_source_spaces_to(m_pSettings->coord_frame,m_mri_head_t,m_spaces,m_iNSpace) != OK) {
2159 return;
2160 }
2161 printf("Source spaces are now in %s coordinates.\n",FiffCoordTransOld::mne_coord_frame_name(m_pSettings->coord_frame));
2162
2163 // Prepare the BEM model if necessary
2164
2165 if (!m_pSettings->bemname.isEmpty()) {
2166 QString bemsolname = FwdBemModel::fwd_bem_make_bem_sol_name(m_pSettings->bemname);
2167 // FREE(bemname);
2168 m_pSettings->bemname = bemsolname;
2169
2170 printf("\nSetting up the BEM model using %s...\n",m_pSettings->bemname.toUtf8().constData());
2171 printf("\nLoading surfaces...\n");
2172 m_bemModel = FwdBemModel::fwd_bem_load_three_layer_surfaces(m_pSettings->bemname);
2173
2174 if (m_bemModel) {
2175 printf("Three-layer model surfaces loaded.\n");
2176 }
2177 else {
2178 m_bemModel = FwdBemModel::fwd_bem_load_homog_surface(m_pSettings->bemname);
2179 if (!m_bemModel) {
2180 return;
2181 }
2182 printf("Homogeneous model surface loaded.\n");
2183 }
2184 if (iNEeg > 0 && m_bemModel->nsurf == 1) {
2185 qCritical("Cannot use a homogeneous model in EEG calculations.");
2186 return;
2187 }
2188 printf("\nLoading the solution matrix...\n");
2189 if (FwdBemModel::fwd_bem_load_recompute_solution(m_pSettings->bemname.toUtf8().data(),FWD_BEM_UNKNOWN,FALSE,m_bemModel) == FAIL) {
2190 return;
2191 }
2192 if (m_pSettings->coord_frame == FIFFV_COORD_HEAD) {
2193 printf("Employing the head->MRI coordinate transform with the BEM model.\n");
2194 if (FwdBemModel::fwd_bem_set_head_mri_t(m_bemModel,m_mri_head_t) == FAIL) {
2195 return;
2196 }
2197 }
2198 printf("BEM model %s is now set up\n",m_bemModel->sol_name.toUtf8().constData());
2199 } else {
2200 printf("Using the sphere model.\n");
2201 }
2202 printf ("\n");
2203
2204 // Try to circumvent numerical problems by excluding points too close our ouside the inner skull surface
2205
2206 if (m_pSettings->filter_spaces) {
2207 if (!m_pSettings->mindistoutname.isEmpty()) {
2208 out = fopen(m_pSettings->mindistoutname.toUtf8().constData(),"w");
2209 if (out == Q_NULLPTR) {
2210 qCritical() << m_pSettings->mindistoutname;
2211 return;
2212 }
2213 printf("Omitted source space points will be output to : %s\n",m_pSettings->mindistoutname.toUtf8().constData());
2214 }
2215 MneSurfaceOrVolume::filter_source_spaces(m_pSettings->mindist,
2216 m_pSettings->bemname.toUtf8().data(),
2217 m_mri_head_t,
2218 m_spaces,
2219 m_iNSpace,out,m_pSettings->use_threads);
2220 if(out != Q_NULLPTR)
2221 fclose(out);
2222 }
2223}
2224
2225//=========================================================================================================
2226
2227void ComputeFwd::calculateFwd()
2228{
2229 int iNMeg = 0;
2230 int iNEeg = 0;
2231
2232 if(m_megcoils) {
2233 iNMeg = m_megcoils->ncoil;
2234 }
2235 if(m_eegels) {
2236 iNEeg = m_eegels->ncoil;
2237 }
2238 if (!m_bemModel) {
2239 m_pSettings->use_threads = false;
2240 }
2241
2242 // check if source spaces are still in head space
2243 if(m_spaces[0]->coord_frame != FIFFV_COORD_HEAD) {
2244 if (MneSurfaceOrVolume::mne_transform_source_spaces_to(m_pSettings->coord_frame,m_mri_head_t,m_spaces,m_iNSpace) != OK) {
2245 return;
2246 }
2247 }
2248
2249 // Do the actual computation
2250 if (iNMeg > 0) {
2251 if ((FwdBemModel::compute_forward_meg(m_spaces,
2252 m_iNSpace,
2253 m_megcoils,
2254 m_compcoils,
2255 m_compData,
2256 m_pSettings->fixed_ori,
2257 m_bemModel,
2258 &m_pSettings->r0,
2259 m_pSettings->use_threads,
2260 *m_meg_forward.data(),
2261 *m_meg_forward_grad.data(),
2262 m_pSettings->compute_grad)) == FAIL) {
2263 return;
2264 }
2265 }
2266 if (iNEeg > 0) {
2267 if ((FwdBemModel::compute_forward_eeg(m_spaces,
2268 m_iNSpace,
2269 m_eegels,
2270 m_pSettings->fixed_ori,
2271 m_bemModel,
2272 m_eegModel,
2273 m_pSettings->use_threads,
2274 *m_eeg_forward.data(),
2275 *m_eeg_forward_grad.data(),
2276 m_pSettings->compute_grad))== FAIL) {
2277 return;
2278 }
2279 }
2280 if(iNMeg > 0 && iNEeg > 0) {
2281 if(m_meg_forward->data.cols() != m_eeg_forward->data.cols()) {
2282 qWarning() << "The MEG and EEG forward solutions do not match";
2283 return;
2284 }
2285 sol->clear();
2286 sol->data = MatrixXd(m_meg_forward->nrow + m_eeg_forward->nrow, m_meg_forward->ncol);
2287 sol->data.block(0,0,m_meg_forward->nrow,m_meg_forward->ncol) = m_meg_forward->data;
2288 sol->data.block(m_meg_forward->nrow,0,m_eeg_forward->nrow,m_eeg_forward->ncol) = m_eeg_forward->data;
2289 sol->nrow = m_meg_forward->nrow + m_eeg_forward->nrow;
2290 sol->row_names.append(m_meg_forward->row_names);
2291 } else if (iNMeg > 0) {
2292 sol = m_meg_forward;
2293 } else {
2294 sol = m_eeg_forward;
2295 }
2296
2297 if(m_pSettings->compute_grad) {
2298 if(iNMeg > 0 && iNEeg > 0) {
2299 if(m_meg_forward_grad->data.cols() != m_eeg_forward_grad->data.cols()) {
2300 qWarning() << "The MEG and EEG forward solutions do not match";
2301 return;
2302 }
2303 sol_grad->clear();
2304 sol_grad->data = MatrixXd(m_meg_forward_grad->nrow + m_eeg_forward_grad->nrow, m_meg_forward_grad->ncol);
2305 sol_grad->data.block(0,0,m_meg_forward_grad->nrow,m_meg_forward_grad->ncol) = m_meg_forward_grad->data;
2306 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;
2307 sol_grad->nrow = m_meg_forward_grad->nrow + m_eeg_forward_grad->nrow;
2308 sol_grad->row_names.append(m_meg_forward_grad->row_names);
2309 } else if (iNMeg > 0) {
2310 sol_grad = m_meg_forward_grad;
2311 } else {
2312 sol_grad = m_eeg_forward_grad;
2313 }
2314 }
2315}
2316
2317//=========================================================================================================
2318
2319void ComputeFwd::updateHeadPos(FiffCoordTransOld* transDevHeadOld)
2320{
2321
2322 int iNMeg = 0;
2323 if(m_megcoils) {
2324 iNMeg = m_megcoils->ncoil;
2325 }
2326
2327 int iNComp = 0;
2328 if(m_compcoils) {
2329 iNComp = m_compcoils->ncoil;
2330 }
2331// // transformation in head space
2332// FiffCoordTransOld* transHeadHeadOld = new FiffCoordTransOld;
2333//
2334// // get transformation matrix in Headspace between two meg -> head transformation
2335// // ToDo: the following part has to be tested
2336// if(transDevHeadOld->from != m_megcoils->coord_frame) {
2337// if(transDevHeadOld->to != m_megcoils->coord_frame) {
2338// qWarning() << "ComputeFwd::updateHeadPos: Target Space not in Head Space - Space: " << transDevHeadOld->to;
2339// return;
2340// }
2341// // calclulate rotation
2342// Matrix3f matRot = m_meg_head_t->rot; // meg->head to origin
2343// Matrix3f matRotDest = transDevHeadOld->rot; // meg->head after movement
2344// Matrix3f matRotNew = Matrix3f::Zero(3,3); // rotation between origin and updated point
2345//
2346// // get Quaternions
2347// Quaternionf quat(matRot);
2348// Quaternionf quatDest(matRotDest);
2349// Quaternionf quatNew;
2350// // Rotation between two rotation matrices is calculated by multiplying with the inverse
2351// // this is computitional easier with quaternions because for inverse we just have to flip signs
2352// quatNew = quat*quatDest.inverse();
2353// matRotNew = quatNew.toRotationMatrix(); // rotation between origin and updated point
2354//
2355// // calculate translation
2356// // translation between to points = vector rsulting from substraction
2357// VectorXf vecTrans = m_meg_head_t->move;
2358// VectorXf vecTransDest = transDevHeadOld->move;
2359// VectorXf vecTransNew = vecTrans - vecTransDest;
2360//
2361// // update new transformation matrix
2362// transHeadHeadOld->from = m_megcoils->coord_frame;
2363// transHeadHeadOld->to = m_megcoils->coord_frame;
2364// transHeadHeadOld->rot = matRotNew;
2365// transHeadHeadOld->move = vecTransNew;
2366// transHeadHeadOld->add_inverse(transHeadHeadOld);
2367// }
2368//
2369// FwdCoilSet* megcoilsNew = m_megcoils->dup_coil_set(transHeadHeadOld);
2370
2371 // create new coilset with updated head position
2372 if (m_pSettings->coord_frame == FIFFV_COORD_MRI) {
2373 FiffCoordTransOld* head_mri_t = m_mri_head_t->fiff_invert_transform();
2374 FiffCoordTransOld* meg_mri_t = FiffCoordTransOld::fiff_combine_transforms(FIFFV_COORD_DEVICE,FIFFV_COORD_MRI,transDevHeadOld,head_mri_t);
2375 if (meg_mri_t == Q_NULLPTR) {
2376 return;
2377 }
2378 if ((m_megcoils = m_templates->create_meg_coils(m_listMegChs,
2379 iNMeg,
2380 m_pSettings->accurate ? FWD_COIL_ACCURACY_ACCURATE : FWD_COIL_ACCURACY_NORMAL,
2381 meg_mri_t)) == Q_NULLPTR) {
2382 return;
2383 }
2384 if (iNComp > 0) {
2385 if ((m_compcoils = m_templates->create_meg_coils(m_listCompChs,
2386 iNComp,
2387 FWD_COIL_ACCURACY_NORMAL,
2388 meg_mri_t)) == Q_NULLPTR) {
2389 return;
2390 }
2391 }
2392 } else {
2393 if ((m_megcoils = m_templates->create_meg_coils(m_listMegChs,
2394 iNMeg,
2395 m_pSettings->accurate ? FWD_COIL_ACCURACY_ACCURATE : FWD_COIL_ACCURACY_NORMAL,
2396 transDevHeadOld)) == Q_NULLPTR) {
2397 return;
2398 }
2399
2400 if (iNComp > 0) {
2401 if ((m_compcoils = m_templates->create_meg_coils(m_listCompChs,
2402 iNComp,
2403 FWD_COIL_ACCURACY_NORMAL,transDevHeadOld)) == Q_NULLPTR) {
2404 return;
2405 }
2406 }
2407 }
2408
2409 // check if source spaces are still in head space
2410 if(m_spaces[0]->coord_frame != FIFFV_COORD_HEAD) {
2411 if (MneSurfaceOrVolume::mne_transform_source_spaces_to(m_pSettings->coord_frame,m_mri_head_t,m_spaces,m_iNSpace) != OK) {
2412 return;
2413 }
2414 }
2415
2416 // recompute meg forward
2417 if ((FwdBemModel::compute_forward_meg(m_spaces,
2418 m_iNSpace,
2419 m_megcoils,
2420 m_compcoils,
2421 m_compData, // we might have to update this too
2422 m_pSettings->fixed_ori,
2423 m_bemModel,
2424 &m_pSettings->r0,
2425 m_pSettings->use_threads,
2426 *m_meg_forward.data(),
2427 *m_meg_forward_grad.data(),
2428 m_pSettings->compute_grad)) == FAIL) {
2429 return;
2430 }
2431
2432 // Update new Transformation Matrix
2433 m_meg_head_t = new FiffCoordTransOld(*transDevHeadOld);
2434 // update solution
2435 sol->data.block(0,0,m_meg_forward->nrow,m_meg_forward->ncol) = m_meg_forward->data;
2436 if(m_pSettings->compute_grad) {
2437 sol_grad->data.block(0,0,m_meg_forward_grad->nrow,m_meg_forward_grad->ncol) = m_meg_forward_grad->data;
2438 }
2439}
2440
2441//=========================================================================================================
2442
2443void ComputeFwd::storeFwd(const QString& sSolName)
2444{
2445 // We are ready to spill it out
2446 // Transform the source spaces back into MRI coordinates
2447 if (MneSourceSpaceOld::mne_transform_source_spaces_to(FIFFV_COORD_MRI,m_mri_head_t,m_spaces,m_iNSpace) != OK) {
2448 return;
2449 }
2450
2451 int iNMeg = m_megcoils->ncoil;
2452 int iNEeg = m_eegels->ncoil;
2453
2454 QString sName;
2455
2456 if(sSolName == "default") {
2457 sName = m_pSettings->solname;
2458 } else {
2459 sName = sSolName;
2460 }
2461
2462 printf("\nwriting %s...",sSolName.toUtf8().constData());
2463
2464 if (!write_solution(sName, /* Destination file */
2465 m_spaces, /* The source spaces */
2466 m_iNSpace,
2467 m_pSettings->mriname,m_mri_id, /* MRI file and data obtained from there */
2468 m_mri_head_t,
2469 m_pSettings->measname,m_meas_id, /* MEG file and data obtained from there */
2470 m_meg_head_t,
2471 m_listMegChs,
2472 iNMeg,
2473 m_listEegChs,
2474 iNEeg,
2475 m_pSettings->fixed_ori, /* Fixed orientation dipoles? */
2476 m_pSettings->coord_frame, /* Coordinate frame */
2477 *m_meg_forward.data(),
2478 *m_eeg_forward.data(),
2479 *m_meg_forward_grad.data(),
2480 *m_eeg_forward_grad.data(),
2481 m_pSettings->compute_grad)) {
2482 return;
2483 }
2484
2485 if (!mne_attach_env(m_pSettings->solname,m_pSettings->command)) {
2486 return;
2487 }
2488 printf("done\n");
2489 printf("\nFinished.\n");
2490}
fiff.h
FIFF class declaration, which provides static wrapper functions to stay consistent with mne matlab to...
FIFF_MNE_SOURCE_SPACE_ID
#define FIFF_MNE_SOURCE_SPACE_ID
Definition fiff_constants.h:348
FIFF_MNE_SOURCE_SPACE_DIST_LIMIT
#define FIFF_MNE_SOURCE_SPACE_DIST_LIMIT
Definition fiff_constants.h:361
FIFF_MNE_SOURCE_SPACE_VOXEL_DIMS
#define FIFF_MNE_SOURCE_SPACE_VOXEL_DIMS
Definition fiff_constants.h:356
FIFF_MNE_COORD_FRAME
#define FIFF_MNE_COORD_FRAME
Definition fiff_constants.h:331
FIFF_MNE_SOURCE_ORIENTATION
#define FIFF_MNE_SOURCE_ORIENTATION
Definition fiff_constants.h:367
FIFF_MNE_SOURCE_SPACE_TYPE
#define FIFF_MNE_SOURCE_SPACE_TYPE
Definition fiff_constants.h:349
FIFF_MNE_SOURCE_SPACE_SELECTION
#define FIFF_MNE_SOURCE_SPACE_SELECTION
Definition fiff_constants.h:344
FIFF_MNE_SOURCE_SPACE_USE_TRIANGLES
#define FIFF_MNE_SOURCE_SPACE_USE_TRIANGLES
Definition fiff_constants.h:354
FIFF_MNE_SOURCE_SPACE_NUSE_TRI
#define FIFF_MNE_SOURCE_SPACE_NUSE_TRI
Definition fiff_constants.h:353
FIFF_MNE_SOURCE_SPACE_INTERPOLATOR
#define FIFF_MNE_SOURCE_SPACE_INTERPOLATOR
Definition fiff_constants.h:357
FIFF_MNE_SOURCE_SPACE_NORMALS
#define FIFF_MNE_SOURCE_SPACE_NORMALS
Definition fiff_constants.h:342
FIFF_MNE_SOURCE_SPACE_NEAREST_DIST
#define FIFF_MNE_SOURCE_SPACE_NEAREST_DIST
Definition fiff_constants.h:347
FIFF_MNE_SOURCE_SPACE_DIST
#define FIFF_MNE_SOURCE_SPACE_DIST
Definition fiff_constants.h:360
FIFF_MNE_SOURCE_SPACE_POINTS
#define FIFF_MNE_SOURCE_SPACE_POINTS
Definition fiff_constants.h:341
FIFFV_REF_MEG_CH
#define FIFFV_REF_MEG_CH
Definition fiff_constants.h:270
FIFF_MNE_SOURCE_SPACE_NTRI
#define FIFF_MNE_SOURCE_SPACE_NTRI
Definition fiff_constants.h:351
FIFF_MNE_SOURCE_SPACE_MRI_FILE
#define FIFF_MNE_SOURCE_SPACE_MRI_FILE
Definition fiff_constants.h:358
FIFF_MNE_ENV_WORKING_DIR
#define FIFF_MNE_ENV_WORKING_DIR
Definition fiff_constants.h:401
FIFF_MNE_ENV_COMMAND_LINE
#define FIFF_MNE_ENV_COMMAND_LINE
Definition fiff_constants.h:402
FIFF_MNE_SOURCE_SPACE_NPOINTS
#define FIFF_MNE_SOURCE_SPACE_NPOINTS
Definition fiff_constants.h:343
FIFF_MNE_SOURCE_SPACE_TRIANGLES
#define FIFF_MNE_SOURCE_SPACE_TRIANGLES
Definition fiff_constants.h:352
FIFF_MNE_SOURCE_SPACE_NUSE
#define FIFF_MNE_SOURCE_SPACE_NUSE
Definition fiff_constants.h:345
FIFF_MNE_FILE_NAME
#define FIFF_MNE_FILE_NAME
Definition fiff_constants.h:336
FIFF_MNE_SOURCE_SPACE_NEAREST
#define FIFF_MNE_SOURCE_SPACE_NEAREST
Definition fiff_constants.h:346
fiff_coord_trans_old.h
FiffCoordTransOld class declaration.
fiff_sparse_matrix.h
FiffSparseMatrix class declaration.
k
int k
Definition fiff_tag.cpp:324
FIFF_NCHAN
#define FIFF_NCHAN
Definition fiff_file.h:453
FIFFB_BEM
#define FIFFB_BEM
Definition fiff_file.h:403
FIFF_COORD_TRANS
#define FIFF_COORD_TRANS
Definition fiff_file.h:475
FIFF_SUBJ_HIS_ID
#define FIFF_SUBJ_HIS_ID
Definition fiff_file.h:573
FIFFB_MRI
#define FIFFB_MRI
Definition fiff_file.h:390
fiff_types.h
Definitions for describing the objects in a FIFF file.
mne_forwardsolution.h
MNEForwardSolution class declaration, which provides the forward solution including the source space ...
mne_ctf_comp_data_set.h
MneCTFCompDataSet class declaration.
mne_source_space_old.h
MneSourceSpaceOld class declaration.
mne_named_matrix.h
MNE Named Matrix (MneNamedMatrix) class declaration.
mne_nearest.h
MneNearest class declaration.
compute_fwd.h
Compute Forward Setting class declaration.
FIFFLIB::FiffCoordTransOld
Coordinate transformation descriptor.
Definition fiff_coord_trans_old.h:81
FIFFLIB::FiffSparseMatrix
Data associated with MNE computations for each mneMeasDataSet.
Definition fiff_sparse_matrix.h:75
FIFFLIB::FiffSparseMatrix::inds
FIFFLIB::fiff_int_t * inds
Definition fiff_sparse_matrix.h:130
FIFFLIB::FiffSparseMatrix::ptrs
FIFFLIB::fiff_int_t * ptrs
Definition fiff_sparse_matrix.h:131
FIFFLIB::FiffSparseMatrix::data
FIFFLIB::fiff_float_t * data
Definition fiff_sparse_matrix.h:129
FIFFLIB::_fiffIdRec
ToDo Old implementation use new fiff_id.h instead.
Definition fiff_types_mne-c.h:214
FIFFLIB::FiffChInfo
Channel info descriptor.
Definition fiff_ch_info.h:75
FIFFLIB::FiffDirEntry
Directory entry description.
Definition fiff_dir_entry.h:68
FIFFLIB::FiffDirEntry::SPtr
QSharedPointer< FiffDirEntry > SPtr
Definition fiff_dir_entry.h:70
FIFFLIB::FiffDirNode::SPtr
QSharedPointer< FiffDirNode > SPtr
Definition fiff_dir_node.h:76
FIFFLIB::FiffId
Universially unique identifier.
Definition fiff_id.h:68
FIFFLIB::FiffId::isEmpty
bool isEmpty() const
Definition fiff_id.h:175
FIFFLIB::FiffId::new_file_id
static FiffId new_file_id()
Definition fiff_id.cpp:85
FIFFLIB::FiffInfo
FIFF measurement file information.
Definition fiff_info.h:85
FIFFLIB::FiffInfoBase::SPtr
QSharedPointer< FiffInfoBase > SPtr
Definition fiff_info_base.h:78
FIFFLIB::FiffStream
FIFF File I/O routines.
Definition fiff_stream.h:105
FIFFLIB::FiffStream::start_file
static FiffStream::SPtr start_file(QIODevice &p_IODevice)
Definition fiff_stream.cpp:1928
FIFFLIB::FiffStream::open_update
static FiffStream::SPtr open_update(QIODevice &p_IODevice)
Definition fiff_stream.cpp:1955
FIFFLIB::FiffStream::SPtr
QSharedPointer< FiffStream > SPtr
Definition fiff_stream.h:107
FIFFLIB::FiffTag::SPtr
QSharedPointer< FiffTag > SPtr
Definition fiff_tag.h:152
FWDLIB::ComputeFwd::updateHeadPos
void updateHeadPos(FIFFLIB::FiffCoordTransOld *transDevHeadOld)
Definition compute_fwd.cpp:2319
FWDLIB::ComputeFwd::m_eeg_forward
QSharedDataPointer< FIFFLIB::FiffNamedMatrix > m_eeg_forward
Definition compute_fwd.h:149
FWDLIB::ComputeFwd::sol
QSharedDataPointer< FIFFLIB::FiffNamedMatrix > sol
Definition compute_fwd.h:144
FWDLIB::ComputeFwd::storeFwd
void storeFwd(const QString &sSolName="default")
Definition compute_fwd.cpp:2443
FWDLIB::ComputeFwd::m_meg_forward_grad
QSharedDataPointer< FIFFLIB::FiffNamedMatrix > m_meg_forward_grad
Definition compute_fwd.h:148
FWDLIB::ComputeFwd::m_meg_forward
QSharedDataPointer< FIFFLIB::FiffNamedMatrix > m_meg_forward
Definition compute_fwd.h:147
FWDLIB::ComputeFwd::sol_grad
QSharedDataPointer< FIFFLIB::FiffNamedMatrix > sol_grad
Definition compute_fwd.h:145
FWDLIB::ComputeFwd::m_eeg_forward_grad
QSharedDataPointer< FIFFLIB::FiffNamedMatrix > m_eeg_forward_grad
Definition compute_fwd.h:150
FWDLIB::ComputeFwd::ComputeFwd
ComputeFwd(ComputeFwdSettings::SPtr pSettings)
Definition compute_fwd.cpp:1822
FWDLIB::ComputeFwdSettings::SPtr
QSharedPointer< ComputeFwdSettings > SPtr
Definition compute_fwd_settings.h:82
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
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::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::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
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
FWDLIB::FwdCoilSet::read_coil_defs
static FwdCoilSet * read_coil_defs(const QString &name)
Definition fwd_coil_set.cpp:416
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
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
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
MNELIB::MneNamedMatrix
Matrix specification with a channel list.
Definition mne_named_matrix.h:85
MNELIB::MneNearest
This is used in the patch definitions.
Definition mne_nearest.h:78
MNELIB::MneSourceSpaceOld
This defines a source space.
Definition mne_source_space_old.h:77
Generated on Sat Feb 14 2026 16:30:24 for MNE-CPP by doxygen 1.9.8